Investigating the genetic profile of dopaminergic neurons in the VTA in response to perinatal nicotine exposure using mRNA-miRNA analyses

Maternal smoking during pregnancy is associated with an increased risk of developmental, behavioral, and cognitive deficits. Nicotine, the primary addictive component in tobacco, has been shown to modulate changes in gene expression when exposure occurs during neurodevelopment. The ventral tegmental area (VTA) is believed to be central to the mechanism of addiction because of its involvement in the reward pathway. The purpose of this study was to build a genetic profile for dopamine (DA) neurons in the VTA and investigate the disruptions to the molecular pathways after perinatal nicotine exposure. Initially, we isolated the VTA from rat pups treated perinatally with either nicotine or saline (control) and collected DA neurons using fluorescent-activated cell sorting. Using microarray analysis, we profiled the differential expression of mRNAs and microRNAs from DA neurons in the VTA in order to explore potential points of regulation and enriched pathways following perinatal nicotine exposure. Furthermore, mechanisms of miRNA-mediated post-transcriptional regulation were investigated using predicted and validated miRNA-gene targets in order to demonstrate the role of miRNAs in the mesocorticolimbic DA pathway. This study provides insight into the genetic profile as well as biological pathways of DA neurons in the VTA of rats following perinatal nicotine exposure.


Results
The VTA was isolated from rat pups treated perinatally with nicotine or saline (control). After dissociating and sorting our samples, an average of 27.9% ± 3.3% of the NeuN-positive neurons were positively stained by TH. These results confirm the proportion of double-stained neurons reported by Guez-Barber et al., who found that 30.2% of all NeuN-positive events were double stained for TH and NeuN from midbrain samples. Another study by Chung et al. followed the same protocol and reported approximately 25% of NeuN-positive cells were also TH-positive 30 . The exact percentage of TH-positive neurons in the VTA is not known, but is speculated to be about 50-60% 31 ; therefore, we considered the estimated 25-30% of TH-positive neurons after fluorescent-activated cell sorting (FACS) to be adequate for the suggested experiment, consistent with other studies 30,32 . There was no significant difference in the number of collected DA neurons per VTA between the treatment groups since the saline-treated and nicotine-treated groups had 27.1% ± 3.7% and 28.5% ± 3.1% double-positive neurons from the total NeuN-positive neurons, respectively. Furthermore, statistical analysis using Student's t-test showed no significant difference (p > 0.5).
After collecting DA neurons using FACS, samples were processed using mRNA and miRNA microarrays in order to find the differential expression profiles of DA neurons after perinatal exposure to nicotine. The transcriptome and miRNome were then used to identify potential points of regulation and enriched pathways following perinatal nicotine exposure in VTA DA neurons. mRNA expression changes following perinatal nicotine exposure. For VTA DA neurons, 2,636 genes were found to be differentially expressed at a q-value < 0.05 using Benjamini & Hochberg (BH) method and an absolute log2 fold change >1. More specifically, 862 differentially expressed genes (DEGs) were upregulated and 1,774 DEGs were downregulated showing that the majority of genes in VTA DA neurons are downregulated due to nicotine exposure. Of these genes, 2,095 were annotated genes (639 upregulated and 1,456 downregulated). Figure 1 shows a heatmap of the top overall results for up and downregulated genes in DA neurons. Details about the top up and downregulated genes are shown in Table 1. miRNA expression changes following perinatal nicotine exposure. For DA neurons, there were 74 differentially expressed miRNAs (DEmiRs) of which 58 miRNAs were found to be upregulated, while 16 were downregulated showing that the majority of miRNAs were upregulated due to perinatal nicotine exposure (see Fig. 1 for heatmap). Applying parameters q-value < 0.05 (BH) and an absolute log2 fold change > 0. 5 29 , we found 11 DEmiRs that were upregulated and 9 DEmiRs that were downregulated. DEmiR results are shown in Table 2.
miRNA-mRNA target predictions. MultimiR 33,34 was used to find validated and predicted miRNA-gene targets using inversely regulated DEmiRs and DEGs. Considering the top 20% of predicted scores and conserved targets sites, 618 unique miRNA-gene target pairs were found with 290 nodes (58 miRNAs and 232 genes) based on conserved prediction sites. The majority of the miRNA-gene pairs were from upregulated DEmiRs and downregulated DEGs. Pair-wise Pearson correlation analysis was performed on expression values of DEGs and DEmiR. Multiple testing was corrected using BH method. Figure 2a shows the predicted network for miRNA-mRNA in VTA DA neuron after perinatal nicotine exposure. Notably, rno-miR-125a-5p was predicted to target the most number of unique genes (Fig. 2b).
Enriched pathway analysis and integrated miRNA-mRNA network. According to DAVID, pathway analysis of the downregulated DA DEGs revealed significant enrichment of many Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways associated with addiction, including neuroactive ligand-receptor interaction (p ≪ 0.001), calcium signaling (p ≪ 0.001), cAMP signaling (p < 0.01), and long-term potentiation (p < 0.05). Additionally, multiple signaling pathways and synapse pathways were enriched as well as drug addiction pathways (see Table 3 for more details). Notably, the DAergic synapse pathway was significantly enriched (p ≪ 0.001). Significant enrichment was determined using a modified hypergeometric test. Upregulated DA DEGs also showed non-significant enrichment of glutamatergic and serotonergic synapse pathways and neuroactive ligand-receptor interaction pathway.
A gene network was created using KEGGgraph 35 incorporating enriched pathways found using DAVID 36,37 . Additional relevant KEGG pathways were included from the literature 18,21 . A subset graph was created removing genes that were not significant in the microarray results in order to visualize gene-gene interactions within our results. Results from miRNA-gene target predictions were included in order to show the relationship between genes and pathways and their corresponding miRNA predictions. Figures 3 and 4 show the mRNA network created from the dopaminergic synapse pathway and neurotrophin signaling pathway, respectively, with miRNA targeting interaction. According to the hypergeometric test, we found two pathways with significant miRNA-gene target prediction with significant miRNA-pathway interaction results: DAergic synapse pathway and rno-miR-30b-5p was (p < 0.01) and the neurotrophin signaling pathway and rno-miR-195-5p and rno-miR-204-5p (p < 0.01).
Validation of miRNA and mRNA microarray results using RT-qPCR. The microarray results were validated for both miRNA and mRNA using the same approach as detailed in Bosch et al. 29 . Candidate genes were selected based on previously reported gene expression profiling of DA neurons in the VTA and biological interest from pathway enrichment and miRNA-gene target prediction. Significant miRNAs and mRNAs were selected for further analysis using RT-qPCR in order to validate the results of the microarrays. For mRNA microarray validation, two previously reported DA neuron markers were chosen, Cck and Gabrg2, as well as the significant DEGs found in the integrated miRNA-gene network, Scn1a, Ntrk2, and Ablim3. The direction of regulation resulting from RT-qPCR was consistent with the microarray results from our experiments, indicating that the data produced by the microarray experiments were valid. Regulation was determined using Student's t-test comparing nicotine and saline treatment groups and corrected for multiple testing using the BH method with false discovery rate of 0.05 (Fig. 5).

Discussion
The current study is the first to conduct a large-scale profile of both miRNA and mRNA expression in DA neurons located in the VTA of rats following perinatal nicotine exposure. In our current work, the offspring was exposed to nicotine over a period equivalent to the three trimesters of human pregnancy 4,13 . Our study sought to identify significant genes and miRNAs in DA neurons modulated by perinatal nicotine exposure. After differential expression analysis, our results showed 1,774 DEGs were downregulated and 862 DEGs were upregulated. Additionally, 16 DEmiRs were downregulated and 58 DEmiRs were upregulated. The gene expression profile of DA neurons was analyzed using DAVID for functional enrichment analysis. Enriched KEGG pathways were analyzed in order to summarize putative processes modulated by perinatal nicotine exposure compared to saline. In addition, target predictions for DEmiRs were compared with DEG results in order to identify miRNA-gene pairs that may be changed after perinatal nicotine exposure. Further analysis of these results revealed potential points of pathway regulation by DEmiRs that may be involved in altered mechanisms following perinatal nicotine exposure. Many gestational or perinatal nicotine exposure studies have focused on nicotinic acetylcholine receptors (nAChRs), which were reduced following gestational nicotine exposure in adolescent rats. One study exposed rats to nicotine during fetal and neonatal brain development (similar to the time of exposure used in our study) 13 . They found that during P7-14, the number of high-affinity nAChRs reach their highest levels in most brain regions 13 , which is the time period of our experiment. They also found decreased expression of multiple nAChR subunits including β4 mRNA subunit (Chrnb4), which we found significantly reduced as well. Additionally, Chrnb4 is involved in the neuroactive ligand-receptor interaction and cholinergic synapse pathways, which were significantly enriched by downregulated DEGs. The β4 subunit is only expressed in about 10% of DA neurons with a much higher expression in GABA neurons, but the time point of our study and equivalent amount of DA neuron total RNA may account for the expression pattern of nAChRs found in this experiment. Further studies focusing specifically on the nAChR subunits are necessary in order to determine if nAChR mRNA expression is decreased due to nicotine exposure or if reduced expression is caused by decreased DA neuron population.
Notably, the nicotine addiction KEGG pathway was enriched from our downregulated DEGs and four genes were found to be significant among our results. The enrichment of the nicotine addiction pathway indicates that there may be common genetic alterations in the brain that occur between gestational nicotine-treated offspring and adults with nicotine dependence. The first two genes, Gabrd and Gabrg2, are subunits of the γ-Aminobutyric acid (GABA) A receptor which were both significantly downregulated in our results. Our microarray tested the α-, β-, and θ-subunits of the GABA-A receptor, but these were not significantly differentially expressed in our results. A recent study by Stojakovic et al. investigated the selective deletion of Gabrg2 on DA neurons. Mice with the deletion showed a high risk of developing alcohol addiction based on behavior and an increased conditioned place preference to alcohol 38 . These results suggest that despite the expression of other subunits of the GABA-A receptor, Gabrg2 expression is important in DA neurons and DA neurotransmission in response to drugs. Another study investigating DA neurons in the VTA following cocaine exposure found GABA-A receptor subunits, including Gabrg2, to be significantly downregulated 17 . While this study did not use a model of gestational drug exposure, it focused on alterations within VTA DA neurons following chronic drug exposure.
The last two genes were downregulated in our results and encode subunits of glutamate receptors. Grin2d is the 2D subunit of N-methyl-D-aspartate (NMDA) receptors, which are a class of glutamate ionotropic receptors. Gria3 is subunit 3 of the alpha-amino-3-hydroxy-5-methyl-4-isoxazole propionate (AMPA) type glutamate ionotropic receptors. A study by Roguski et al. investigated intra-VTA NMDA receptors and concluded that a subpopulation of NMDA receptors were responsible for increased vulnerability to DA release in NAc. In addition, previous studies indicated that the expression of NMDA subunits was unchanged due to chronic ethanol exposure followed by withdrawal, but the function of the receptor was changed 11 . Further investigation needs to be conducted on all receptor subunits in order to determine if GABA and glutamate receptors on DA neurons are significantly modulated by perinatal nicotine exposure.
To date, a vast body of research has been conducted investigating the implication of miRNAs in addiction using different drugs of abuse. However, our present study investigates, for the first time, the large-scale effects of miRNAs in VTA DA neurons after perinatal nicotine exposure. We recently reported the upregulation of miR-140-5p after perinatal nicotine exposure. From our miRNA expression results, miR-140-3p and miR-140-5p were both significantly differentially expressed and have both been implicated in nicotine addiction 26 . Huang et al. found that miR-140-3p expression was increased by nicotine, which then repressed the expression of dynamin 1 (Dnm1) by direct base-pairing to the 3′-untranslated region of the gene 39 . In accordance with their results, miR-140-3p was significantly upregulated and Dnm1 was downregulated in our results, although not significantly. In our previous study, we investigated the modulation of addiction-related miRNAs by nicotine and found upregulation of miR-140-5p, which we confirmed with our current results 26 . Further studies should be conducted on  miR-140-3p and miR-140-5p and their potential gene targets in order to elucidate their relationship to nicotine exposure. The miRNA-gene network indicated that miR-125a-5p was the most highly predicted miRNA based on DEG targets (Fig. 2b). A study by Bosch and colleagues identified miR-125a-5p as an addiction-related miRNA in a study investigating the changes in RNA expression in the VTA due to the self-administration of methamphetamine in rats 29 . Their study found an overall upregulation of miR-125a-5p in the VTA 29 , which confirms the upregulation of miR-125a-5p in VTA DA neurons in our results. The significance of miR-125a-5p in the VTA in both of these studies implicates its role in regulating pathways related to addictive drugs, including methamphetamine and nicotine. Further study into the genes targeted by miR-125a-5p could provide a clearer picture of the precise nature of this miRNA in pathways related to substance abuse.
Our integrated network of DEGs and DEmiRs was created from the enriched KEGG pathways with the addition of significantly upregulated predicted miRNA in order to visualize putative points of pathway regulation by miRNAs. These pathways were enriched by the downregulated DEGs indicating important putative pathways for miRNA regulation as tested by hypergeometric testing with correction for multiple testing. Significantly upregulated predicted miRNAs were added to this network in order to visualize points of pathway regulation by miRNA. In the DAergic synapse pathway, miR-30b-5p was predicted to target Scn1a, a protein-coding gene involved in voltage-gated sodium channels (Fig. 3). As a whole interaction network, the probability of the involvement of miR-30b-5p in the DAergic synapse pathway was significant (p < 0.01). The enrichment of this pathway supports the hypothesis that the mesocorticolimbic DA pathway is involved in the reinforcing effects of substance abuse. The DAergic synapse pathway describes the release of DA neurotransmitter from DA neurons to a postsynaptic neuron. According to the major hypothesis of drug reinforcement, the reinforcing effect of addiction is believed to be conveyed through the activation of the mesocorticolimbic DA system. Stimulation of VTA DA neurons via nicotine administration results in the release of DA in the NAc, describing the role of the DA synapse pathway in the reinforcing effect 7 . The enrichment of this pathway in our results demonstrates that perinatal nicotine exposure leads to genetic alterations in VTA DA neurons that are in accordance with addiction mechanisms. In addition, we have highlighted a potential gene-miRNA target interaction, which may contribute to the observed alterations via post-transcriptional gene regulation. The neurotrophin signaling was enriched in our results at a significance p < 0.1. Our analysis revealed that Ntrk2 (also known as TrkB) was predicted to be targeted by miR-204-5p and is involved in the neurotrophin signaling pathway and miR-195-5p was predicted to target Ikbkb (Fig. 4). As a whole interaction network, the probability of the involvement of both miRNAs in the neurotrophin signaling pathway was significant (p < 0.01). Many studies have been conducted linking nicotine and the neurotrophin signaling pathway, but do not include the VTA DA neurons 28 . TrkB, a transmembrane receptor with a high affinity for tyrosine kinaseB, is activated by Bdnf (brain-derived neurotrophic factor), which has been reported to potentiate the effects of addictive substances through the mesocorticolimbic DA pathway 40,41 . Additionally, TrkB has been identified as a susceptibility gene for psychiatric disorders, such as schizophrenia and other mood and anxiety disorders. TrkB plays a role in synaptic plasticity and neurotransmitter release and may provide additional support for the alterations induced by nicotine in DA neurons in the VTA 42 .
In summary, we focused on the response of VTA DA neurons to perinatal nicotine exposure and investigated transcriptional and post-transcriptional regulation disruptions. We also concentrated on the identification of biological pathways modulated by nicotine within the DA neurons of the VTA. Our study suggested that the dopaminergic synapse pathway was significantly altered by perinatal nicotine exposure as well as significant interactions with miRNAs. Although our results indicate that perinatal nicotine exposure alters the expression of miRNAs and genes, demonstrating the involvement of several biological pathways in the mechanism of addiction in DA neurons of the VTA, we did not investigate the effect of gender on DA neurons in the VTA following perinatal nicotine exposure in the current study. Such an investigation merits further exploration as there is evidence indicating nicotine exposure differentially affects gene expression changes in many different regions of the brain depending on gender 43 . Additionally, further investigation needs to be performed to develop an interactive model to evaluate the effect of miRNAs on biological pathways. Our results predicted that Ntrk2 was targeted by miR-204-5p (dark blue). Additionally, rno-miR-195-5p (dark blue) was predicted to target Ikbkb. The predicted interaction of both miRNAs on the neurotrophin signaling pathway was significant (p < 0.01).Light blue nodes are significantly downregulated DEGs in the results and yellow nodes are genes that were not detected or not significantly in the results. In addition to significance, location of the mRNA has been indicated by color (pink for membrane, orange for nucleus, otherwise cytoplasm).

Materials and Methods
Animal treatment. All experiments were performed in accordance with the protocols and surgical procedures that were approved by the Institutional Animal Care and Use Committee (IACUC) and the University of Houston Animal Care Operations (ACO). Pregnant female Sprague-Dawley (SD) rats (Charles River, Wilmington, MA, USA) were maintained on a 12-h light/12-h dark schedule at a temperature of 22 ± 2 °C and 65% humidity. Access to standard food and water was ad libitum. Rats were acclimated to the animal facility for 72 hours before they received treatment via an osmotic minipump (Alzet, Cupertino, CA) that was implanted subcutaneously containing either nicotine hydrogen tartrate (Sigma-Aldrich, St. Louis, MO, USA) released at a rate of 6 mg/kg/day (moderate to heavy smokers 5,12 ), or an equal volume of saline vehicle for the control. Nicotine continues to be released from the osmotic minipump for 4 weeks from gestational day 6 to postnatal day 14.
Seven-to-fourteen day-old pups (male and female) were anesthetized with isoflurane before decapitation. On a VT1200 semiautomatic vibrating blade microtome (Leica, Nussloch, Eisfeld, Germany), horizontal slices containing VTA were cut at a thickness of 1000 μm. Brain punches containing the VTA were collected bilaterally using a 1 mm biopsy punch (Integra Miltex, VWR, Radnor, PA, USA) and placed in 1 mL of Hibernate A (Gibco, Thermo Fisher Scientific, USA) on ice to maintain cell viability. Brain punches were pooled so that one litter resulted in one litter. A total of four samples for both saline and nicotine treated groups were processed for RNA extraction and microarray processing.
Cell Dissociation, FACS, and RNA extraction. Tissue was dissociated into a single cell solution and sorted by FACS as reported by Guez-Barber et al. 32 . Briefly, tissue punches were placed in 1 mL of Accutase (Gibco, Thermo Fisher Scientific, USA) and shaken for 30 minutes at 4 °C in order to dissociate the tissue punches. Accutase is a mixture of proteolytic and collagenolytic enzymes, which showed to produce the most single cells with the least cell damage 32 . Cells were pelleted at 425 × g and resuspended in Hibernate A. To further dissociate cells, cells were gently pipetted with increasingly smaller pipette tips. The supernatant containing single cells was collected until all cells were collected. Then, to remove debris and cell clusters, the cell suspension was serially filtered through pre-wetted 100 µm and 40 µm cell strainers. The strained cell suspension was added to a three-density step gradient made using Percoll (GE Healthcare, VWR, USA) and centrifuged at 430 × g for 3 minutes in order to further remove debris. The cloudy top layer containing debris was removed and remaining solution was centrifuged at 550 × g for 5 minutes to pellet the cells. The cells were fixed for immunolabeling by resuspending them in equal parts Hibernate A and 100% cold ethanol, gently vortexed, and kept on ice for 15 minutes. Cells were co-labeled with conjugated primary antibodies neuronal marker, NeuN/Alexa Fluor 488 (NeuN/ AF488, ab190195, Abcam, Cambridge, MA, USA), and tyrosine hydroxylase/phycoerythrin (TH/PE, ab209921, Abcam, Cambridge, MA, USA). Using conjugated antibodies simplifies the staining procedure by eliminating the need for a secondary antibody, which can bind unspecifically and reduces the time it takes to stain cells.
An Influx (BD Biosciences, San Jose, CA, USA) instrument was used to sort the cells at the Flow Cytometry and Cellular Imaging Core Facility (MD Anderson -South Campus, Houston, TX, USA). Samples were sorted based on double-positive NeuN + /TH + staining. Once sorted, cells were centrifuged into a pellet at 2650 × g for 8 minutes at 18 °C. Total RNA was isolated using miRNeasy Micro Kit (Qiagen, Hilden, Germany) following manufacturer's instructions, including DNAse treatment and used in subsequent microarray and RT-qPCR validation experiments. RNA quality and quantity was established according to the optical density (OD) of each sample at 260 nm and 280 nm determined using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). Only total RNA samples with A260/280 ratio of 1.9 or greater were used in subsequent experiments. mRNA and miRNA expression microarrays. All gene and miRNA expression reagents and kits were purchased from Agilent (Santa Clara, CA, USA) unless otherwise stated. The gene expression was profiled using a SurePrint G3 Rat Gene Expression v2 8 × 60 K microarray (ID: 074036) with 30,584 unique genes. With 25 ng of total RNA, samples were prepared using the One-Color Low Input Quick Amp Labeling kit with RNA Spike-Ins according to manufacturer's instructions. First, total RNA was amplified and labeled with Cyanine-3 (Cy3). The amplified cRNA was purified and quantified using a NanoDrop 2000 spectrophotometer. Next, 600 ng of Cy 3-labeled cRNA with required specific activity (≥6 pmol Cy3/µg cRNA) were fragmented and prepared for hybridization using the Gene Expression Hybridization kit. Then, slides were hybridized for 17 hours at 65 °C and washed in Gene Expression Wash Buffers according to manufacturer's protocol.
For miRNA expression, an 8 × 15 K Rat miRNA Microarray, Release 21.0 (ID: 070154) was used containing 758 mature miRNAs. For labeling and hybridization, miRNA Complete Labeling and Hyb kit with RNA Spike-Ins was used according to manufacturer's instructions. RNA Spike-In was included to test quality control of the microarray amplification, labeling, and hybridization. Briefly, 100 ng of total RNA containing miRNAs was dephosphorylated and labeled with Cy3-pCp. Samples were purified on a Micro Bio-Spin P-6 gel column (Bio-Rad), dried, and hybridized at 55 °C for 20 hours. Microarray slides were washed using Gene Expression Wash Buffers according to manufacturer's protocol.
The gene and miRNA expression slides were scanned using G4900DA SureScan Microarray Scanner using G3_GX_1color and AgilentG3_miRNA scan protocols, respectively. Microarray data was extracted from the TIFF result images using Feature Extraction (FE) Software v12.0.1.1 and the corresponding protocol, GE1_1200_ Jun14 FE protocol for gene expression and miRNA_1200_Jun14 protocol for miRNA expression.
Data Analysis. All pre-processing, normalization, and statistical analyses were performed using several Bioconductor packages in R version 3.4.2 44 . The quality of microarray data was assessed using arrayQualityMetrics 45 package. Raw mRNA expression data was imported using limma 46 package. Before normalization, probes expressing Agilent flags for feature and background outliers were removed from the analysis. Raw median intensity values were background corrected using "normexp" method and quantile normalized. Control probes were filtered out as well as low expressed probes, which were defined by expression intensity less than 10% brighter than 95% intensity of negative controls. Replicated probes were averaged and a total of 27,791 genes were selected for analysis. Quality assessment was performed to check for outliers among microarrays. A linear fit model was used to calculate fold changes and standard errors for each gene of interest. To increase the power of the data, standard errors were moderated using the eBayes function, which applies a simple empirical Bayes model and computes a log-odds of differential expression for each contrast. P-values were corrected for multiple testing using BH method. DEGs were identified using a threshold of q < 0.01 and absolute log fold change >1.
For miRNA microarray data, raw intensity data was imported and processed using AgiMicroRna 47 package. Quality of the microarray was assessed using arrayQualityMetrics package as well as quality control functions within the AgiMicroRna package. Raw miRNA data was preprocessed using the robust multi-array average (RMA) algorithm without background correction. Data was filtered according to the following criteria: (1) positive for flag IsGeneDetected, (2) expressed in at least 50% of the experimental samples, and (3) signal intensity greater than the mean value of the negative control + 1.5 standard deviations. After filtering, 329 miRNAs remained for analysis. Linear model was fitted to the miRNA expression data and moderated statistics calculated using eBayes, similar to mRNA data analysis. Differential expression was identified using q-value threshold <0.05 (calculated using BH method) 28,29 . Integrated analysis of miRNA-mRNA. Functional enrichment analysis was performed on differentially expressed genes using DAVID v6.8 36,37 . Enriched KEGG pathways were identified and analyzed. Further analysis of KEGG pathways was done using DEgraph 48 and KEGGgraph 35 . MultiMiR 33,34 was used to identify predicted and validated miRNA-gene target pairs based on inversely correlated regulation of DEmiRs and DEGs. multiMiR is a comprehensive compilation of predicted and validated miRNA-gene target interactions from 14 external databases. Hypergeometic testing was performed on target DEGs negatively correlated with DEmiRs using pairwise Pearson correlation analysis and then BH corrected for multiple comparisons 49 . Remaining significant miRNA-gene target pairs were integrated into gene network models based on enriched KEGG pathways. Following the procedure used in the miRPathDB 50 , we performed a hypergeometric test with multiple testing correction using BH method and an FDR < 0.05 for each miRNA-gene predicted interaction. miRPathDB is a tool to measure the association between miRNAs and putative target pathways. Quantitative RT-PCR validation of microarray data. All reagents and kits for quantitative reverse transcription polymerase chain reaction (RT-qPCR) were purchased from Applied Biosystems (Thermo Fisher Scientific, Carlsbad, CA, USA) unless otherwise stated. Seven genes and five miRNAs were chosen for validating microarray data using RT-qPCR. Total RNA was isolated from FACS samples for each experimental group as described above. For gene expression validation, cDNA was prepared using High Capacity cDNA Reverse Transcription Kit according to manufacturer's instructions. For miRNA expression validation, cDNA was prepared and preamplified using TaqMan Advanced miRNA cDNA Synthesis Kit according to manufacturer's instructions. For gene validation and miRNA validation, quantitative PCR (qPCR) was carried out using TaqMan Fast Advanced Master Mix and corresponding TaqMan Assay (TaqMan Gene Expression Assay or TaqMan Advanced miRNA Assay) on a StepOnePlus Real-Time PCR System according to manufacturer's instructions using the following parameters: 2 min at 50 °C, 2 min at 95 °C, 40 cycles of 1 sec at 95 °C and 20 sec at 60 °C. Each reaction was prepared in triplicate for both validation sets. Comparative Ct method was used to find the relative quantity of the target genes or miRNAs. Student's t-test was performed comparing nicotine and saline treatment groups (n = 3) and corrected for multiple testing using the BH method with false discovery rate of 0.05.