Comparison between dopaminergic and non-dopaminergic neurons in the VTA following chronic nicotine exposure during pregnancy

Exposure to nicotine during pregnancy through maternal smoking or nicotine replacement therapy is associated with adverse birth outcomes as well as several cognitive and neurobehavioral deficits. Several studies have shown that nicotine produces long-lasting effects on gene expression within many brain regions, including the ventral tegmental area (VTA), which is the origin of dopaminergic neurons and the dopamine reward pathway. Using a well-established rat model for perinatal nicotine exposure, we sought to investigate altered biological pathways using mRNA and miRNA expression profiles of dopaminergic (DA) and non-dopaminergic (non-DA) neurons in this highly-valuable area. Putative miRNA-gene target interactions were assessed as well as miRNA-pathway interactions. Our results indicate that extracellular matrix (ECM) receptor interactions were significantly altered in DA and non-DA neurons due to chronic nicotine exposure during pregnancy. They also show that the PI3K/AKT signaling pathway was enriched in DA neurons with multiple significant miRNA-gene targets, but the same changes were not seen in non-DA neurons. We speculate that nicotine exposure during pregnancy could differentially affect the gene expression of DA and non-DA neurons in the VTA.

Maternal smoking during pregnancy is associated with adverse birth outcomes including cognitive and neurobehavioral deficits 1-3 . The results from several studies have suggested that the exposure to nicotine, the major addictive and psychoactive component of tobacco, from nicotine replacement therapy may carry risks to the offspring similar to maternal smoking during pregnancy. To investigate these underlying mechanisms, we used an established rat model of perinatal nicotine exposure using a moderate to heavy dose of nicotine 4,5 . To date, several studies have used this animal model to elucidate the genetic, biochemical, behavioral, and electrophysiological alterations that occur in response to nicotine exposure at different time points of development in different regions of the brain and by gender 3,6,7 . Due to the complexity of these contrasts, we have specifically focused on the category of neuron type: dopaminergic (DA) and non-dopaminergic (non-DA) neurons in the VTA following perinatal nicotine exposure without separating for gender differences in order to investigate more universal alterations.
Due to the addictive nature of nicotine, the brain regions of interest are part of the mesocorticolimbic pathway also known as the reward pathway. This pathway is a network of dopaminergic (DA) neurons connecting the ventral tegmental area (VTA) to the striatum, nucleus accumbens (NAc), and prefrontal cortex (PFC). In the VTA, there are three main neuron types: DA, GABA, and glutamate, which interact with each other to inhibit and activate the release of DA in other regions of the brain. Under normal behavior, DA neurons are inhibited by input from GABA neurons and activated by input from glutamate neurons. After perinatal nicotine exposure, this relationship is altered, as shown by nicotine-induced decreases in dopamine release in the NAc in adolescent rats 8 . Nicotine modulates nicotinic acetylcholine receptors (nAChRs), which are ligand-gated ion channels distributed through the central and peripheral nervous system. In the VTA, nicotine directly activates DA neurons through these receptors as well as indirectly via glutamate and GABA neurons [9][10][11] . Following perinatal nicotine exposure, studies have found decreased expression of nAChR subunits in the VTA as well as a decrease in the number of DA neurons 5,12,13 . These findings demonstrate specific changes that occur during development in response to chronic nicotine exposure during development. While previous studies investigating perinatal nicotine exposure have investigated biological pathways that are modulated by nicotine, none to date have explored the VTA 14 While extensive research has been conducted on epigenetic, genetic, and behavioral changes in response to perinatal nicotine exposure, less information exists on alterations that occur in the miRNome. Although there are several studies that have investigated microRNA (miRNA) changes in response to other drugs of abuse including cocaine and methamphetamine [17][18][19] , there are very few studies investigating miRNome expression in response to chronic nicotine exposure during pregnancy. Therefore, we previously conducted a small-scale study on specific addiction related miRNAs in DA and ND neurons in the VTA 20 . MicroRNAs (miRNAs) are short, non-coding RNA sequencing spanning about 22 nucleotides that act to regulate the expression of mRNAs by binding to the 3′ untranslated region (3′UTR) which results in translational inhibition or transcriptional repression 17,21 . They are of particular interest because a single miRNA can regulate many factors affecting gene regulatory networks, thus affecting a large downstream response 21,22 . MiRNAs have been found to modulate the effects of abusive substances and to be involved in synaptic plasticity, developing connections, and other behaviors related to addition 17,23,24 .
In the present study, we investigated RNA expression profiles of DA and non-DA neurons in the VTA following perinatal nicotine exposure due to the importance of the VTA in the mesocorticolimbic DA reward pathway as well as the changes in miRNA and mRNA expression that are modulated by nicotine. We determined significantly differentially expressed mRNAs and miRNAs as well as significant miRNA-mRNA target interactions in order to identify potential gene regulatory networks and biological pathways that result from perinatal nicotine exposure in DA and non-DA neurons in the VTA.

Results
Nicotine-treated DA and non-DA neurons were compared on the basis of transcriptome and miRNome profiles in order to find and compare the differential expression between the two neuron groups after perinatal exposure to nicotine. Rat offspring were treated with nicotine or saline (control) for up to 28 days from G6 to the second postnatal week. This period of exposure in the rat model is developmentally equivalent to the three trimesters of human pregnancy 3,25 . The VTA was then isolated from the offspring and separated into DA and non-DA neuron groups based on expression of TH (a common dopaminergic neuron marker) and NeuN (a neuronal marker) using fluorescent-activated cell sorting (FACs) (Fig. 1). Additional detail about FACs can be found in Supplemental Figure 1. The ratio of DA and non-DA neurons we collected was confirmed by reported yields from Guez-Barber et al. 26 and Chung et al. 27 .
To investigate the differential expression profile between DA and non-DA neurons after perinatal nicotine or saline (control) exposure, we used Agilent's SurePrint microarrays to compare gene and miRNA expression profiles. Figure 1b Between the DA and non-DA neuron groups, 33 significant DEGs were regulated in the same direction. These genes may indicate similar changes that occur in all cells in the VTA in response to nicotine exposure during pregnancy. Additionally, our data found 45 significant DEGs that were oppositely regulated when comparing the DA and non-DA neuron group results (Table 1).

Functional enrichment analysis of DEG lists.
Enriched GO biological processes were investigated using the Cytoscape v3.6.0 28 (http://cytoscape.org/) plugin ClueGO v2.5.1 29 (Fig. 2). For the DA neuron group, the most enriched biological processes were axon development and neuron ensheathment (p ≪ 0.0001) for the downregulated DEGs and ion transport (p < 0.01) for the upregulated DEGs. For upregulated non-DA neurons, the generation of neurons was the most enriched biological process (p ≪ 0.0001), while for downregulated non-DA neurons, synapse organization (p < 0.001) was the most significantly altered biological process. In addition, we found that the expression of up-regulated DEGs in the non-DA neurons showed more dramatic changes than the changes in the DA neurons.
Next, we analyzed the DEG lists using SPIA, which compared the log fold change results from our DA and non-DA groups with signaling pathway topology in order to find significantly altered KEGG pathways 30 (see Table 2). Two pathways were implicated in both DA and non-DA neuron results: the ECM-receptor interaction pathway was inhibited by DA DEGs and activated by non-DA DEGs; while the GABAergic synapse pathway was inhibited in both groups. Notably, the results from non-DA group showed the morphine and amphetamine addiction pathways to be inhibited.
Differential miRNA expression analysis of dopamine and non-dopamine neurons following perinatal nicotine treatment. Applying a similar differential analysis scheme for miRNA expression analysis, a total of 329 annotated mature miRNAs were tested on the microarray for differential expression for each cell type. Using a q-value <0.05, we found 58 upregulated differentially expressed miRNAs (DEmiRs) and 16  DEmiRs for non-DA neurons. Interestingly, DEmiRs were more dysregulated in the DA group compared to the non-DA group, while the opposite interpretation was observed for DEGs.

Integrated analysis of DEmiRs and DEGs to find enriched pathways.
We processed the lists of DEmiRs and DEGs for DA and non-DA in order to find putative miRNA-gene pairs prevalent to perinatal nicotine exposure. Using multiMiR, DEmiRs and DEGs with opposite regulation were paired and compared to several databases with validated and predicted miRNA-gene target interaction scores. A total of 644 and 47 putative miRNA-gene pairs with opposite regulation were found for DA and non-DA neurons, respectively. In order to further assess target interactions and to eliminate weak correlations, pair-wise Pearson correlation was applied with BH correction. Figure 3 shows the results from our analysis, which found 125 and 42 unique miRNA-gene pairs for DA and non-DA neurons, respectively. Among the remaining significant DEmiRs, rno-miR-29b-3p and rno-miR-30b-5p were significantly enriched in both groups. Using miRPath v.3, enriched pathways were found from the results of DEmiR-DEG pairs using the input DEmiRs and DEGs lists. Table 3 shows the pathway enrichment results from DEmiR-DEG network for DA and non-DA neurons. Due to the significance of the ECM-receptor interaction pathway in both the DA and non-DA groups, we visualized the network with putative miRNA interaction (Fig. 4). MiRNA interaction with the ECM-receptor interaction pathway was p ≪ 0.0001 for ND DEmiRs. For the DA neuron, the PI3K/AKT signaling pathway was significantly enriched by DEGs (p < 0.1) with significant miRNA interaction (p ≪ 0.0001) (Fig. 5). Validation of microarray data by RT-qPCR. The microarray results were validated for both miRNA and mRNA using the same approach detailed in Keller et al. 31 . Genes and miRNAs were selected based on validated miRNA-gene target interaction identified using multiMiR as well as biological interest from functional enrichment analysis. For validation of the microarrays, we chose 14 mRNAs and 10 miRNAs. All mRNAs and miRNAs tested by RT-qPCR showed the same direction of regulation, which confirmed our microarray results. The results from RT-qPCR validation are shown in Fig. 6. Reference samples used were saline-control and reference was GAPDH for mRNA testing and U6 snRNA for miRNA testing. Significance was evaluated using Student's t-test (n = 3) and corrected for multiple comparisons using BH procedure with false discovery rate of 0.05.

Discussion
In this study, we analyzed the differential expression of mRNA and miRNA from VTA neurons in order to compare the function enrichment analyses of DA and non-DA neurons from the VTA with respect to perinatal nicotine treatment. To investigate the effect of chronic nicotine exposure during pregnancy, we used a well-established animal model for perinatal nicotine exposure that mimics maternal smoking or nicotine replacement therapy during a full term human pregnancy 25 . We employed FACS in order to sort VTA neurons based on the expression of TH, a DA neuron marker, and NeuN, a neuronal marker, and confirmed our results against previous studies 26,27 . In summary, we found 711 downregulated and 183 upregulated DEGs for DA neurons and 925 downregulated and 499 upregulated DEGs for non-DA neurons based on q-value < 0.01 (BH-corrected) and absolute log2 fold change >1. For miRNA differential expression, we found 16 downregulated and 58 upregulated DEmiRs for DA neurons and 4 upregulated and 4 downregulated DEmiRs for non-DA neurons based on q-value <0.05 (BH-corrected). Enriched GO biological processes were analyzed to expose altered pathways in response to perinatal nicotine exposure and to reveal how nicotine alters RNA expression differently based on neuron type. Lastly, we visualized putative miRNA-gene target interactions that may highlight gene regulatory networks that are significantly modulated by exposure to nicotine during development.
Our functional enrichment analysis of DA and non-DA DEGs revealed multiple GO biological processes related to cell death and survival pathways to be significantly enriched, results that are consistent with a study by Wei et al., who investigated the apoptotic effects of perinatal nicotine using a pathway-focused microarray. In that study, which investigated only a panel of 29 genes using RT-qPCR, they concluded that, in adolescents, perinatal nicotine modulates gene expression in cell death and survival-related pathways, could be classified into three subgroups: growth factor, death receptor, and kinase cascade, and these alterations differed by brain region 14 33 . Therefore, significant downregulation of genes involved in neuron ensheathment may result in dysregulated behavior along the DA reward pathway. The upregulated DEGs in the DA neuron were most significantly enriched in the ion transport process, which was a significantly enriched process when analyzing a list of nicotine addiction-related genes compiled by Liu et al. 34,35 . This implicates that the genetic alterations found in nicotine addiction studies relating to ion transport may be more significant in DA neurons of the VTA. It was found that synapse organization was significantly altered by downregulated non-DA DEGs, which is highly important in the arrangement and microenvironment of neurons during neurodevelopment. The third trimester of human gestation and the first twelve postnatal days for rats are characterized by rapid brain growth and is a time period when nAChRs are highly involved in the formation of cortical and limbic circuitry that is essential for synaptic formation and neuronal signaling 25 . Furthermore, upregulated non-DA DEGs were highly implicated in many processes related to the generation of neurons, including morphogenesis, neurogenesis, development, and differentiation. Neuronal proliferation and death are intrinsic to normal neurodevelopment; therefore any alteration of the balance between proliferation and death may have severe results 14 . Perinatal nicotine exposure has been shown to decrease the total number of neurons in the brain and certain sub-regions including the VTA 7,12 .
Results from our miRNA analysis suggest that miRNAs are more involved in dysregulation of DA neurons compared to non-DA neurons in the VTA. The common miRNA, rno-miR-30b-5p, was implicated in several pathways, including mucin type O-glycan biosynthesis, gap junction, cGMP-PKG signaling pathway, and glutamatergic synapse pathways.
For the DA neuron, the PI3K/AKT signaling pathway was significantly enriched by DEGs with five putative DEmiR-DEG interactions. This pathway is involved in regulating cell proliferation and growth, and inhibiting cell death. Nicotine is known to trigger specific immune-related signaling pathways including PI3K/AKT signaling, which is manifested in the neuroprotective behavior attributed to nicotine 14,16,36 . Cui et al. found that immune pathways including PI3K/AKT signaling were modulated via Chrna7, the alpha-7 subunit of nAChR and exhibited neuroprotection against Alzheimer's disease in animal models 16 . The putative miRNA-gene targets may shed light on possible mechanisms to induce neuroprotection against degenerative disorders like Alzheimer's and Parkinson's diseases. Furthermore, other abusive substances such as cocaine have been shown to modulate the PI3K/AKT signaling pathway, showing similar characteristics between perinatal nicotine exposure and drug addiction 37 .
In the pathway analysis, the ECM-receptor interaction pathway was implicated in both DA and non-DA groups, which highlights the importance of the pathway with respect to nicotine exposure and how the variation of different genes can cause similar patterns of alteration. ECM receptors are important during neurodevelopment-playing a key role in axon guidance, dendrite projection, synapse formation, and connection establishment between cells 38 . Many of these processes were significantly enriched in our analysis of GO biological processes. In summary, we investigated significantly altered biological pathways determined by DEGs and DEmiRs. Analyzing the relationships between miRNAs and genes revealed putative miRNA-gene target interactions that may be important in regulation of biological pathways in DA and non-DA neurons in the VTA. Our results suggest that ECM-receptor interactions are significantly altered in DA and non-DA neurons due to chronic nicotine exposure during pregnancy. Additionally, the PI3K/AKT signaling pathway was enriched in DA neurons with multiple significant miRNA-gene targets, but the same changes were not seen in non-DA neurons. Further investigation of the miRNA-gene target interaction pairs needs to be conducted. The development of an interactive model to evaluate the role of these miRNAs from our predictive model in regulating gene networks and biological pathways is also necessary.

Materials and Methods
Animal treatment. All experiments were performed in accordance with the protocols 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 in 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 implanting a subcutaneous osmotic minipump (Alzet, Cupertino, CA) containing either nicotine hydrogen tartrate salt (Sigma-Aldrich, St. Louis, MO, USA) released at a rate of 6 mg/kg/day (moderate to heavy smoker), or an equal volume of saline for the control. Over a period spanning from gestational day 6 (G6) to delivery (around G21-22), the offspring received nicotine via the placenta, which readily allows nicotine to cross the placental barrier 4,39 . After birth, the offspring continue to receive nicotine via the milk which has a milk-to-plasma ratio of 2.9 ratio 4 .
Ten-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), 1 mm horizontal slices containing VTA were sliced and 1 mm biopsy punch (Integra Miltex, VWR, Radnor, PA, USA) was used to collect the VTA bilaterally. The brain punches were placed on ice in Hibernate A (Gibco, Thermo Fisher Scientific, USA) to maintain tissue viability. Brain punches from eight to ten pups from one litter were pooled into one sample with an equal number of male and female pups. We collected and analyzed four samples for each treatment group (saline and nicotine) for RNA extraction and microarray processing.  Cell Sorting and RNA extraction. Brain punches were dissociated and sorted by FACS as reported by Guez-Barber et al. 26 . Briefly, brain punches were dissociated in Accutase (Gibco, Thermo Fisher Scientific, USA). Cellular debris was removed by serial filtration through decreasing cell strainers and then density centrifugation   Microarray Preparation, Labeling, and Hybridization. All gene and miRNA microarray 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 using a starting quantity of 25 ng total RNA. Samples were prepared using the One-Color Low Input Quick Amp Labeling kit with RNA Spike-Ins according to manufacturer's instructions. Then, amplified cRNA was fragmented and prepared for hybridization using the Gene Expression Hybridization kit according to manufacturer's instructions. Then, the slides were hybridized for 17 hours at 65 °C.
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 using starting quantity of 100 ng total RNA containing miRNAs. Samples were purified on a Micro Bio-Spin P-6 gel column (Bio-Rad), dried, and hybridized at 55 °C for 20 hours. After hybridization, the 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 appropriate protocols. Raw microarray data was extracted from microarray images using Feature Extraction (FE) Software v12.0.1. Data Analysis. Data analysis was performed as described previously in Keller et al. 31 . All pre-processing, normalization, and statistical analyses were performed in R version 3.4.2 40 using the following packages: arrayQuali-tyMetrics 41 to assess quality of raw microarray data, limma 42 to analyze mRNA microarrays, and AgiMicroRna 43 to analyze miRNA microarrays. For the mRNA microarray experiment, feature and background outliers were removed from the data based on Agilent flags. Then, the data was background corrected using "normexp" method and quantile normalized. For the miRNA microarray experiment, raw data was preprocessed using the robust multi-array average (RMA) algorithm without background correction 44 . For mRNA and miRNA analysis, differential expression was found using 2 × 2 factorial design comparing treatment (nicotine vs. saline) and cell type (DA vs. non-DA). A linear model was fitted with a factor containing the four conditions (nicotine treated DA = NDA, saline treated DA = SDA, nicotine treated non-DA = NND, and saline treated non-DA = SND). Differential expression between two conditions was calculated using a moderated t-test. Correction for multiple testing for all microarray analyses was performed using Benjamini-Hochberg correction using false discovery rate of 0.05.
Functional enrichment analysis of DEGs. Functional enrichment analysis was performed on differentially expressed genes using DAVID v6.8 45,46 and the ClueGO v2.5.1 29 plugin for Cytoscape v3.6 28 . Gene ontology (GO) term and Kyoto Encylopedia of Genes and Genomes (KEGG) pathways were both included in the analysis. ClueGO performed single cluster analysis to create network of functionally-related terms biological processes using the list of significantly differentially expressed genes. To find significance, the p-value was calculated by performing a two-sided hypergeometric test and corrected for multiple testing using Benjamini-Hochberg method.
In order to identify enriched pathways, we used SPIA, which identifies the most relevant KEGG pathways by comparing a genelist of DEGs and their log fold change with signaling pathway topology 30 . Results from this analysis were corrected for multiple testing using the BH method. miRNA-mRNA integrated analysis. As previously described in Keller et al. 31 , multiMiR 47,48 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. For each predicted miRNA-gene pair, we performed pair-wise Pearson correlation analysis on expression levels of the gene and miRNA in order to assess probability of the significance. We removed pairs that did not meet the following parameters: r < −0.5 and q < 0.05, according to Mach et al. 49 . From these results we compiled a miRNA-gene network containing all DEmiRs and DEGs with significant correlation. The miRPath v.3 50 tool was used to perform functional enrichment analysis on the miRNA-gene network with BH correction for multiple comparisons. Figure 6. Validation of (a,b) mRNA and (c,d) miRNA microarray results by RT-qPCR. Fourteen DEGs and 10 DEmiRs were tested for DA and non-DA to assess the validity of microarray experiments. The results are shown as ΔΔCt values relative to saline-counterpart and reference primer, GAPDH for mRNA and U6 snRNA for miRNA. Significance was evaluated using Student's t-test (n = 3) and corrected for multiple comparisons using Benjamini-Hochberg procedure with false discovery rate of 0.05 (*p < 0.05, **p < 0.01, ***p < 0.001).
Scientific REPORTS | (2019) 9:445 | DOI:10.1038/s41598-018-37098-1 Quantitative RT-qPCR 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. As previously described in Keller et al. 31 , 14 genes and 10 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 on a StepOnePlus Real-Time PCR System 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. For each sample, there were three biological samples (n = 3) and each reaction was prepared in triplicate. Comparative Ct method was used to find the relative quantity of the target genes or miRNAs compared to saline counterpart. GAPDH and U6snRNA were used as the reference gene and miRNA, respectively. Significance was computed using Student's t-test with BH-correction for multiple testing with false discovery rate of 0.05.