Animal hosts and drug treatments were performed as described previously (Wright et al., 2015). Briefly, adult male C57/BL6 mice (20–30 g) were housed under a 12 h light/dark cycle supplemented with food and water ad libitum. All animal experiments were approved by the Institute Ethics Committee of the Institute of Psychology, CAS.

During the methionine (Sigma, United States) experiments, the mice were injected subcutaneously with 1 g/kg (6.6 mmol/kg) L-methionine twice daily for 10 consecutive days using the CPP procedure (including pre-test, training, and post-test days). During training, the mice were injected with methionine 1 h before each behavioral experiment. The animals were divided into four groups: 1) saline + saline (SS), 2) MET + saline (MS), 3) MET + cocaine (MC), and 4) cocaine + saline (CS).

An unbiased conditioning protocol was used as described previously (Walters and Blendy, 2001). In brief, the mice were placed into a middle chamber and allowed to roam freely between the two side chambers for 15 min on day 1, and then arranged into control and experimental groups with equivalent pre-test scores. After experimental manipulation, the mice were paired up for 8 days. The saline group received saline in both sides of the chambers, while the drug groups were injected with cocaine (20 mg/kg, i.p., Qinghai Pharmaceutical Group Co. LTD., China) and saline on one side and saline only on the opposite side (Romieu et al., 2008). After each manipulation, the mice were confined to the corresponding conditioning chambers for 30 min before being returned to their cages. On the test day, the mice were allowed access to all chambers and a place preference score (CPP score) was assigned by subtracting the time spent in the drug-paired chamber from the time spent in the saline-paired chamber.

Two hours after the (day 10) CPP test, the animals were sacrificed and their mPFCs were surgically excised. The tissue was stored in liquid nitrogen immediately after extraction and then transferred to −80°C. Total RNA was extracted from the frozen tissues using the TIANamp DNA/RNA Isolation Kit (TIANGEN), including additional treatment with RNase-free DNase I (Ambion) for 30 min at 37°C to remove contaminating DNA.

A total of 1 µg of RNA was taken from a pool of two animals (500 ng each) with three biological replicates per experimental condition to generate an mRNA sequencing library as follows: briefly, after poly-A addition to the mRNA molecules, the mRNA was fragmented. The cleaved RNA fragments were then copied into first-strand cDNA before second-strand cDNA synthesis. A single ‘A’ base was added to the cDNA, to which an adaptor was ligated and enriched by PCR amplification. Single-strand DNA circles (ssDNA circles) were constructed and DNA nanoballs (DNBs) were generated. These DNBs were loaded into the patterned nanoarrays and pair-end reads of 100 bp were read on the BGISEQ-500 platform at Beijing Genomics Institute (BGI; Shenzhen, China).

Raw sequencing reads were filtered for clean reads using SOAPnuke (Mak et al., 2017; Chen et al., 2018) (v1.5.2, parameters −l 15, −q 0.2, −n 0.05) (https://github.com/BGI-flexlab/SOAPnuke). The HISAT pipeline (Kim et al., 2015) was applied to the align reads to the reference genome (mm10). StringTie (Pertea et al., 2015) was then used for transcriptional reconstruction. Subsequently, Cuffcompare (Trapnell et al., 2012) (Cufflinks tools) was utilized to compare the reconstructed transcripts and estimate the expression levels (FPKM, Fragments Per Kilobase of exon model per Million mapped fragments) of all detected isoforms. The coding potentials of the novel transcripts were predicted by CPC (Kong et al., 2007). The identification of DEGs (differentially expressed genes) was based on the negative binomial distribution of the DEseq2 package (Love et al., 2014). The cutoff of DEGs was p< 0.1.

A total of 1 µg of RNA from a pool of two animals (500 ng each) with three biological repeats per experimental condition was used to generate a miRNA sequencing library. RNA segments 18–30 nt in length were separated and recovered by PAGE. The subsequent steps were the same as those used to generate the mRNA sequencing library from the poly-A addition. The miRNA DNBs (DNA nanoball) were loaded into the patterned nanoarrays and 50-bp single-end reads were read through on the BGISEQ-500 platform at Beijing Genomics Institute (BGI; Shenzhen, China).

The raw sequencing reads were filtered with FASTQ to remove reads with low quality, adapter contamination, and lengths <16 nt (Cock et al., 2010). The clean reads were mapped to the reference genome using AASRA (Tang et al., 2017). Reads matching rRNAs and tRNAs were excluded. The remaining reads were aligned against the miRBase (v21) (Griffiths-Jones et al., 2008) using Bowtie allowing one mismatch. Unaligned sequences were pooled to identify novel miRNAs using miRDeep (v2.0.0.5) (An et al., 2013) with the default parameters. To investigate the expression profiles of miRNAs, the frequency of miRNA counts was normalized to TPM (tags per million) using the following formula: normalized expression = actual miRNA read count/total clean read count × 106.

Differentially expressed miRNAs between the paired groups were analyzed using DEGseq (Wang et al., 2010). The p-values calculated for each gene were adjusted to Q-values for multiple testing corrections using two alternative strategies (Hochberg, 1995; Storey and Tibshirani, 2003). To improve the overall accuracy of the DEGs results, a gene was defined as a DE-miRNA (differentially expressed miRNA) when the Q-values was ≤ 0.05 and the |log2fold_change| was >0.5.

The potential miRNA targets were identified using miRanda (v3.3a, parameters -en -20 -strict) (Enright et al., 2003; Betel et al., 2008) and TargetScan (v6.0, parameters-c 4) (Agarwal et al., 2015). Sequences predicted by both miRanda and TargetScan were considered miRNA targets. Furthermore, if there are more than 100 target genes predicted for a miRNA, the top 100 target genes were selected according to the scoring system. To reduce the false identification rate, the gene expression profiles of putative targets were required to be negatively correlated with miRNA profiles.

A STRING web-based tool database (https://string-db.org) was used to construct the PPI network of DEGs affected by cocaine and reversed by MET. The PPI network was constructed by mapping the genes to the STRING database with confidence scores of 0.4. Cytoscape (v3.9.0) was used to visualize and analyze the PPI networks (Shannon et al., 2003).

To identify intersecting clusters in the PPI networks, we used the Molecular Complex Detection plugin (MCODE) Cytoscape plugin with identification criteria including a degree cutoff of 2, a node score cutoff of 0.2, a k-core of 2, and a maximum depth of 100. Significant modules were those with MCODE scores of ≥4 and nodes ≥6, which constrained the cluster size for co-expressing networks (Bader and Hogue, 2003).

The Transcriptional Regulatory Relationships Unraveled by Sentence-based Text mining (TRRUST, http://www.grnpedia.org/trrust/) database (Han et al., 2018) was used to predict TFs regulating the DEGs based on existing literature. Subsequently, TransmiR 2.0 (Tong et al., 2019) was used to predict the potential target transcription factors of miRNA with validated information. Finally, a TF-miRNA-mRNA coregulation network was constructed to show the potential molecular mechanisms of cocaine initiation and extinction. Cytoscape was used to visualize the interactions in the TF-miRNA-mRNA coregulation network.

Functional annotation of GO and KEGG pathway enrichment was performed using the web-based DAVID v6.8 tool (Huang da et al., 2009). Gene and miRNA expression regulation and KEGG visualized with heatmap, volcano, bar plot, and bubble plot were plotted by https://www.bioinformatics.com.cn, an online platform for data analysis and visualization.

DEGs were validated by quantitative reverse-transcription PCR (RT-qPCR) using the Maxima SYBR Green qPCR Master Mix kit (Fermentas) according to the manufacturer’s instructions on an ABI Prism 7500 Sequence Detection System machine (Applied Biosystems Inc.). All real-time RT-qPCR data were normalized to SS expression level (see the additional data file for primer information).

Data are expressed as means ± SEM. Statistical analysis of qPCR data was performed using an unpaired t test with two tailed distributions using GraphPad Prism version 6.0. The results were considered statistically significant when p < 0.05.

We previously showed that L-methionine inhibits cocaine-induced behaviors (Zhang et al., 2021). Briefly, comparisons before and after drug treatment (pre-test and pro-test) for each group showed significant differences only in the cocaine-induced CPP group (p < 0.001)(Supplementary Figure S1). t-Test analysis of the pre-test data revealed that methionine administration significantly attenuated cocaine-CPP expression (CS vs. MC p < 0.05).

To determine the gene expression changes in the whole genome underlying the effects of cocaine CPP and methionine in response to cocaine CPP, we first performed mRNA-sequencing of the extracted PFCs of mice from each of the four treatment groups (SS: saline+saline, CS: cocaine+saline, MS: MET+saline, and MC: MET+cocaine). A total of approximately 50 million 100-bp, paired-end reads were obtained per sample after removing low-quality reads. Approximately 90% of the clean reads mapped to the genome and 55% of the clean reads uniquely mapped to the genome (Supplementary Table S2).

Analysis of the mRNA sequencing data in the CS group compared to the SS group revealed 655 DEGs (Figure 1A), including 333 upregulated genes and 322 downregulated genes. KEGG pathway analysis of the upregulated and downregulated DEGs identified certain pathways enriched in the upregulated genes that were functionally associated with drug addiction, including mmu04721: synaptic vesicle cycle (Luscher and Ungless, 2006) and mmu04080: neuroactive ligand-receptor interaction and AMPK signaling pathway (Gao et al., 2017) (Figure 1B). Moreover, the results of the GO analysis in the upregulated genes showed enrichment for GO terms that directly corresponded with the characteristics of cocaine addiction, including response to cocaine, locomotor behavior, and regulation of synaptic plasticity (Supplementary Table S3).

As shown in Supplementary Figure S1, MET did not alter behavior in the saline-treated group. To determine how MET treatment modified the mouse transcriptome, we simultaneously performed mRNA sequencing of the MS group. A total of 1531 DEGs were identified in the MS group compared to the SS group, including 204 upregulated DEGs and 674 downregulated DEGs (Figure 1C). KEGG pathway analysis of the upregulated and downregulated DEGs revealed that the addiction-related pathways; for example, the calcium signaling pathway, glutamatergic synapse, the MAPK signaling pathway, the cAMP signaling pathway, and dopaminergic synapse, were only enriched in the presence of the downregulated genes (Figure 1D). As the MET treatment did not induce CPP, we speculated that the downregulated genes enriched in drug addiction-related pathways did not induce CPP.

Next, we applied the same analyses to study the effect of L-methionine on cocaine CPP. We first identified induced by repeated cocaine treatment and inhibited by MET treatment. We performed a differential expression analysis between the MC and CS groups to obtain a list of DEGs affected by MET treatment in the context of cocaine treatment. A total of 644 DEGs were upregulated and 352 DEGs were downregulated in the MC-treated mice compared to the CS-treated mice (Figure 1E). The KEGG analysis plot showed that most of the pathways enriched with downregulated genes were associated with drug addiction, including those related to amphetamine addiction, glutamatergic synapse, nicotine addiction, GABAergic synapse, cocaine addiction, morphine addiction, etc. (Figure 1F). Combine with the results observed for MET treatment only, these findings strongly indicated that MET inhibited cocaine CPP by downregulating these addiction-related pathways.

We showed above that the MET treatment downregulated addiction-related pathways (Figures 1D, F). To further investigate the genes induced by cocaine and inhibited by MET, we plotted the log2fold_change of DEGs in the CS, MS, and MC groups in a heatmap. The resulting hierarchical clustering indicated that many genes exhibited fold changes in opposite directions in the MC mice compared to those in the CS and MS mice (Figures 2A–C). These genes were clustered into two groups as marked in black: genes inhibited by cocaine and induced by L-methionine (CS-down and MC-up) and the reverse: genes induced by cocaine and inhibited by L-methionine (CS-down and MC-up). The KEGG pathways of these two groups were analyzed separately to show that the pathways of the CS-down and MC-up genes were enriched for neurodegenerative diseases while the CS-up and MC-down genes were enriched for substance dependence (morphine addiction and nicotine addiction) and nervous function (glutamatergic synapse, cholinergic synapse, and GABAergic synapse) pathways (Figure 2D).

Posthoc analysis of the log2fold changes in genes related to substance dependence and nervous function in a less stringent manner showed a significant negative relationship between the CS and MC groups (p = 0.0000614 for morphine addiction, 0.0069 for nicotine addiction, 0.0000211 for glutamatergic synapse, 0.000000539 for GABAergic synapse, 0.000000305 for cholinergic synapse, and 0.00097 for circadian entrainment) (Figure 3A). This further supported the role of MET in reversing the addiction pathways modified by cocaine treatment by reversing the activation of genes in these pathways. We also plotted the top 10 most significant genes in the CS-down and MC-up group and the CS-down with MC-up group to identify key genes involved in drug addiction (Figures 3B, C). Among them, Fosb, Hrh3, Ins2, and Nfkbia were involved in drug addiction-associated pathways.

STRING provides original reliable protein data for consequent analysis. DEGs induced by cocaine and reversely regulated by L-methionine were mapped to generate PPI networks containing 1,040 nodes and 5,946 edges (Figure 4A). The top 10 hub genes of the total DEGs PPI network were Rps27a, Hras, Cdh1, Mapt, Acta2, Grin1, Isg15, Smad3, Nfkbia, and Trrap (Supplementary Table S4).

The overall PPI network of the DEGs was surveyed to identify functional modules in the network (Figure 4A; Table 1; Supplementary Table S5). Nine modules were identified in the PPI network, with MCODE scores ≥4 and nodes ≥6. To further investigate the functions of these modules, we performed pathway and GO analysis. In some modules, most genes were enriched in the same pathway categories. For example, most genes were enriched in ribosomes in module M1, while module M5 was enriched for genes related to the spliceosome. Meanwhile, module M4 was enriched for pathways associated with neurodegenerative diseases such as Parkinson’s disease, Alzheimer’s disease, and Huntington’s disease.

The pathways in module M8 were involved in substance dependence and nervous organismal systems. Morphine addiction, GABAergic synapse, and circadian entrainment were represented in this module. Module M7 was enriched in the cardiac muscle contraction and MAPK signaling pathways, which was a signal transduction pathway and reportedly involved in drug addiction and possibly required for the establishment of NAc amphetamine-produced CPP (Gerdjikov et al., 2004; El Rawas et al., 2020). Moreover, module M9 was enriched in the circadian rhythm pathway, which is reported as the common pathway underlying drug abuse (Li et al., 2008). However, though no KEGG pathways were assigned to Module M6, the GO terms for the biological processes of response to cAMP were enriched in this module. cAMP is a second messenger involved in the molecular network of drug addiction (Li et al., 2008).

The GO and KEGG analysis of the modules showed that the functions of modules M6–M9 were involved in drug addiction. Pathway networks were visualized using the ClueGO plugin (Figure 4B). Additionally, subnetworks of the DEGs enriched in the drug addiction-related pathways from modules M6–M9 were generated (Figure 4C). The subnetworks mainly comprised calcium channel (cacna1e, cacna1i, cacnb1, cacng4, and cacng8), glutamate receptor (grin3a, grin1, and grm1), G protein (gng10, gng11, and gng13), protein kinase C (prkcb and prkcg), and immediate early (Egr1 and Fosb) genes.

MicroRNAs (miRNAs) are key mediators of the silencing of post-transcriptional gene expression. We performed small RNA-seq to reveal the miRNome of the MET response to cocaine addiction. The miRNAs were to obtain approximately 12 million reads for each sample, of which ∼15% of reads were filtered out for low quality or short length (<15 nt) (Supplementary Table S6).

Compared to the SS group, 297 differentially expressed miRNAs (DE-miRNAs) were identified using DEseq2 in the cocaine-treated group (CS). Among these, 257 miRNAs were significantly highly induced by cocaine treatment, while 40 miRNAs showed significantly low expression (Figures 5A, D; Supplementary Table S7). To predict their target genes, we filtered the DEGs mediated by differentially expressed miRNA (DEGs-miRNA) according to the miRNA target prediction rules. Our findings showed enrichment of the cholinergic synapse and FoxO signaling pathways with upregulated DEGs-miRNA (Figure 5E). Meanwhile, compared to the SS group, the MS group showed 222 significantly high expressed miRNAs and 56 significantly lowly expressed miRNAs following CPP (Figures 5B, D; Supplementary Table S7). KEGG pathway analysis of miRNA targeting prediction showed that the axon guidance and calcium signaling pathway were enriched in downregulated predicted DEGs (Figure 5E).

The same analysis was applied to the MC group. We compared the MC group to the CS group to assess the effect of MET on response to cocaine. We identified 153 upregulated miRNAs and 76 downregulated miRNAs (Figures 5C–E). KEGG pathway analysis showed that most pathways enriched in downregulated-predicted DEGs were associated with drug addiction, including the GABAergic synapse, cholinergic synapse, glutamatergic synapse, circadian entrainment, axon guidance, and calcium signaling pathways, consistent with the pathways enriched by DEGs and indicating that MET acted on miRNAs to regulate gene expression.

To further investigate miRNA regulation in the MET treatment of cocaine CPP, we first plotted the log2fold_changes of DE-miRNAs in the CS, MS, and MC conditions as a heat map. As expected, a significant number of miRNAs exhibited fold changes in opposite directions in the MC mice compared to the CS mice, which was confirmed by hierarchical clustering marked in black (Figure 6A). We then plotted the fold-changes between the CS and MC groups and between the CS and MS groups, which showed a positive relationship (R = 0.42) between CS and MS and a negative relationship (R = −0.49) between CS and MC, consistent with the relationship of DEGs, indicating that MET inhibited cocaine CPP behavior by reversing the DEGs and DE-miRNAs regulated by cocaine (Figures 6B,C).

KEGG pathway analyses of DEG-miRNAs (DEGs regulated by DE-miRNAs) in the CS upregulation with MC downregulation subgroup and CS downregulation with MC upregulation subgroup were performed separately. The pathways associated with drug addiction were only enriched in the CS downregulation with MC upregulation subgroup, including the GABAergic synapse, glutamatergic synapse, circadian entrainment, axon guidance, and calcium signaling pathways (Figure 7A). These pathways were top-ranked pathways enriched with DEGs induced by cocaine and downregulated by L-methionine, indicating that the expression of these genes located in these pathways might be regulated by miRNAs in responses to cocaine and MET.

We then selected genes in the substance dependence, synapse, and other addiction-related pathways in addition to their source miRNAs to create a network. A total of 32 miRNAs targeted 22 genes, including those related to calcium channels (Cacna1c, Cacna1e, Cacna1g, and Cacng8), ephrin receptors (Ephb6 and Epha8), and ryanodine receptors (Ryr1 and Ryr2). Calcium channel genes were also present in the addiction-associated DEG PPI network shown in Figure 4C, indicating that the calcium channel gene network is a core gene network that plays an important role in the process of L-methionine counteracting cocaine CPP (Figure 7B).

MiRNAs and transcription factors (TFs) can work cooperatively with synergistic and antagonist actions as essential mediators of gene expression. Analysis of 148 DE-miRNAs mediated by cocaine and reversely regulated by L-methionine predicted potential significantly differential TFs in the TransmiR v2.0 database. Subsequently, the TRRUST database was used to search the target DEGs to construct transcriptional regulatory networks. Thus, a TF-miRNA-mRNA interaction network was constructed based on the TRRUST and TransmiR analysis and visualized using Cytoscape (Figure 8). The network consisted of 97 nodes and 173 edges. We identified 11 TFs that were significantly regulated by 24 DE-miRNAs. Among them, the TF of Nrf1 targeted by mmu-miR-15b, mmu-miR-486a, mmu-miR-532, mmu-miR-592, mmu-miR-615, and mmu-miR-7054 could regulate the expression of Grin1, which is a hub gene in the addicted associated DEG network interacting with calcium channel genes, glutamate receptor genes, and Egr1 shown in Figure 4C. The TFs of Crebbp targeted by mmu-miR-7054 and Nr3c1 targeted by mmu-miR-15b, mmu-miR-182/3, mmu-miR-195a, mmu-miR-2137, and mmu-miR-7054/8 could regulate Fosb and Nfkbia, the most significantly highly expressed genes affected by cocaine and reversely regulated by L-methionine. Fosb is one of the IEGs following cocaine experience and a transcription factor that regulates reward in the brain (Zhang et al., 2014), suggesting that Fosb may be one of the earliest affected genes by cocaine through transcription factors via miRNAs.

The top-ranked DEGs by fold-change expression were Grin1 (glutamate receptor ionotropic, NMDA 1) and Cacna1e (voltage-dependent R-type calcium channel subunit alpha-1E), which are involved in a variety of calcium-dependent processes including neurotransmitter release, and Fosb (FBJ osteosarcoma oncogene B) play a role in neurogenesis in the hippocampus and in learning and memory-related tasks. An RT-qPCR assay confirmed that the expression of these three genes was significantly altered by cocaine and reversed by MET, consistent with the results of RNA sequencing (Figure 9).