- Research article
- Open Access
Methamphetamine induces alterations in the long non-coding RNAs expression profile in the nucleus accumbens of the mouse
BMC Neurosciencevolume 16, Article number: 18 (2015)
Repeated exposure to addictive drugs elicits long-lasting cellular and molecular changes. It has been reported that the aberrant expression of long non-coding RNAs (lncRNAs) is involved in cocaine and heroin addiction, yet the expression profile of lncRNAs and their potential effects on methamphetamine (METH)-induced locomotor sensitization are largely unknown.
Using high-throughput strand-specific complementary DNA sequencing technology (ssRNA-seq), here we examined the alterations in the lncRNAs expression profile in the nucleus accumbens (NAc) of METH-sensitized mice. We found that the expression levels of 6246 known lncRNAs (6215 down-regulated, 31 up-regulated) and 8442 novel lncRNA candidates (8408 down-regulated, 34 up-regulated) were significantly altered in the METH-sensitized mice. Based on characterizations of the genomic contexts of the lncRNAs, we further showed that there were 5139 differentially expressed lncRNAs acted via cis mechanisms, including sense intronic (4295 down-regulated and one up-regulated), overlapping (25 down-regulated and one up-regulated), natural antisense transcripts (NATs, 148 down-regulated and eight up-regulated), long intergenic non-coding RNAs (lincRNAs, 582 down-regulated and five up-regulated), and bidirectional (72 down-regulated and two up-regulated). Moreover, using the program RNAplex, we identified 3994 differentially expressed lncRNAs acted via trans mechanisms. Gene ontology (GO) and KEGG pathway enrichment analyses revealed that the predicted cis- and trans- associated genes were significantly enriched during neuronal development, neuronal plasticity, learning and memory, and reward and addiction.
Taken together, our results suggest that METH can elicit global changes in lncRNA expressions in the NAc of sensitized mice that might be involved in METH-induced locomotor sensitization and addiction.
It has been reported that repeated exposure to drugs of abuse results in long-lasting behavioural changes such as locomotor sensitization that are thought to be due to structural and functional changes in associated brain regions, particularly the NAc [1,2]. Previous studies have indicated that addictive drugs induce persistent and dynamic cellular and molecular modifications accompanied by distinct processes of drug addiction [3-6]. The complexity of drug-induced stable changes suggests that synchronized programs of gene regulation might be executed during drug addiction ; however, the precise mechanism remains unclear.
Previous research has shown that addictive drugs induce aberrant expression of non-coding RNAs (ncRNAs) [8,9], the function of which are thought to be among the most important mechanisms underlying gene regulation [10,11]. For example, increasing evidence has demonstrated that microRNAs play an important role in modulating the potency of addictive drugs by mediating the expressions of target genes [7,12,13]. However, the potential effects of lncRNA on drug addiction remain largely unknown.
LncRNAs are defined as transcripts longer than 200 nt that lack the ability to encode protein products . As transcriptional modulators and epigenetic regulators, lncRNAs have been found to regulate the expressions of proximal and distal protein-coding genes through cis- and trans-acting mechanisms . Emerging evidence has implicated lncRNAs in neuroplasticity, brain development, neurodegenerative, and neuropsychiatric disorders [15-18], together, this evidence suggests a meaningful role of lncRNAs in brain diseases including drug addiction. Using microarrays, recent studies have revealed that cocaine and heroin induce widespread alterations of lncRNAs in the NAc of cocaine-conditioned mice and heroin addicts [9,19],which suggests that lncRNAs might play an important role in the regulation of drug addiction. However, all of these studies targeted small numbers of candidate lncRNAs and therefore could not identify unknown lncRNAs and did not provide a complete spectrum of drug-induced changes in lncRNA levels.
To investigate the expression profiles of lncRNAs and their potential effects on METH-induced locomotor sensitization in the current study, we examined the alterations in lncRNAs expression profiles in the NAc of METH-sensitized mice via the transcriptomics-based approach, ssRNA-seq. We found that METH elicited global changes in lncRNAs expression in the NAc of mice and that predicted cis- and trans-associated genes were significantly enriched during neuronal development, neuronal plasticity, learning and memory, and reward and addiction. Our results suggest that lncRNAs might be involved in the regulation of expression of associated genes and thus contribute to METH-induced locomotor sensitization and addiction.
Complementary DNA samples generated from RNA that was extracted from the NAc lysates of saline and METH-treated mice were measured with ssRNA-seq. A total of 49.62 million and 50.33 million clean reads were obtained from the saline and METH groups of mice, respectively. The clean reads of each group were then separately aligned to the mouse genome (UCSC mm9)  and 84.02% (saline) and 80.82% (METH) reads were mapped to the reference genome, which included 28.65 million (saline) and 27.08 million (METH) perfectly matched reads (Table 1). Additionally, among the total mapped reads, there were 36.76 million (74.09%) and 33.14 million (65.84%) uniquely matched reads that were obtained from the saline and METH groups of mice, respectively (Table 1). All of the mapped reads were then assembled and annotated. Consequently, 25677 and 23579 lncRNAs were obtained from the saline and METH groups of mice, respectively, by alignment to the database of non-coding RNAs (NONCODE v3.0)  (Table 2). The analysis of the relative expression levels (RPKM) of these known lncRNAs revealed that the vast majority of the lncRNAs (71.31% and 78.35% of the total known lncRNA transcripts in saline and METH groups of mice respectively) were expressed at the level of RPKM < 5 (Table 2). Moreover, we also identified 17860 and 15965 novel lncRNA candidates in the saline and METH groups, that could not matched to any sequences that correspond to known lncRNAs or protein-coding transcripts.
METH- induced aberrant expressions of lncRNAs in the NAc of mice
To further identify the differentially expressed lncRNAs in the NAc of the METH-sensitized mice, the RPKM ratios of the lncRNAs in each group were subjected to a log-2 transform to produce fold changes and threshold based on a combination of statistical significance (P < 0.001, FDR ≤ 0.0001) and the absolute value of the fold change (≥1.25) was set. We found that 6246 known lncRNAs exhibited significantly altered expressions in the NAc of the METH-sensitized mice that included 31 up-regulations and 6215 down-regulations (Additional file 1). A volcano plot illustrated the variance in the lncRNAs numbers at different P-values and fold changes (Figure 1).
Additionally, analysis of the novel lncRNA candidates revealed that 8442 novel lncRNA candidates exhibited significantly different expressions (P < 0.001, FDR ≤ 0.0001 and absolute value of the fold change ≥ 1.25) that included 34 increases and 8408 decreases (data were not shown).
Validation of lncRNAs expression by quantitative real-time PCR
To validate the METH-induced changes in lncRNA expression that were detected by ssRNA-seq, 12 differentially expressed lncRNAs were randomly selected, and their expressions were then examined with quantitative real-time PCR (qPCR). As shown in Figure 2, 10 of the selected lncRNAs were significantly changed (eight down-regulated and two up-regulated) in METH-sensitized mice as detected by both ssRNA-seq (Figure 2A, #P < 0.001, FDR ≤ 0.0001 and an absolute value of the fold change ≥ 1.25 compared to the saline group of mice) and qPCR (Figure 2B, *P < 0.05 compared to the saline group of mice, n = 11-15 per group). Although the expression of AK036791 and AK080587 detected by qPCR were not significantly regulated, they showed similar expression trends in qPCR and ssRNA-seq. Moreover, a strong agreement across the two methods was observed in that the results of the qPCR were similar to those obtained from the ssRNA-seq analyses (Figure 2C, r = 0.89, P < 0.05). These data indicated the good reproducibility of the observed expression changes in the lncRNAs based on an independent method.
Genomic characterization of differentially expressed known lncRNAs
To predict the potential role of lncRNAs in the regulation of the expressions of protein-coding genes, we next investigated the genomic context of the differentially expressed known lncRNAs. LncRNAs can be classified into intergenic lncRNAs, sense-overlap lncRNAs and antisense-overlap lncRNAs based on their locations relative to protein-coding genes. Here, we identified 5141 sense-overlap (5128 down-regulated and 13 up-regulated), 172 antisense-overlap (164 down-regulated and eight up-regulated) and 933 intergenic (923 down-regulated and 10 up-regulated) lncRNAs among the differentially expressed known lncRNAs that we detected (Figure 3).
LncRNAs in cis regulation
LncRNAs can regulate the expressions of genes that are located on that same chromosome, and such regulation is called cis regulation . Based on their genomic localization relative to nearby protein-coding genes, lncRNAs can be further classified as sense intronic, overlapping, NAT, lincRNA and bidirectional, and these classes of lncRNAs have been reported to regulate their protein-coding host genes in cis manners [22,23]. In the present study, we subjected the differentially expressed known lncRNAs to cis analysis, and we found that 82% (5139 of 6246) of the differentially expressed known lncRNAs could act in a cis manner, including 4296 sense intronic lncRNAs, 26 overlapping lncRNAs, 156 NATs, 587 lincRNAs and 74 bidirectional lncRNAs (Figure 4A).
Sense intronic lncRNAs
Sense intronic lncRNAs originate from long introns that are transcribed from the same strand as the associated protein-coding genes. The sense intronic lncRNAs are biologically significant because they have been found to be both co-expressed with their host protein-coding gene and independently expressed, particularly in the mouse brain . In the present study, we identified 4296 known lncRNAs that were located in the introns of protein-coding genes that include 4295 down-regulated and one up-regulated lncRNAs (Figure 4A, Additional file 2).
Overlapping lncRNAs are lncRNAs contain a protein-coding gene and are transcribed in the same direction as that gene. This type of lncRNA can regulate downstream transcription by opening the chromatin structure, depositing histone marks , and cis-acting promoter competition . Here, we identified 26 significantly changed lncRNAs that overlapped with protein-coding genes, including 25 decreased and one increased lncRNAs (Figure 4A, Additional file 3).
NATs are lncRNAs that are transcribed from the antisense strand of a gene locus, and are overlapping with the RNA that transcribed from the sense strand. LncRNAs of this type have been discovered to be widespread in the mammalian genome and work through multiple mechanisms to regulate the expressions of their sense partners [27,28]. In the present study, we identified 156 differentially expressed NATs that included 148 that were down-regulated and eight that were up-regulated (Figure 4A, Additional file 4). For example, the potassium voltage-gated channel, subfamily Q, member 1, opposite strand/antisense transcript 1 (Kcnq1ot1) [Genebank: NR_001461], and the zinc finger homeobox 2, antisense (Zfhx2as) [Genebank: AK032589] were found to be significantly down-expressed in the METH-treated mice and have previously been found to modulate the expression of their sense partners [29,30].
LincRNAs are found more than 10 kb away from any nearby protein-coding locus [31,32]. A possible working model of the role of lincRNAs in gene regulation involves their actions as enhancers that activate transcriptional promotion and chromatin looping [33,34]. We identified 587 lincRNAs that were significantly altered in the METH-sensitized mice, including 582 that were down-regulated and five that were up-regulated (Figure 4A, Additional file 5). Several lincRNAs, including nuclear-enriched abundant transcript 2 (Neat2) [Genebank: AY722410], nuclear enriched abundant transcript 1 (Neat1) [Genebank: GQ859163], and myocardial infarction associated transcript (Miat) [Genebank: NR_033657], that were significantly down-regulated by METH have previously been characterized as possessing neurological functions [15,16].
Bidirectional lncRNAs are oriented head-to-head with a protein-coding gene within 1 kb, but are transcribed in the opposite direction. Bidirectional lncRNAs have been shown to affect the cis regulation of the nearby protein-coding genes potentially via promoter competition or the maintenance of an open chromatin structure [35,36]. We identified 74 aberrantly altered lncRNAs that formed bidirectional pairs with protein-coding genes in the current study (Figure 4A, Additional file 6) that included 72 down-regulated and two up-regulated lncRNAs.
LncRNAs in trans regulation
LncRNAs can work in a trans manner when they affect genes on other chromosomes . Previous studies have shown that lncRNAs can interact with associated mRNAs via the formation of complementary hybrids [37,38]. Therefore, using the RNAplex program , we subjected the differentially expressed known lncRNAs to trans-analysis and found that 64% (3994 of 6246) of the differentially expressed known lncRNAs were capable of acting in a trans manner, and 2386 of the associated genes have been found. We further investigated the networks formed by the trans-acting lncRNAs and their associated genes, which are termed the ‘many-to-many’ type; i.e., one lncRNA can have one or more associated genes. As shown in Figure 4B, 403 down-regulated lncRNAs exhibited one-to-one trans-regulation relationships with protein-coding genes (Additional file 7). In contrast, over 90% lncRNAs (3591 of 3994) might have more than one trans-associated gene. Additionally, the cis- and trans-acting lncRNAs also overlapped, and 3838 known lncRNAs were identified as having both cis- and trans-associated genes (Figure 4C).
Functional analyses of the cis- and trans-associated genes
To further investigate the potential effects of known lncRNAs on METH-induced locomotor sensitization, we subjected the cis- and trans-associated genes to GO and KEGG pathway analyses. Using P < 0.05 and FDR < 0.05 as cutoff for significance, we found that the predicted cis and trans-associated genes were significantly enriched in axon guidance, ubiquitin mediated proteolysis, neuron projection, the MAPK signalling pathway, long-term potentiation (LTP), long-term depression (LTD), calcium signalling pathway, dopaminergic synapse, and glutamatergic synapse; these processes are generally linked to neuronal development, neuronal plasticity, learning and memory, and reward and addiction (Figure 5, see full list in Additional files 8 and 9).
Previous studies have found that the expressions of lncRNAs are aberrantly altered in the NAc of cocaine-conditioned mice and heroin addicts, which suggests an important role of lncRNAs in drug addiction [9,19]. Nevertheless, the expression profiles of lncRNAs and their potential effects on METH-induced locomotor sensitization are largely unknown. Here, we used high-throughput ssRNA-seq technology to examine the alteration in the lncRNAs expression profile in the NAc of METH-sensitized mice. Using a stringent threshold for statistical significance (P < 0.001, FDR ≤ 0.0001 and an absolute value of the fold change ≥ 1.25), we identified numerous (6246) METH-regulated lncRNAs (Additional file 1), and 125 of these lncRNAs were also significantly altered in the NAc of cocaine-conditioned mice . Interestingly, the differentially expressed lncRNAs were less likely to be up-regulated and more likely to be down-regulated in the METH-sensitized mice. These data are consistent with the pattern of lncRNA expression in the NAc of cocaine-treated mice . Although the precise regulatory mechanism remains unclear, these results suggest that METH might reduce the expressions of lncRNAs. Further studies are needed to investigate the potential role of these differentially expressed lncRNAs. Nevertheless, our findings suggest that lncRNAs might be involved in METH addiction and provide new insight into the molecular mechanisms of METH abuse. To our knowledge, this is the first description of global lncRNAs expression profiling in the context of METH-induced behavioural sensitization. Moreover, our results provide numerous of METH-responsive lncRNA candidates for further functional research.
It has been shown that lncRNAs originate from complex loci that contain interlaced networks of long non-coding and protein-coding transcripts [40,41]. Analysis of the genomic characterization of lncRNAs is helpful in predicting their regulatory effects at the biological level. Indeed, a number of previously characterized lncRNAs have been proven to regulate the expressions of protein-coding genes that share genomic loci with the lncRNAs [30,42,43]. In the present study, we found that lncRNAs were associated with protein-coding genes in a variety of manners that included sense intronic, overlapping, lincRNA, NAT and bidirectional. From the perspective of the lncRNAs in their genomic context, our findings suggest the potential functions of these lncRNAs in terms of METH sensitization. Furthermore, lncRNAs have been reported to be involved in regulation of gene expression through trans-acting pathways in which they affect genes on other chromosomes. Here, we showed that numerous of METH-regulated lncRNAs interacted with associated protein-coding genes in trans manners, which suggests that these lncRNAs might be biologically meaningful. Notably, the trans-acting lncRNAs were found to form a “many-to-many” network with their associated genes, which reflects the complexity of the mechanisms of the regulation of METH-regulated lncRNAs. Interestingly, five lncRNAs (Kcnq1ot1, Zfhx2as, Neat1, Neat2, and Miat) that were found to be regulated by METH in our study have been reported to interact in cis or trans manners with targeted loci. For example, the antisense transcripts Kcnq1ot1 and Zfhx2as were found to regulate the expressions of their sense partners [29,30], which have been reported to be involved in the modulation of LTP in the hippocampus  and behavioural abnormalities , respectively. It has been reported that Neat1, Neat2 and Miat function as cofactors for pre-mRNA splicing by interacting with the splicing factors [46-48], and significantly, Neat2 appears to regulate neuronal plasticity by modulating the expressions of multiple synaptic genes , which suggests that lncRNAs-related nuclear modification might play a role in METH addiction through rapid post-transcriptional changes in gene expression. Notably, a preliminary examination of a published dataset based on heroin abusers revealed up-regulations of NEAT1, NEAT2 and MIAT , which is suggestive of differential responses across different drugs of abuse. Although the precise regulatory mechanism remains unclear, these well-characterized lncRNAs might play roles in terms of METH abuse. Taken together, the identification of cis and trans-acting lncRNAs suggests the potential functional implications of METH-regulated lncRNAs that might control the expressions of proximal or distal associated genes and thus contribute to METH-induced locomotor sensitization and addiction.
A major function of lncRNAs appears to be the control of the gene expression via cis- and trans-acting pathways . Thus, functional analyses of cis- and trans-associated genes are helpful for predicting the potential effects of lncRNAs on METH-induced locomotor sensitization. In the present study, we found that the predicted cis- and trans-associated genes were significantly enriched during neuronal development, neuronal plasticity, learning and memory, and reward and addiction (Figure 5). Previous studies have demonstrated that the neuronal plasticity that occurrs as a consequence of exposure to drugs of abuse plays a critical role in the modulation of persistent addictive behaviours [49-51]. Similarly, it has been reported that chronic exposure to drugs of abuse modulates learning and memory, which are thought to underlie rewarding and addictive behaviours [52-54]. Therefore, our findings suggest that the lncRNAs that were modified by METH in this study might influence the expressions of genes that are involved in neuronal plasticity, learning and memory and thus contribute to METH addiction.
Moreover, we unexpectedly found that sense intronic lncRNAs comprised the largest portion of the cis-lncRNAs. Previous studies have demonstrated that intronic lncRNAs can either be transcriptional segments of processed mRNAs or independent transcripts that are simply located within intron-annotated genomic regions, and intronic lncRNAs that originated from pre-mRNAs are thought to be the main actors in the regulation of gene expression [55,56]. However, current evidence indicates that intronic lncRNAs that exhibit independent transcription might also be biologically significant [24,57-59]. Although the ssRNA-seq technique used in this study cannot determine whether these intronic lncRNAs are alternative splicing products of a pre-mRNA or independent transcripts, METH-regulated sense intronic lncRNAs were found to originate from the intronic regions of the corresponding protein-coding genes, such as calcium/calmodulin-dependent protein kinase IV (Camk4), cAMP response element-binding protein 1(Creb1), CREB-binding protein (Crebbp), glutamate receptor, ionotropic, AMPA1 (alpha 1) (Gria1) and mitogen-activated protein kinase 10(Mapk10), which have been suggested to be responsible for synaptic transmission and specific signal transduction in long-term drug-induced neuroadaptation. Further experiments are needed to investigate the precise natures of these sense intronic lncRNAs. These results suggest that the sense intronic lncRNAs might play an important role in the regulation of METH-induced locomotor sensitization.
In summary, we reported a transcriptional profiling of lncRNAs in the NAcs of METH-sensitized and control mice. We have identified a number of METH-responsive lncRNAs. The predicted cis- and trans-associated genes of these METH-regulated lncRNAs were significantly enriched during the cellular and molecular events that contribute to reward and addiction. Although further experiments are needed to investigate the distinct function and the precise regulatory mechanism of each candidate lncRNA, our data suggest that exposure to METH elicits a global alterations in lncRNA expression in the NAc of sensitized mice that might be involved in METH-induced locomotor sensitization and addiction.
Adult wild-type C57BL/6 mice (7–8 weeks old, male, 20-25 g), purchased from Beijing Vital River Laboratory Animal Technolxogy Co. Ltd were used for these experiments. The mice were kept maintained in a regulated environment (23 ± 1°C, 50 ± 5% humidity) on a 12-h light/dark cycle (lights on from 7:00 am to 7:00 pm) and were handled in accordance with the Institutional Animal Care and Use Committee of Xi’an Jiaotong University. All efforts were made to minimize the number of animals used and to reduce stress from handling during the injections.
The METH hydrochloride used for the tests was purchased from the National Institute for the Control of Pharmaceutical and Biological Products (Beijing, P.R. China), and was dissolved in 0.9% physiological saline. The volume of the intraperitoneal (i.p.) injections was 10 ml/kg.
Procedure of METH exposure
The treatment regimens used in the current study have been shown to produce robust locomotor sensitization in our previous studies [2,6]. Briefly, the mice were given once-daily injections of saline for two consecutive days (day 1–2), after which they were randomly divided into two groups. The groups of mice were then given once-daily injections of METH (2 mg/kg) or saline for five consecutive days (day 3–7) followed by two injection-free days (day 8–9). On day 10, the mice were given a challenge injection of either 2 mg/kg METH or saline. Horizontal locomotor activity was performed on all drug treatment days for 60 minutes before and after the injections. The injections were performed in the open field test apparatus and during the light phase of the light/dark cycle. On the drug injection days, the mice were brought into the behaviour room 60 min prior to the beginning of the experiments to acclimate to the new environment. For all experiments, the mice were sacrificed 24 h after the final injection.
NAc sample preparation and RNA isolation
After the mice were sacrificed, the brains were micro-dissected and the NAcs were harvested and immediately frozen in liquid nitrogen. NAc lysates from eight mice from each group were pooled for total RNA isolation, following the instructions of the manufacturer of TRIzol (Invitrogen, USA). One saline and one METH sample were prepared from the RNAs that were extracted from the pooled NAc lysates and these samples were referred to as the saline and METH samples. The RNA qualities were evaluated with an Agilent 2100 BioAnalyzer (Agilent Technologies, USA), and all samples exhibited an RIN > 8.
Strand-specific cDNA library construction and sequencing
After the total RNA passed the RNA quality control for deep sequencing, we prepared to construct the strand-specific cDNA libraries. Briefly, the total RNA (5 μg) from each sample was fragmented into ~200 base pair (bp) units using a Covaris-S2 system after the removal of the rRNA without preselecting the mRNA. Next, the RNA fragments were used to generate double-stranded cDNA. The first cDNA strand was synthesized using random hexamers, and the second strand of cDNA was synthesized using deoxy-UTP instead of deoxy-TTP with DNA polymerase I. Then the double-stranded cDNAs were end-repaired after purification with a QiaQuick PCR kit, and the adapters were ligated. Subsequently, the uridine-containing strand was destroyed by uracil-N-glycosylase, which enabled the identification of the transcript orientation. Subsequently, to acquire the sequencing library products, the single-stranded adapted cDNA fragments of 200 bp were recovered and purified with agarose gel electrophoresis and then enriched by PCR for 12 cycles. The purified cDNA library products were evaluated using the Agilent 2100 BioAnalyzer and then sequenced on an Illumina HiSeq 2000.
SsRNA-seq data analyses
After sequencing, the raw reads that were generated by sequencers, were saved in the fastq format. To obtain reliable clean reads, the dirty raw reads were filtered according to four criteria: reads with sequence adaptors were removed; reads with more than 5% ‘N’ bases were removed; low-quality reads, in which more than 50% of the QA were ≤ 15 bases were removed; and ribosomal RNA sequences that were obtained from the ribosomal RNA database SILVA  by the software SOAP v2.2.0  were removed based on an allowance of no more than three mismatched bases. All subsequent analyses were based on clean reads. The clean reads of the saline and METH groups were separately aligned to the mouse genome, UCSC mm9  using the software TopHat v2.0.4 . Mismatches of no more than 5 bp were allowed in the alignment of each read. The resulting alignment data from TopHat were then assembled into transcripts by the assembler Cufflinks v2.0.0 . The assembled transcripts that corresponded to known lncRNAs were determined by perfect sequence matching to the NONCODE v3.0 database of known non-coding RNA . Furthermore, the assembled transcripts were aligned to protein databases, including KEGG Orthology , non-redundant protein database , COG  and UniProtKB/Swiss-Prot , to obtain protein-coding transcripts. To obtain novel non-coding transcript candidates, the remaining transcripts that were not matched to any known lncRNAs or protein-coding sequences were then used to predict their abilities to encode proteins using the Coding Potential Calculator (CPC) , which based on framefinder. Thus, the transcripts that lacked the ability to encode proteins were considered as novel lncRNA candidates.
The clean reads that were uniquely mapped to lncRNAs were used to calculate the expression levels. The relative expression levels of the lncRNAs in the saline and METH groups were measured as the number of uniquely mapped reads per kilobase per million mappable reads (RPKM). The formula was defined as follows: RPKM = 106 × C/(NL × 10−3), where C was the number of reads that uniquely mapped to the given transcript, N was the number of reads that uniquely mapped to all transcripts, and L was the total length of the given transcript. The RPKM method eliminates the influences of different transcript lengths and sequencing discrepancies on the calculation of expression. Therefore, the RPKM value was directly used to compare the differences in lncRNA expressions between samples. The fold change from the normalized expression was calculated as log2 (RPKMMETH/RPKMsaline) to assess the levels. Because we only had one replicate per group, the variances of regulated levels were directly estimated from the RPKM values using the Poisson distribution, and P-values were calculated . Therefore, to compensate for false-positive findings at each significance threshold, we calculated a false discovery rate (FDR; Benjamini-Hochberg) for each lncRNA and applied it for genome-wide corrections . We identified lncRNAs that were differentially regulated between the METH and saline groups based on the following criteria: P < 0.001, FDR ≤ 0.0001 and absolute value of the fold change ≥ 1.25.
Cis and trans analyses
Identification of the genes associated with differentially expressed lncRNAs via cis-or trans-regulation might provide insight into the potential functions of lncRNAs. We subjected the significantly changed known lncRNAs to cis and trans analyses. For the cis analyses, we classified the differentially expressed lncRNAs into the following 5 categories according to their genomic contexts relative to protein-coding genes: sense intronic, overlapping, NAT, lincRNA and bidirectional. For the trans analyses, we predicted the trans-associated genes of the differentially expressed lncRNAs with RNAplex v0.2, which is a fast tool for RNA-RNA interaction searches by neglecting intramolecular interactions and by using as lightly simplified energy model . The RNAplex parameters were set as –e < −20 in the current study to identify the trans-associated genes , and genes that were found to be located on that same chromosome as the lncRNA were excluded.
GO and KEGG pathway enrichment analyses
For the GO and KEGG pathway enrichment analyses, the cis- and trans-associated genes of the lncRNAs that were significantly modified in the METH-sensitized mice were analysed with the functional annotation tool Blast2GO . First, all associated genes were mapped to GO terms and pathways in the gene ontology  and KEGG pathway  databases by calculating the gene numbers for every term and pathway. Then the P values were calculated via hypergeometric tests and go through a correction with FDR. P < 0.05 and FDR < 0.05 were used as thresholds for defining significantly enriched GO terms and pathways. Finally, the associated genes that corresponded to specific biological functions were filtered.
QPCR was performed on the RNAs isolated from the NAcs of individual mice (n = 11-15 per group). First-strand cDNA was synthesized using the Thermo Scientific RevertAid First Strand cDNA Synthesis Kit (Thermo Scientific, USA). Based on the manufacturer’s instructions and suggested parameters (25°C 5 min, 42°C 60 min, and 70°C 5 min), 500 ng of total RNA from each sample was utilised. QPCR was performed on Bio-Rad iQ™5 system real-time detection instrument (Bio-Rad, USA) using SYBR Premix Ex Taq II (TaKaRa Biotechnology, Japan) in the following conditions: 95°C for 30 sec; 95°C for 10 sec and 60°C for 1 min, which were 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. Fold changes in expression were calculated with a log 2 transform. Independent-sample t-tests were used to test for significant differences (SPSS v17.0, SPSS Inc., USA). The primer pairs are shown in Table 3. Pearson’s coefficient analysis was also performed with SPSS (version 17.0, USA).
Availability of supporting data
The raw sequences have been deposited in the ArrayExpress database (www.ebi.ac.uk/arrayexpress) under accession number E-MTAB-2843.
Calcium/calmodulin-dependent protein kinase IV
Coding potential calculator
cAMP response element-binding protein 1
False discovery rate
Glutamate receptor, ionotropic, AMPA1 (alpha 1)
Potassium voltage-gated channel, subfamily Q, member 1
Potassium voltage-gated channel, subfamily Q, member 1, opposite strand/antisense transcript 1
Long intergenic non-coding RNA
Long non-coding RNA
Mitogen-activated protein kinase 10
Myocardial infarction associated transcript
Natural antisense transcripts
Nuclear enriched abundant transcript 1
Nuclear-enriched abundant transcript 2
quantitative real-time PCR
Reads per kilobase per million mappable reads
strand-specific complementary DNA sequencing technology
Zinc finger homeobox 2
Zinc finger homeobox 2, antisense
TE Robinson KB. Structural plasticity associated with exposure to drugs of abuse. Neuropharmacology. 2004;47 Suppl 1:33–46.
Zhu J, Chen Y, Zhao N, Cao G, Dang Y, Han W, et al. Distinct roles of dopamine D3 receptors in modulating methamphetamine-induced behavioral sensitization and ultrastructural plasticity in the shell of the nucleus accumbens. J Neurosci Res. 2012;90(4):895–904.
Chen L, Xu M. Dopamine D1 and D3 receptors are differentially involved in cue-elicited cocaine seeking. J Neurochem. 2010;114(2):530–41.
Kong H, Kuang W, Li S, Xu M. Activation of dopamine D3 receptors inhibits reward-related learning induced by cocaine. Neurosci. 2011;176:152–61.
Ren Z, Sun WL, Jiao H, Zhang D, Kong H, Wang X, et al. Dopamine D1 and N-methyl-D-aspartate receptors and extracellular signal-regulated kinase mediate neuronal morphological changes induced by repeated cocaine administration. Neurosci. 2010;168(1):48–60.
Zhao N, Chen Y, Zhu J, Wang L, Cao G, Dang Y, et al. Levo-tetrahydropalmatine attenuates the development and expression of methamphetamine-induced locomotor sensitization and the accompanying activation of ERK in the nucleus accumbens and caudate putamen in mice. Neurosci. 2014;258:101–10.
Hollander JA, Im H, Amelio AL, Kocerha JK, Bali P, Lu Q, et al. Striatal microRNA controls cocaine intake through CREB signaling. Nature. 2010;466:197–202.
Eipper-Mains J, Kiraly DD, Palakodeti D, Mains RE, Eipper BA, Graveley BR. microRNA-Seq reveals cocaine-regulated expression of striatal microRNAs. RNA. 2011;17(8):1529-1543.
Bu Q, Hu Z, Chen F, Zhu R, Deng Y, Shao X, et al. Transcriptome analysis of long non-coding RNAs of the nucleus accumbens in cocaine-conditioned mice. J Neurochem. 2012;123:790–9.
Chekulaeva M, Filipowicz W. Mechanisms of miRNA-mediated post-transcriptional regulation in animal cells. Curr Opin in Cell Biol. 2009;21(3):452–60.
Kornienko AE, Guenzl PM, Barlow DP, Pauler FM. Gene regulation by the act of long non-coding RNA transcription. BMC Biol. 2013;11:59.
Im H, Hollander JA, Bali P, Kenny PJ. MeCP2 controls BDNF expression and cocaine intake through homeostatic interactions with microRNA-212. Nat Neurosci. 2010;13(9):1120–7.
Chandraskear V, Dreyer J. Regulation of miR-124, let-7d, and miR-181a in the accumbens affects the expression, extinction, and reinstatement of cocaine-induced conditioned place preference. Neuropsychopharmacology. 2011;36:1149–64.
Kapranov P, Cheng J, Dike S, Nix DA, Duttagupta R, Willingham AT, et al. RNA maps reveal new RNA classes and a possible function for pervasive transcription. Science. 2007;316:1484–8.
Bernard D, Prasanth VK, Tripathi V, Colasse S, Nakamura T, Xuan Z, et al. A long nuclear-retained non-coding RNA regulate synaptogenesis by modulating gene expression. EMBO J. 2010;29:3082–93.
Mercer TR, Qureshi IA, Gokhan S, Dinger ME, Li G, Mattick JS, et al. Long noncoding RNAs in neuronal-glial fate specification and oligodendrocyte lineage maturation. BMC Neurosci. 2010;11:14.
Faghihi MA, Modarresi F, Khalil AM, Wood DE, Sahagan BG, Morgan TE, et al. Expression of a noncoding RNA is elevated in Alzheimer’s disease and drives rapid feed-forward regulation of beta-secretase. Nat Med. 2008;14:723–30.
Chubb JE, Bradshaw NJ, Soares DC, Porteous DJ, Millar JK. The DISC locus in psychiatric illness. Mol Psychiatr. 2008;13:36–64.
Michelhaugh SK, Lipovich L, Blythe J, Jia H, Gregory K, Bannon MJ. Mining Affymetrix microarray data for long noncoding RNAs: altered expression in the nucleus accumbens of heroin abusers. J Neurochem. 2011;116(3):459–66.
UCSC. mm9. http://genome.ucsc.edu/.
NONCODE V3.0. http://www.noncode.org/.
Huang Y, Liu N, Wang JP, Wang YQ, Yu XL, Wang ZB, et al. Regulatory long non-coding RNA and its functions. J Physiol Biochem. 2012;68:611–8.
Knauss JL, Sun T. Regulatory mechanisms of long noncoding RNAs in vertebrate central nervous system development and function. Neurosci. 2013;235:200–14.
Mercer TR, Dinger ME, Sunkin SM, Mehler MF, Mattick JS. Specific expression of long noncoding RNAs in the mouse brain. Proc Natl Acad Sci U S A. 2008;105:716–21.
Wagner EJ, Carpenter PB. Understanding the language of Lys36 methylation at histone H3. Nat Rev Mol Cell Biol. 2012;13:115–26.
Conte C, Dastugue B, Vaury C. Promoter competition as a mechanism of transcriptional interference mediated by retrotransposons. EMBO J. 2002;21:3908–16.
Katayama S, Tomaru Y, Kasukawa T, Waki K, Nakanishi M, Nakamura M, et al. Antisense transcription in the mammalian transcriptome. Science. 2005;309:1564–6.
Werner A. Biological functions of natural antisense transcripts. BMC Biol. 2013;11:31.
Lewis A, Green K, Dawson C, Redrup L, Huynh KD, Lee JT, et al. Epigenetic dynamics of the Kcnq1 imprinted domain in the early embryo. Development. 2006;133:4203–10.
Komine Y, Nakamura K, Katsuki M, Yamamori T. Novel transcription factor zfh-5 is negatively regulated by its own antisense RNA in mouse brain. Mol Cell Neurosci. 2006;31:273–83.
Guttman M, Amit I, Garber M, French C, Lin MF, Feldser D, et al. Chromatin signature reveals over a thousand highly conserved large non-coding RNAs in mammals. Nature. 2009;458:223–7.
Ponting CP, Oliver PL, Reik W. Evolution and functions of long noncoding RNAs. Cell. 2009;136:629–41.
Kim TK, Hemberg M, Gray JM, Costa AM, Bear DM, Wu J, et al. Widespread transcription at neuronal activity-regulated enhancers. Nature. 2010;465:182–7.
Ørom UA, Shiekhattar R. Long non-coding RNAs and enhancers. Curr Opin Genet Dev. 2011;21:194–8.
Trinklein ND, Aldred SF, Hartman SJ, Schroeder DI, Otillar RP, Myers RM. An abundance of bidirectional promoters in the human genome. Genome Res. 2004;4:62–6.
Wei W, Pelechano V, Järvelin AI, Steinmetz LM. Functional consequences of bidirectional promoters. Trends Genet. 2011;27:267–76.
Yoon JH, Abdelmohsen K, Srikantan S, Yang X, Martindale JL, De S, et al. LincRNA-p21 suppresses target mRNA translation. Mol Cell. 2012;47:648–55.
Gong CG, Maquat LE. LncRNAs transactivate Staufen1-mediated mRNA decay by duplexing with 3'UTRs via Alu elements. Nature. 2011; 470:284–288.
Tafer H, Hofacker IL. RNAplex: a fast tool for RNA-RNA interaction search. Bioinformatics. 2008;24:2657–63.
Derrien T, Johnson R, Bussotti G, Tanzer A, Djebali S, Tilgner H, et al. The GENCODE v7 catalog of human long noncoding RNAs: analysis of their gene structure, evolution, and expression. Genome Res. 2012;22:1775–89.
Engström PG, Suzuki H, Ninomiya N, Akalin A, Sessa L, Lavorgna G, et al. Complex loci in human and mouse genomes. PLoS Genet. 2006;2:e47.
Feng J, Bi C, Clark BS, Mady R, Shah P, Kohtz JD. The Evf-2 noncoding RNA is transcribed from the Dlx-5/6 ultraconserved region and functions as a Dlx-2 transcriptional coactivator. Gene Dev. 2006;20:1470–84.
Martianov I, Ramadass A, Barros AS, Chow N, Akoulitchev A. Repression of the human dihydrofolate reductase gene by a non-coding interfering transcript. Nature. 2007;445:666–70.
Petrovic MM, Nowacki J, Olivo V, Tsaneva-Atanasova K, Randall AD, Mellor JR. Inhibition of post-synaptic Kv7/KCNQ/M channels facilitates long-term potentiation in the hippocampus. PLoS One. 2012;7:e30402.
Komine Y, Takao K, Miyakawa T, Yamamori T. Behavioral abnormalities observed in Zfhx2-deficient mice. PLoS One. 2012;7:e53114.
Hutchinson JN, Ensminger AW, Clemson CM, Lynch CR, Lawrence JB, Chess A. A screen for nuclear transcripts identifies two linked noncoding RNAs associated with SC35 splicing domains. BMC Genomics. 2007;8:39.
Tollervey JR, Curk T, Rogelj B, Briese M, Cereda M, Kayikci M, et al. Characterizing the RNA targets and position-dependent splicing regulation by TDP-43. Nat Neurosci. 2011;14:452–8.
Tsuiji H, Yoshimoto R, Hasegawa Y, Furuno M, Yoshida M, Nakagawa S. Competition between a noncoding exon and introns: gomafu contains tandem UACUAAC repeats and associates with splicing factor-1. Genes Cells. 2011;16:479–90.
Li Y, Kauer JA. Repeated exposure to amphetamine disrupts dopaminergic modulation of excitatory synaptic plasticity and neurotransmission in nucleus accumbens. Synapse. 2004;51:1–10.
Gerdeman GL, Partridge JG, Lupica CR, Lovinger DM. It could be habit forming: drugs of abuse and striatal synaptic plasticity. Trends Neurosci. 2003;26:184–92.
Dietz DM, Dietz KC, Nestler EJ, Russo SJ. Molecular mechanisms of psychostimulant-induced structural plasticity. Pharmacopsychiatry. 2009;42:S69–78.
Whitlock JR, Heynen AJ, Shuler MG, Bear MF. Learning induces long-term potentiation in the hippocampus. Science. 2006;313:1093-1097.
Mitsushima D, Sano A, Takahashi T. A cholinergic trigger drives learning-induced plasticity at hippocampal synapses. Nat Commun. 2013;4:2760.
Steven E, Hyman MD. Addiction: a disease of learning and memory. Am J Psychiatry. 2005;162:1414–22.
Mattick JS. Introns: evolution and function. Curr Opin Genet Dev. 1994;4:823–31.
Mattick JS, Gagen MJ. The evolution of controlled multitasked gene networks: the role of introns and other noncoding RNAs in the development of complex organisms. Mol Biol Evol. 2001;18:1611–30.
Martone R, Euskirchen G, Bertone P, Hartman S, Royce TE, Luscombe NM, et al. Distribution of NF-kappa B-binding sites across human chromosome 22. Proc Natl Acad Sci U S A. 2003;100:12247–52.
Cawley S, Bekiranov S, Ng HH, Kapranov P, Sekinger EA, Kampa D, et al. Unbiased mapping of transcription factor binding sites along human chromosomes 21 and 22 points to widespread regulation of noncoding RNAs. Cell. 2004;116:499–509.
Euskirchen G, Royce TE, Bertone P, Martone R, Rinn JL, Nelson FK, et al. CREB binds to multiple loci on human chromosome 22. Mol Cell Biol. 2004;24:3804–14.
Li R, Yu C, Li Y, Lam TW, Yiu SM, Kristiansen K, et al. SOAP2: an improved ultrafast tool for short read alignment. Bioinformatics. 2009;25(15):1966–7.
Kim D, Pertea G, Trapnell C, Pimentel H, Kelley R, Salzberg SL. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 2013;14(4):R36.
Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, van Baren MJ, et al. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol. 2010;28(5):511–5.
Non-redundant protein database. ftp://ftp.ncbi.nih.gov/blast/db/FASTA/.*.
Coding Potential Calculator. http://cpc.cbi.pku.edu.cn/.
Audic S, Claverie JM. The significance of digital gene expression profiles. Genome Res. 1997;7(10):986–95.
Benjamini Y, Yekutieli D. The control of the false discovery rate in multiple testing under dependency. Ann Stat. 2001;29(4):1165–88.
Gene ontology. http://www.geneontology.org/.
KEGG pathway. http://www.genome.jp/kegg/pathway.html.
This research was supported by a grant from the National Natural Science Foundation of China (No. 81172913) to TC and the Fundamental Research Funds for the Central Universities. We thank Na Zhao, Nan Dong, and Jiaqi Li for stimulating discussions and for providing animal care.
The authors declare that they have no competing interests.
TC initiated the project; LZ, JZ and TC designed the experiments; YL, YC, YL, LH, SC, TL, and YD performed the experiments; YL, LH, and SC performed the computational analyses; LZ, JZ, and TC participated in the experimental analyses; LZ, JZ, and TC wrote the manuscript. All authors have read and approved this manuscript for publication.
Li Zhu and Jie Zhu contributed equally to this work.
Full list of the METH-regulated lncRNAs.
Full list of the sense intronic lncRNAs that were significantly regulated by METH.
Full list of the overlapping lncRNAs that were significantly regulated by METH.
Full list of the NATs that were significantly regulated by METH.
Full list of the lincRNAs that were significantly regulated by METH.
Full list of the bidirectional lncRNAs that were significantly regulated by METH.
Full list of the one-to-one trans -acting lncRNAs that were significantly regulated by METH.
Full list of the GO and pathways of the cis- acting lncRNAs and their associated genes.
Full list of the GO and pathways of the trans -acting lncRNAs and their associated genes.