mRNA changes in nucleus accumbens related to methamphetamine addiction in mice

Methamphetamine (METH) is a highly addictive psychostimulant that elicits aberrant changes in the expression of microRNAs (miRNAs) and long non-coding RNAs (lncRNAs) in the nucleus accumbens of mice, indicating a potential role of METH in post-transcriptional regulations. To decipher the potential consequences of these post-transcriptional regulations in response to METH, we performed strand-specific RNA sequencing (ssRNA-Seq) to identify alterations in mRNA expression and their alternative splicing in the nucleus accumbens of mice following exposure to METH. METH-mediated changes in mRNAs were analyzed and correlated with previously reported changes in non-coding RNAs (miRNAs and lncRNAs) to determine the potential functions of these mRNA changes observed here and how non-coding RNAs are involved. A total of 2171 mRNAs were differentially expressed in response to METH with functions involved in synaptic plasticity, mitochondrial energy metabolism and immune response. 309 and 589 of these mRNAs are potential targets of miRNAs and lncRNAs respectively. In addition, METH treatment decreases mRNA alternative splicing, and there are 818 METH-specific events not observed in saline-treated mice. Our results suggest that METH-mediated addiction could be attributed by changes in miRNAs and lncRNAs and consequently, changes in mRNA alternative splicing and expression. In conclusion, our study reported a methamphetamine-modified nucleus accumbens transcriptome and provided non-coding RNA-mRNA interaction networks possibly involved in METH addiction.

Scientific RepoRts | 6:36993 | DOI: 10.1038/srep36993 accumbens (NAc) 12,16 , suggesting a potential role of lncRNAs in mediating drug addiction. We recently analysed METH-induced ncRNAs changes in the NAc of mice by using RNA sequencing technologies 17,18 and identified numerous METH-responsive miRNAs and lncRNAs. However, the potential post-transcriptional regulations of those ncRNAs on METH-induced mRNA changes have not been studied. Changes in expression of mRNA, specifically various mRNAs isoforms differentially spliced from a single gene can also be regulated post-transcriptionally by alternative splicing. Evidences both from in vitro non-nervous tissues and in vivo cocaine-treated mouse model suggested a regulatory role of alternative splicing for specific genes 19,20 . These changes in alternative splicing in response to cocaine expanded and diversified the proteome modulating drug addiction-related changes 21 . However, the contribution of alternative splicing to METH action is still unclear.
Although gene expression changes and their post-transcriptional regulation are implicated to play important roles in drug addiction, there is still insufficient data for METH-mediated changes in the NAc. Therefore, we sought to study genome-wide METH-induced mRNA expression in the NAc using ssRNA-Seq. The post-transcriptional regulation of mRNA by miRNAs, lncRNAs and alternative splicing were analyzed to reveal the potential regulatory consequences in response to METH.

Results
Overview of the ssRNA-Seq data. To examine transcriptional changes in response to METH, we used a METH administration regimen that has been previously demonstrated to produce robust locomotor sensitization. Animals were randomly allocated to the METH (METH1 and METH2) or saline-treated (saline 1 and saline 2) groups. Mice were injected daily with 2 mg/kg of METH or saline through the intraperitoneal route for 5 consecutive days followed by two injection-free days and a final challenge injection (Fig. 1A). These mice were sacrificed 24 h after the final injection and samples prepared for ssRNA-Seq to examine for transcriptome changes in the NAc of mice. An average of 50.5 million clean reads per sample passed filtering and we observed no difference in the alignment reads between groups (Fig. 1B). Assembling and annotation of the mapped reads have shown that an average of 15497 and 15396 mRNAs were expressed (reads > 1.0) in both samples from the saline and METH groups respectively. Given the nearly equivalent values for each mRNA within the samples of each group ( Fig. 2A,B), we used average reads of the shared mRNAs from each group to represent their expression levels in saline and METH subjects in the following analyses.

METH treatment induced aberrant expression of mRNAs in the NAc of mice.
Under a combination of statistical significance (P < 0.0001, FDR ≤ 0.00001), the differential expression analyses showed 2171 mRNAs exhibited differentially expressed level in the NAc of METH-treated mice, including 1059 up-regulated and 1112 down-regulated mRNAs (full list shown in Supplementary Table S1). Changes in mRNA expression identified differentially expressed genes (DEGs) and several of these DEGs have been implicated in the regulation of neuroplasticity and drug addiction, such as Arc (2.06-fold), Junb (1.68-fold), Fos (1.89-fold) and Ntrk1 (1.61-fold) (Supplementary Table S1).
In order to reveal the molecular and cellular functional characteristics affected by METH exposure, we subjected all of the DEGs to Gene Ontology (GO) analysis. As shown in Table 1, we observed significant effects of METH exposure on protein-coding genes (pcGenes) involved in synaptic transmission, synaptic plasticity, mitochondrial inner membrane functions, and oxidative phosphorylation. In addition, METH treatment also altered the expression of pcGenes that are responsible for viral infectious cycle.
Non-coding RNA-mRNA networks revealed biological functions in response to METH. METH have been shown to regulate expression levels of miRNAs and lncRNAs in the NAc of mice 17,18 , suggesting potential gene regulatory consequences following METH exposure. To investigate the post-transcriptional regulatory relationships between DEGs and ncRNAs, we performed regulatory analysis by aligning these DEGs to the putative targets for METH-responsive miRNAs and lncRNAs that we previously identified 17,18 . We found 309 DEGs (168 up-regulated and 141 down-regulated) that are putative targets for 32 METH-responsive miRNAs (31 down-regulated and one up-regulated, Table 2). These miRNAs-DEGs networks exhibited both positively and negatively correlated interactions (Fig. 2C). Positively correlated interactions refer to down-regulated miRNA and down-regulated DEGs, while negatively correlated interactions refer to down-regulated miRNA and up-regulated DEGs. Examples of negatively correlated interactions include miR-212-3p, miR-338-3p and miR-138-5p and their potential targets involved in neural functions, Arc, Fos and Ntrk1 respectively. These neural function-related DEGs were upregulated while their interacting miRNA were downregulated upon METH treatment 17 . We also found 586 DEGs (109 up-regulated genes and 477 down-regulated genes) associated with genes for 3988 differentially expressed lncRNAs in response to METH. The lncRNAs-DEGs networks also exhibited both negatively and positively correlated interactions (Fig. 2D).
We next carried out pathway analysis for those DEGs that are associated with changes in miRNA or lncRNA to find potential functions of these associations. This pathway analysis was done using Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.kegg.jp/), a bioinformatics tool based on external pathway data sources. All significantly enriched pathways (P < 0.05and EF > 2.0) were shown in Table 3. Expectedly, both miRNA targets and lncRNA associated genes showed significant enrichment in pathways that were verified to be responsible for drug addiction, such as amphetamine addiction, glutamatergic synapse, dopaminergic synapse, and MAPK signalling pathway (Table 3). Presented here are pathways of amphetamine addiction (Fig. 3) and glutamatergic synapse ( Supplementary Fig. S1), regulated by differentially expressed miRNAs and lncRNAs associated with changes in DEGs in the NAc of mice following METH treatment. Several significantly altered DEGs involved in vibrio cholerae infection and mTOR signalling pathway were also predicted to be regulated by both METH-induced changes in miRNAs and lncRNAs ( include metabolic pathways, such as fatty acid elongation, arginine and proline metabolism, oxidative phosphorylation, biosynthesis of secondary metabolites. Several other signalling pathways such as ErbB signaling pathway, phosphatidylinositol signaling system and wnt signaling pathway were predicted to be modulated by METH-dependent lncRNAs-DEGs networks. Interestingly, we also observed DEGs that function in ribosome, spliceosome and protein export are potentially modulated by miRNAs, suggesting an alternative splicing regulation following METH exposure.

METH-induced variation on alternative splicing.
To determine the splicing differences in response to METH, we analyzed alternative splicing events (ASEs) on all clean data of each group. The ASEs analysis included exon skipping, intron retention and alternative 3′ and 5′ splicing sites (A3′ SS and A5′ SS). The numbers of ASEs corresponding to each type in both saline and METH-treated groups were identified and calculated. As observed previously, the most frequent event type in both groups was exon skipping, followed by A3′ SS, A5′ SS and intron retention 22 . However, the numbers of ASEs in METH-treated group were about 10-40% less than those in saline-treated group (Fig. 4A-D). As compared to all other ASEs, the intron retention events in METH-treated group showed the biggest decrease of about 40%. Further analysis considering group-specific alternative splicing The mice were handled daily for 1 week to adapt to the experimenter before treatment. Pre-test (day 1-2): mice were given once daily injections of saline for two consecutive days. Development (day 3-7): After pre-test, mice were divided into two groups randomly, and were given once daily injections of METH (2 mg/kg) or saline for five consecutive days. Transfer (day 8-9): two groups of mice were experienced two injection-free days. Expression (day 10): mice were given a challenge injection of either 2 mg/kg dose of METH or saline. Horizontal locomotor activity was performed 60 minutes before and after the injections on all drug treatment days. (B) Statistical alignment of the sequencing data from the saline and METH groups. Unique mach represents the mapped reads that aligned to only one position in the mouse genome. Multi-position mach represents the mapped reads that aligned to more than one positions in the mouse genome. Total unmapped indicates reads that did not match to the mouse genome. events (gsASEs), defined as a splicing event identified only in one group but not in the other group, showed a total of 1450 events that were saline-specific splicing variants while a total of 818 events were METH-specific (Fig. 4E).
Validation of changes in gene expression using qPCR. To validate changes in gene expression from the ssRNA-Seq analysis, confirmatory qPCR was conducted on 12 DEGs. Genes were selected based on fold change, functional enrichment, and potential role in drug addiction. Seven of these genes (Junb, Ntrk1 Fos, Dicer1, Rps23, Avpr1a and Gbp4) were significantly changed in METH-treated mice, consistent with the ssRNA-Seq data. The other five genes examined (Ptprv, Mrpl52, Rgs11, Rreb1, Ligp1) also showed increasing or decreasing trend, that are consistent with our ssRNA-Seq data. Thus, these data verified the consistency between the two independent techniques used for the study. The differential expression of genes known to be involved in drug addiction (Junb, Ntrk1, Ptprv, Fos) and miRNA regulation (Dicer1) are presented in Fig. 5A-E with their paired ssRNA-Seq results. The other seven genes are presented in the Supplementary Figure S2. The comparison of the fold change obtained using qPCR and ssRNA-Seq indicated a strong correlation between both techniques (Fig, 5F, r = 0.896, P < 0.0001). In summary, these confirmations by qPCR indicated a strong reproducibility of changes in genes expression by using another independent technique.

Discussion
High-throughput genome-wide approaches have been widely utilized in animal models to provide an overview of global gene regulation pattern, specifically in relation to various diseases. In this study using ssRNA-Seq, we identified 2171 DEGs in the NAc of METH-treated mice. Many of these DEGs are known to be critical for brain functions. GO analysis showed significant enrichment of METH-responsive genes involved in neuronal plasticity, mitochondrial energy metabolism and immune response, based on GO-defined changes in molecular and The number of uniquely mapped reads per kilobase per million mappable reads. Interaction networks formed by miRNA-mRNAs and lncRNAs-mRNAs were shown by the fold change of miRNAs (C, x-axis) and lncRNAs (D, x-axis) identified in NAc of mice following METH exposure 17,18 and their potential target genes with significant changes (y-axis) in response to METH (current study). cellular functions and biological processes. Previous studies demonstrated that exposure to abused drugs results in neural dysfunction that is critical for modulating addictive behaviors [23][24][25] . Thus, it is not surprising that GO analysis show relatively higher number of genes involved in neural functions.
Given that METH is a powerful psychostimulant that can produce robust rewarding effects through modulating the dopaminergic synapse, we expectedly found several DEGs (such as Comt, Akt, Gsk3a/b, PP2A, Slc family and PLC) that are presented in the dopaminergic synapse pathway (KEGG pathway database). These genes were either related to dopamine re-uptake (Comt, Slc family) or dopamine-regulated downstream signals (Akt, Gsk3a/b, PP2A, and PLC). Interestingly, we found some overlaps in METH-and cocaine-responsive genes. 136 DEGs in NAc of METH-treated mice were also changed in the NAc of mice treated with cocaine (Supplementary Figure S3) 12 . Although cocaine can also produce rewarding effects through interaction with the mesolimbic dopaminergic system, we did not find any dopamine-related changes among these overlaps. Therefore, the observed gene changes from the treatment with these two different drugs suggest the presence of common and differential mechanisms in mediating METH and cocaine addiction. Further studies are needed to investigate these common and differential pathways and mechanisms in drug addiction.
METH treatment also altered the expression of pcGenes that are involved in other non-neural functions. Notably, we found a large number of mitochondrial energy metabolism changes in response to METH, suggesting a neurotoxic response may have been induced in the brain by the METH regimen used in the current study. This result is consistent with the previous studies. Based on extant literatures, either a single administration or repeated injections of METH could causes energetic metabolism impairment in the amygdala, prefrontal cortex, hippocampus and striatum of rats accompanied with a significant behavioral sensitization 26,27 . Mitochondrial energy metabolism provides energy source responsible for brain functions and mitochondrial energy metabolism disruptions are known to play important roles in psychiatric disorders such as schizophrenia, bipolar disorder, and major depression disorders 28,29 . Thus, the effects of METH exposure on mitochondrial functions suggest that the molecular adaptations of cells responsible for a large extent of brain energy metabolism may play a role in METH addiction. In addition, we found significant changes in ionotropic glutamate receptor genes, such as Grin2a, Grin2b and Grin2d. Considering that mitochondria in nerve terminals are capable of maintaining cellular calcium (Ca 2+ ) homeostasis 30,31 and neuronal mitochondrial Ca 2+ processing supports Ca 2+ influx via ionotropic glutamate receptors 32 , these ionotropic glutamate receptor alterations may possibly involved in the dysregulation of mitochondrial calcium load contributing to METH addiction. Exposure to METH also causes inflammatory response that may play a role in METH-induced neuronal injury 33 . Furthermore, there is increasing evidence showing METH (both at low and high doses) induces a heightened inflammatory environment in the brain, that affects the regulation of cellular immune responses 34,35 . In line with these literatures, we found that METH treatment altered the expression of pcGenes that are involved in immune response. Therefore, how METH having divergent impacts on other responses in addition to neural functions warrants further investigations.
In addition to the mRNA expression patterns, the underlying molecular regulatory networks are also important to understand how various mRNA changes in response to METH can regulate potential functions. Previously, we have identified a subset of METH-responsive miRNAs and lncRNAs in the NAc of mice and predicted their target genes using bioinformatics tools. Thus, in order to determine how non-coding RNAs are involved in METH-induced various mRNA changes, we performed a comparison between the DEGs and the potential target genes of miRNAs and lncRNAs. We identified numerous potential METH-responsive miRNA-mRNA and  lncRNA-mRNA networks, which suggested changes in the complex gene regulatory networks upon METH exposure. Since miRNA is known to suppress mRNA expression and considering the large down-regulation pattern of miRNAs 17 , our data revealed a larger negatively correlated networks (down-regulated miRNAs with up-regulated mRNAs) and a smaller positively correlated networks (down-regulated miRNAs with down-regulated mRNAs). This suggested that the up-regulation of these relevant mRNAs maybe due to the down-regulated miRNAs in response to METH. The transcription data also revealed coordinative changes between mRNAs and miRNAs while the putative target genes of METH-regulated miRNAs did not show global changes. These data indicated  that the effects of METH exposure on genes transcription were not only regulated by miRNAs but may also be resulting from direct effects of METH or indirect effects of other regulators. Besides miRNA-mRNA networks, we also identified complex interaction networks formed by lncRNAs-mRNAs. LncRNAs are demonstrated to be transcription and epigenetic modulators, and are involved in various mechanisms to regulate genes expression. Unlike miRNAs, lncRNAs are able to increase or decrease the expression of pcGenes. Thus, we observed a more evenly distribution of both negatively and positively correlated interactions with mRNAs. This study identifying potential METH-responsive lncRNA-mRNA interaction networks, suggests a novel regulatory role of lncRNAs following METH exposure.
The altered ncRNAs-mRNAs networks underscored the biological importance of induction of ncRNAs targeting in response to transcriptional regulation induced by METH. Therefore, we performed KEGG pathway analysis on METH-responsive gene lists for miRNAs and lncRNAs, and observed a significant enrichment of miRNA-mRNA and lncRNA-mRNA on pathways involved in drug addiction (Table 3). For example, predicted METH-related targets include a range of amphetamine addiction-and synaptic plasticity-related receptors and downstream signals, such as glutamatergic receptors (Grin2a, Grin2b, Gria1, Gria3, Gria4), Gnas, Jun and Fos ( Fig. 3 and Supplementary Figure S1). In addition, we also identified that up-regulated genes implicated in the regulation of neuroplasticity such as Arc and Ntrk1 36,37 were potential targets for miR-212-3p and miR-138-5p respectively. Both miR-212-3p and miR-138-5p were significantly down-regulated in response to METH and were previously demonstrated as important regulators of cocaine addiction and neuronal plasticity 7,38 . These results suggested that ncRNAs are involved in the modulation of METH addiction. Moreover, we also observed enrichment of miRNA-mRNA and lncRNA-mRNA on immune response related pathways, such as vibrio cholerae infection, suggesting an involvement of ncRNAs in the METH-induced immune dysfunction. These results indicated the effects of METH on neural and other non-neural systems were, at least partially, dependent on changes of miRNAs and lncRNAs.
Alternative splicing is known as another kind of the post-transcriptional mechanisms for generating diverse gene products. Abnormal changes in splicing can result in disease consequences 39 . Interestingly, we found that Dicer1, an endonuclease responsible for the generation of miRNAs from their precursors 40 , was significantly down-regulated in response to METH. This finding implicated dysregulation of splicing of miRNAs in the NAc of METH-treated mice 17 . We also identified contribution of changes in alternative splicing induced in NAc of METH-treated mice (Fig. 4). Although there are a few recent studies on splicing regulation of genes transcription in addiction models 20,41,42 , our data is the first comprehensive analysis of METH-induced splicing regulation in NAc of mice. These data underscored the importance of discrimination of isoform differences on total gene transcription and also expanded the degree to which METH modified the transcriptome in NAc. In addition, the impact of spliceosome on pre-mRNA splicing requires interactions of pre-mRNA with splicing factors. Pathway analysis of DEGs that were putative target genes of differentially expressed miRNA and lncRNAs was carried out to assess their putative biological functions ( Table 2). As examples, amphetamine addiction pathway composed of DEGs and their corresponding miRNAs and lncRNAs was shown. Genes with red font suggested an up-regulation while genes with blue font indicated a down-regulation in response to METH. Our previous work revealed that lncRNAs such as Neat1, Neat2 and Miat, function as cofactors for pre-mRNA splicing through interaction with splicing factors [43][44][45] , were significantly down-regulated in response to METH. Moreover, we identified that miRNAs might modulate DEGs that function in ribosome and spliceosome, critical for pre-mRNA splicing. These results indicated that alternative splicing regulation of genes transcription in drug addiction may require ncRNAs-related nuclear modification.
Several limitations should be considered in the present study using bioinformatics tools to relate our DEG data to the ncRNA data. It is important to note that changes in non-coding RNA and mRNA do not always directly translate to the expression and functions of proteins. Although we currently observed a large number of mRNAs variants, but variations of protein products and changes of the localization and ultimate function of these variant protein products requires more intensive investigations. Subsequent verifications are then required to determine if manipulating the expression of the ncRNAs can recapitulate the effects of METH and/or alter METH-induced DEGs. Another limitation is the currently available methods and technologies to efficiently isolate pure core and shell region as well as different neuron subpopulations within each region. The core and shell subregions of the NAc are components of distinct microcircuits and exhibit different drug-induced changes 46 . However, the core and shell regions of the NAc are not distinguished in the present study. Nonetheless, this is the first description of the global mRNAs expression profiling and overview of post-transcriptional regulation in the NAc by miR-NAs and lncRNAs in the context of METH-induced behavioral sensitization. Future studies include identifying the precise regulatory mechanisms of each of these ncRNAs-DEGs networks, and also the distinct changes and functions between the different neuron subpopulations in different regions of the NAc. This may further strengthen our understanding of the relationships between the ncRNAs-DEGs networks and METH addiction.
In summary, we reported a METH-modified transcriptome in the NAc and provided an overview of post-transcriptional regulation by miRNAs, lncRNAs and alternative splicing upon METH exposure. Functional analysis of all DEGs and the DEGs correlated with changes in miRNAs and lncRNAs revealed that the effects of METH on neural and other different systems were, at least in part, dependent on changes of miRNAs and lncRNAs. A large number of METH-mediated changes in mRNA variants and the identification of alternative splicing changes expanded the current view on how METH can modify the transcriptome. Our results provided important information concerning the basic mechanisms and potential consequences of post-transcriptional regulation in response to METH in vivo.

Methods
Animals. Adult wild-type C57BL/6 mice (7-8 weeks old, male, 20-25 g), purchased from Beijing Vital River Laboratory Animal Technology Co. Ltd, were housed in cages (4 per cage) within regulated environment (23 ± 1 °C, 50 ± 5% humidity) with 12-h light/dark cycle (lights on from 7:00 am to 7:00 pm). The mice had ad libitum access to food and water and were handled daily for 1 week for adaptation to the experimenter before treatment. The experiments involving animals were approved and all treatment procedures were carried out in accordance with guidelines by the Institutional Animal Care and Use Committee of Xi'an Jiaotong University.
METH treatment and NAc dissection. METH hydrochloride used in this study was purchased from the National Institute for the Control of Pharmaceutical and Biological Products (Beijing, PR China) and was dissolved in sterile physiological saline. The METH injection paradigm used in this study has been demonstrated in previous studies to produce robust locomotor sensitization 47,48 . Mice were given once-daily injections of saline for two consecutive days (day 1-2), after which they were randomly grouped into two. The groups of mice were then received once-daily intraperitoneal (i.p) injections of METH (METH group, 2 mg/kg/injection, n = 16) dissolved in sterile physiological saline or saline (saline group, n = 16) for five consecutive days (day 3-7). After five injected days, the mice were housed in cages for two injection-free days (day 8-9). On day 10, the mice were received a challenge injection of either METH or saline. Before the beginning of the experiments, the mice were brought into the experiment room for 1 h to adapt to the experimental environment. Meanwhile, on all drug treatment days, horizontal locomotor activity was performed in the open field apparatus (43 cm × 43 cm × 43 cm) 1 h before and after the injections. Previous studies including our own works have shown changes in gene expression at 24 hr in drug-induced models to allow for evaluation of the steady state response 12,17,18,49,50 . For comparison and consistency, the mice used in the present study were also sacrificed 24 hr after the final injection. The brains were removed quickly and the NAc (+ 1.70 mm from Bregma) harvested according to the brain atlas of Paxinos and Franklin 51 were immediately frozen in liquid nitrogen until use. Sixteen mice per group were used to generate two independent samples of eight mice each and NAc lysates of eight mice from each group were pooled as one sample for total RNA isolation. RNA extraction. Total RNA was extracted using TRIzol following instructions of the manufacturer (Invitrogen, USA). Given two independent samples were prepared for each group, the RNA samples of each group were referred as saline1, saline2 and METH1, METH2. Samples were eluted in RNase-free H 2 O, and quantified using Agilent 2100 BioAnalyzer (Agilent Technologies, USA). All samples were sufficient to construct the strand-specific cDNA libraries (> 180 ng/μ l total RNA) and a RIN > 8.0 were accepted for RNA sequencing analysis.
ssRNA-Seq. ssRNA-Seq were carried out by the Beijing Genomics Institute Shenzhen (BGI Shenzhen, China). Total RNA (5 μ g) from each sample was used for construction of strand-specific cDNA libraries, which were prepared as previously described 18 . The cDNA libraries were evaluated using the Agilent 2100 BioAnalyzer and then sequenced on an Illumina HiSeq 2000. The sequencing reads alignment and assembling were done using TopHat v2.0.4 52 and Cufflinks v2.0.0 53 packages after filtering off dirty reads. The mRNA transcripts were annotated by perfect sequence mapping to the databases, including KEGG Orthology, non-redundant protein database(ftp://ftp.ncbi.nih.gov/blast/db/FASTA/),COG(http://www.ncbi.nlm.nih.gov/COG/) and UniProtKB/ Swiss-Prot 54 . For relative expression analysis, the number of uniquely mapped reads per kilobase per million mappable reads (RPKM) for each pcGene was calculated 55 . RPKM across samples of each group were then averaged because the values for each transcript were very consistent. The fold change in expression between METHand saline-treated mice was obtained by dividing RPKM METH by RPKM saline . The alternative splicing analyses was performed with SOAPsplice 56 and the alternative splicing events of each sample were counted in each group. Only those events appeared in both samples of each group were included for further analysis.
Functional enrichment. The differentially expressed miRNAs and lncRNAs data used for all comparison analysis in this study were previously identified by our published works 17,18 . The potential target genes of these miRNAs and lncRNAs were mapped to the differentially regulated pcGenes identified here. For functional annotation, we used Blast2GO 57 to analyze the GO and KEGG pathway enrichment for the DEGs and calculated the enrichment factor (EF) for each term or pathway. The P values were calculated via hypergeometric tests and go with a correction. Only those GO terms or KEGG pathways with enrichment of corrected P < 0.05 and EF > 2.0 were included. qPCR assay. cDNAs were synthesized from total RNA of the NAc of individual mice using Thermo Scientific RevertAid First Strand cDNA Synthesis Kit (Thermo Scientific, USA). Based on the manufacturer suggested parameters (25 °C 5 min, 42 °C 60 min, and 70 °C 5 min), 500 ng of total RNA from each sample was utilized for each reaction. qPCR was performed on Bio-Rad iQ5 system (Bio-Rad, USA) using SYBR Premix Ex Taq II (TaKaRa Biotechnology, Japan) under the following conditions: 95 °C for 30 sec; 95 °C for 10 sec and 60 °C for 1 min, repeated for 40 cycles. Gapdh was used as an endogenous control for the qPCR, and the relative expression levels were determined by the 2 −△△Ct method 58 . Primer pairs were chosen to keep the melt temperature (Tm) between 56 °C and 62 °C and products in the100-200bp size range (Supplementary Table S2).
Statistical analyses. The variances in sequencing assays were directly estimated by the Poisson distribution 59 and a false discovery rate (FDR) 60 for each pcGene was calculated to compensate for false-positive findings at each significance threshold. The following criteria: P < 0.0001, FDR ≤ 0.00001 was used to set up differentially regulated pcGenes. For the significance of qPCR assay, we performed an independent-sample t-test and data with significance P < 0.05 were considered to be differentially expressed (SPSS v17.0, SPSS Inc., USA). Pearson's coefficient analysis was also performed with SPSS (version 17.0, USA).