Integrated miRNA-mRNA analysis in the habenula nuclei of mice intravenously self-administering nicotine

A considerable amount of evidence suggests that microRNAs (miRNAs) play crucial roles in the neuroadaptation of drug addiction. Habenula (Hb), one of the critical brain regions involved in reward and addiction, can be divided into two anatomically and transcriptionally distinct regions: medial habenula (MHb) and lateral habenula (LHb) nuclei. However, very few studies have compared the functional roles of these regions. Here, by using mirConnX integrator and KEGG pathway mapping, we simultaneously analysed the differential expression patterns of miRNAs and messenger RNA (mRNA) within MHb and LHb under nicotine addiction. Significantly altered miRNAs and mRNAs were found in the Hb of mice intravenously self-administering nicotine. Interestingly, some miRNAs were oppositely regulated between the MHb and the LHb, and their potential targets included various genes of cell signalling pathways related to the degeneration of fasciculus retroflexus (FR). This study provides an improved insight into the differential regulation of habenular transcripts in nicotine addiction, as well as the potential functions of miRNAs in several biological pathways involved in the nicotine addiction.

in both aversion-and reward-related behaviors, while its projection to the nucleus accumbens (NAcc) has been controversial because the functions of the shell and the core regions of the NAcc have not been clearly resolved 7,8 . The VTA regulation has also been studied with regard to the MHb. It is noteworthy that the projection from the MHb to the IPN, which is called fasciculus retroflexus (FR), has been identified to play a major role in nicotine addiction. IPN sends its GABAergic inhibitory projection to the VTA, and the excitatory signals from the MHb to the IPN via glutamate and acetylcholine signalling could inhibit the VTA functionally [9][10][11] . Accordingly, previous studies on the function of the Hb complex under the state of nicotine addiction have focused on how nAChRs play a role in these circuitries. However, molecular dynamics other than nAChR functions have been undervalued and still poorly understood. In this study, the expression changes of miRNAs and mRNAs in both the MHb and the LHb nuclei were analysed using microarrays after the intravenous self-administration of nicotine. Microarray experiments have been widely used for the systemic analysis of various target molecules, including miRNAs or mRNAs, by measuring the intensity of their expressions. With multiple microarrays from different experimental groups, detailed molecular profiles can be extracted and comparatively analysed. Therefore, in this study, the integrated miRNA-mRNA expression analysis was adopted to compare the transcriptional networks of the MHb versus the LHb nuclei and show the various molecular regulations and the functional implications of the Hb in nicotine addiction.

Results
Mice intravenously self-administering nicotine showed the general drug-seeking behavior. Prior to the nicotine infusion, mice underwent 10 days of food training in the operant self-administration (SA) chamber during which the mice learned to obtain the food pellet reward when they pressed the 'active' lever instead of the 'inactive' lever ( Fig. 1a,b). Subsequently, the reward was changed from the food pellets to the nicotine infusion (0.03 mg/kg per infusion). Mice consumed ~40 nicotine rewards during their first SA session, but rapidly adjusted their level of responding over the next 2-3 sessions to obtain ~11 rewards during the consecutive 1 h session (Fig. 1c). Mice stably pressed the 'active' lever significantly more than 'inactive' lever throughout the nicotine infusion experiment (Fig. 1d). When the unit dose of nicotine available was varied (0, 0.01, 0.25, 0.1 and 0.4 mg/kg per infusion), mice received nicotine infusions according to an inverted U-shaped dose-response curve (Fig. 1e). Compared to the active lever, inactive lever responses remained low across all nicotine doses (Fig. 1f). Interestingly, mice continued to increase drug intake throughout all unit doses of nicotine (0.01, 0.25, 0.1 and 0.4 mg/ kg per infusion) (Fig. 1g).
Microarray profiling of the habenular transcripts showed differentially altered expressions in mice intravenously self-administering nicotine. The MHb and the LHb exhibited differential expression profiles of numerous miRNAs and mRNAs after intravenous nicotine SA (Fig. 2a,b). Among 30434 miRNAs analysed, a total of 44 miRNAs were found to be significantly altered. In the MHb, 8 miRNAs increased more than two-fold and 7 miRNAs decreased under half compared to the drug-naïve group. In the LHb, 10 miRNAs increased more than two-fold, and 8 miRNAs decreased under half. Five miRNAs increased in both the MHb and the LHb. It is interesting to note that 18 miRNAs were altered oppositely with different degrees in the MHb and the LHb ( Table 1). The alterations of some miRNAs were additionally confirmed by qRT-PCR assay (Fig. 3a,b). Alterations of mmu-miR-3078-5p, mmu-miR-200c-3p, mmu-miR-496a-5p, mmu-miR-412-5p, and mmu-miR-323-5p were tested in the MHb and the LHb with the same samples used in the microarrays, and their altered expression patterns were well-matched with the array data (Fig. 3).
A mirConnX analysis of the miRNA-mRNA network identified various targets of forty-four miRNAs altered by the mice self-administering nicotine. To confirm the correlations between the array of the miRNAs and the mRNAs, 44 altered miRNAs were analysed using the mirConnX web interface (http://mirconnx.csb.pitt.edu) with the Pearson's correlation as a measurement of association. The representative networks of miRNA-mRNA interactions were derived from the miRNA target predictions, and the interaction patterns were analysed within the total network (Fig. 4). It is noteworthy that both mmu-miR-3078 and mmu-miR-467e were associated within a single network regulating four common targets, Ush1c, Olfr152, Lrrc31 and Cts8. (Fig. 4a). In addition, mmu-miR-669c, mmu-miR-467b and mmu-miR-376b were associated with other transcription factors for the regulation of the common targets within each regulation network.
Altered miRNAs have putative targets overlapping with the mRNA targets related to the various functions of the MHb and the LHb. The 44 miRNAs could target ~3000 genes based only on the sequence complementarity in the miRNA target database such as DIANA database and microRNA. org. Among these genes, 615 mRNAs were found to be altered in the mRNA microarray profiling. These mRNAs had reverse patterns of alteration with the corresponding 44 miRNAs, which could be interpreted as the repression of the mRNAs by the miRNAs (Suppl . Table S8). It is notable that there were 10 oppositely regulated mRNAs between the MHb and the LHb whose corresponding miRNAs were also oppositely altered (Table 2). Next, various categories of KEGG pathways were analysed through DAVID (e,f) Dose-response test of nicotine SA showed inverted U-shaped dose-response curve with the highest intake at the 0.1 mg/kg per infusion. Each dose was tested for 5 days and the reinforcement numbers of the final 3 days were averaged. (g) Total quantity of nicotine infused at each dose was calculated. All data are presented as mean ± SEM of nicotine infusions or lever presses (* P< 0.05, **P < 0.01, ****P < 0.0001).  Continued functional annotation tool (Fig. 5). Among the miRNA targets within the MHb, the neurotrophin signalling pathway and the neuroactive ligand-receptor interaction were found to be altered by intravenous nicotine SA. Among the miRNA targets within the LHb, various pathways including retinol metabolism and drug metabolism were found to be changed. For the analysis of the targets of oppositely regulated miRNAs in the MHb versus the LHb, the MAPK signalling pathway and the calcium signalling pathway were found to be altered by intravenous nicotine SA.

Discussion
This is the first study to show the characteristic expression patterns of transcripts (miRNAs and mRNAs) for the MHb and the LHb in a mouse model of nicotine addiction. We analysed three patterns (up-regulated, down-regulated and oppositely altered) of the miRNA expression profiles in the MHb and the LHb compared to the drug-naïve group for the functional study of nicotine addiction. By exploring the three patterns of miRNA expression profiles, we found significant differences in terms of gene ontology features. In addition, we found that the expression of up-regulated miRNAs in the LHb showed more dramatic changes than the changes in the MHb. Furthermore, we discovered an interesting aspect of the oppositely altered miRNAs in the MHb versus the LHb showing that the targeted mRNA networks are mainly involved in the development of nicotine addiction.  Table 1. Forty-four nicotine-responsive miRNAs in Hb. Forty-four microRNAs selected from the array data were shown. Each genomic locus within the chromosome and sequence of the miRNAs were represented. Fold-changes and detection p-values of the array were listed. In addition to the up-or down-regulated miRNAs in MHb and LHb, miRNAs that were oppositely regulated between MHb and LHb were also listed. The p cutoff value was P < 0.01 except for the opposite pattern. The p cutoff value of the opposite pattern was P < 0.05 except for the 2 microRNAs with the p-value of 0.053792 and 0.053524, respectively.  At present, the MHb is known to regulate various behavioral functions, including stress responses, depression and addiction 12 . The function of MHb in nicotine addiction has been elucidated from the perspectives of the neurotoxicity and the nAChR functions. A previous report on nicotine-induced neurotoxicity revealed that the neuronal projection from the MHb to the IPN, which is called FR, is degenerated by nicotine treatment 13 . Consequently, the nicotine-induced FR degeneration has been interpreted as a loss of forebrain control by drug abuse, possibly leading to the addiction of the nicotine 14 . Additionally, several recent studies have revealed that highly localised nAChR expression in the MHb regulates the condition of nicotine addiction 15,16 . According to our KEGG pathway analyses, altered miRNA-mRNA networks including 15 miRNAs and 156 mRNAs in the MHb appeared to regulate immune-response pathways including cytokine-cytokine receptor interaction and T cell receptor signalling pathways and neuroactive ligand-receptor interaction between GABA and GABA receptor (Suppl. Figs S1-S3). In general, it has been suggested that the increased incidence of various diseases in smokers may be due to altered immune responses rising from the chronic inhalation of chemicals (e.g. nicotine) in cigarette smoke 17 . In this respect, nicotine is known to induce T-cell anergy and immunosuppression 18 . Therefore, the alteration of these pathways observed in our study could be interpreted as missing links between the nicotine neurotoxicity and the signalling pathway of nicotine which were found to be regulated by the nicotine responsive miRNA-mRNA in our study 19,20 . In this regard, it is also possible that mmu-miR-200c-3p, which was up-regulated in the MHb by intravenous nicotine SA, possibly regulates phospholipase C, gamma 1 (Plcg1), which is shown to regulate immune-response pathways, namely the T-cell receptor signalling pathway and neurotrophin signalling pathway (Suppl. Figs S2 and S3).
Multiple studies of the LHb nucleus have revealed its regulatory function on the dopaminergic pathways of reward circuit in the brain, particularly concerning the negative prediction errors in the behavioral studies 7,8,21 . In intravenous nicotine SA, negative states of addiction could be induced in the brain because the mice experienced both reinforcement and withdrawal during the period of nicotine SA. In KEGG analyses of the LHb mRNA profile, the gene regulation induced by the PPAR signalling was involved with nicotine addiction (Suppl. Fig. S4). A previous report found that nicotine withdrawal could be modulated by the PPARα ligands and the PPARα agonists 22 , while another study using the inhibitor of PPARγ reported that the nicotine-induced gene expression could be regulated by the PPAR signalling pathway through upregulation of CD36 signaling 23 .
Opposite alteration patterns of the miRNA-mRNA networks between the MHb and the LHb also provided an interesting interpretation regarding the role of the Hb in mice intravenously self-administering nicotine. Our data revealed that several KEGG pathways, including MAPK signalling, calcium signalling and neurotrophin signalling, were oppositely regulated between the MHb and the LHb by the oppositely altered miRNAs (Suppl. Figs S5 and S6). It is noteworthy that these pathways have been studied in respect to the effects of nicotine [24][25][26] . On the other hand, mmu-miR-467c-3p is also oppositely regulated in the MHb vs the LHb and targets the brain-derived neurotrophic factor (Bdnf), a signalling protein in the central nervous system that is well-known as a positive regulator of substance abuse such as cocaine 27 . In addition, mmu-miR-721 was found to target the mitogen-activated protein kinase 1 (Mapk1), a well-known kinase acting on the intersecting point of various biochemical signalling pathways including differentiation, transcription and development. It is noteworthy that nuclear factor-kB (NF-kB) and   cyclic AMP response element binding protein (CREB) were shown to function in downstream of MAPK signalling (Suppl. Fig. S5) and these factors have been extensively studied in the drug abuse and addiction 28 . Furthermore, mmu-miR-669c-3p regulated the sortilin 1 (Sort1), one of the receptors involved in cell-death signalling, and this regulation may underlie the mechanisms of FR degeneration 29 . Finally, all these targets-Bdnf, NF-kB, CREB and Sort1-were critical regulators of the neurotrophin pathway, and various aspects of the relationship between nicotine addiction and the neurotrophin pathway have been studied in dozens of publications especially dealing with tyrosine kinase (Trk) receptor A, postnatal rat hippocampus and Trk receptor signalling [30][31][32] . Altogether, the bioinformatic analysis of the habenular miRNA-mRNA network in mice intravenously self-administering nicotine revealed several useful interpretations about the role of Hb in the nicotine addiction and the neuroadaptive changes by nicotine. In future, numerous in vivo experiments will be required to elucidate these interrelations linking those biological processes and miRNA-mRNA networks.
For both the nicotine-administering group and the age-matched drug naïve control group, mice were individually housed for 12 hours with a reversed light-dark cycle in a laboratory breeding room at the Korea Institute of Science and Technology (KIST). For the nicotine-administering group, water was freely available but foods were mildly restricted throughout the test (about 85-90% of their free-feeding body weight). For age-matched drug naïve control group, water and foods were provided ad libitum. Nicotine SA was conducted in the operant SA chamber enclosed by the sound attenuating cubicle (MED-307A-CT-D1, Med Associates, Inc., St. Albans, USA) in the dark room. All procedures regarding the use and the handling of the animals were conducted as approved by the Institutional Animal Care and Use Committee of the KIST.
Intravenous self-administration procedure. Surgery. Operant nicotine SA for mice was performed according to the method described previously 33 . After completing the ten-day food training with the food pellets in the SA chamber, the mice were anaesthetised by intraperitoneally injecting a mixture of ketamine (120 mg/kg) and xylazine (0.6 mg/kg). Then a catheter was implanted in the jugular vein of the mice and the end of the catheter was placed on the backs of the mice. To prevent infection, instruments used in the surgery were sterilised and the mice were injected daily with gentamycin if necessary.
Apparatus. Operant chambers (29.5 cm × 32.5 cm × 23.5 cm) with two retractable levers were used in the tests. The left lever was used as an active lever to confer the reward. Cue light was located above the lever and programmed to be turned on when the lever was pressed correctly. For nicotine infusion, nicotine syringe pumps were connected to the outlet on the backs of the mice via metal spring-covered tubes. The operant chambers were maintained and controlled by MED-PC software (Med Associates).
Procedure. Mice were trained to discriminate between the levers and to press the active lever instead of the inactive lever in order to deliver the food reward in the training program prior to the SA of nicotine. The fixed-ratio (FR) was escalated from FR1 to FR5 gradually when the mice received more than 25 food pellets for two consecutive days. After finishing the FR5 food training, the mice had intravenous catheter surgery into the right jugular vein to allow for the nicotine delivery, followed by three days of recovery. At the nicotine SA session, the mice obtained a 0.03 mL nicotine solution as a reward for every active lever press, with FR5 (delivered over three seconds) in a one-hour session. When they were infused with the nicotine solution, a cue light above the active lever was turned on for 20 seconds. When the mice showed a stable infusion in the unit doses of 0.03 mg/kg per infusion, 0.1 mg/kg per infusion was tested and used as a 'training dose' between the change of each dose. Dose-response test was done using unit doses of 0, 0.01, 0.25 and 0.4 mg/kg per infusion, with each dose tested over five days and the average of the final three days calculated for the reward number. A latin-square cross-over design was used in the test to control any effect of test order. After dose-response test, 0.03 mg/kg per infusion dose of nicotine was infused to the mice during additional 10 days for optimal induction of nicotine addiction prior to sacrifice. Nicotine solutions were prepared with saline and (-)-nicotine hydrogen tartrate salt (Sigma, St. Louis, MO, USA).

RNA extraction from the MHb and the LHb and microarray experiments.
Mice were decapitated on the next day of the final SA of nicotine. Frozen whole brains were cut into 100 μ m slices with the coronal plane on a cryostat (CM3050S, Leica) at − 20 °C. The MHb and the LHb were micro-dissected separately (Suppl. Fig. S7), and the total RNA was isolated using the RNA STAT-60 (Amsbio, Abingdon, UK) for the microarray experiments according to the manufacturer's instructions. The concentration and the purity of the RNA samples were checked using a Nanodrop spectrophotometer (Thermo Scientific, Waltham, MA, USA). The miRNA expression profiling was performed using Affymetrix GeneChip ® miRNA 4.0 arrays (Affymetrix, Santa Clara, CA, USA) containing 30434 mature miRNAs. The mRNA expression profiling was performed using the Affymetrix GeneChip ® Mouse Gene 2.0 ST Array containing 770,317 probes. With regard to miRNA expression profiling, the RNA was labelled using the FlashTag Biotin HSR (Genisphere, Hatfield, PA, USA) and was then hybridised to the miRNA array. After hybridization, staining and washing were performed according to the user guide. For mRNA expression profiling, the RNA was reverse transcribed to double-stranded cDNA, fragmented and labelled using the Biotin labelling kit, and then hybridized to the gene array as recommended. Standard Affymetrix array cassette staining, washing and scanning were then performed. In the last step for scanning the signals and analysing the data, Affymetrix ® Expression Console Software (version 1.2.1) was used for the microarray analysis. Raw data (CEL files) were normalised at the transcript level by using a robust multi-average (RMA) method. Data analysis. Normalised data were analysed from each array. The fold change and the detection of p-values were used to screen the miRNAs and the mRNAs that showed significantly different expressions. Hierarchical clustering of the miRNAs and mRNAs with significantly different expressions was performed using the Cluster 3.0 software, and was visualised using MultiExperiment Viewer (www.tm4.org).
Integrated analysis of the miRNA targets. DIANA database (http://www.microrna.gr/ microT-CDS) 34,35 and microRNA.org 36 were used to predict the putative targets of differentially expressed 44 miRNAs. To improve the accuracy of the target prediction, we combined the microarray data of differentially expressed mRNA with the predicted targets of the differentially expressed 44 miRNAs. The intersecting gene set was subjected to the subsequent bioinformatic analysis.

Bioinformatic analysis.
To identify the functional pathways of intersecting genes and to uncover the miRNA-gene regulatory network on the basis of biological processes and molecular functions, we applied the Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.genome.ad.jp/kegg) enrichment analysis 37,38 . The list of the intersecting genes was matched onto the Database for Annotation, Visualisation and Integrated Discovery (DAVID) v6.7 39,40 , followed by the listing of the KEGG pathways with the p-value of the analysis. mirConnX. The web interface mirConnX 41 (http://www.benoslab.pitt.edu/mirconnx) was used to generate mRNA-miRNA interaction networks. The normalised miRNA and mRNA microarray expression data were used as input files. We selected mouse-mm9 (NCBI37)) _20110812 for organism type. Gene symbol and ID were selected for the gene ID and the miRNA ID, respectively. In the analysis options section, we selected Pearson's correlation as the association measure of choice. The regulation threshold for the minimum integrated regulation score was set at 0.5.