Investigating the influence of perinatal nicotine and alcohol exposure on the genetic profiles of dopaminergic neurons in the VTA using miRNA–mRNA analysis

Nicotine and alcohol are two of the most commonly used and abused recreational drugs, are often used simultaneously, and have been linked to significant health hazards. Furthermore, patients diagnosed with dependence on one drug are highly likely to be dependent on the other. Several studies have shown the effects of each drug independently on gene expression within many brain regions, including the ventral tegmental area (VTA). Dopaminergic (DA) neurons of the dopamine reward pathway originate from the VTA, which is believed to be central to the mechanism of addiction and drug reinforcement. Using a well-established rat model for both nicotine and alcohol perinatal exposure, we investigated miRNA and mRNA expression of dopaminergic (DA) neurons of the VTA in rat pups following perinatal alcohol and joint nicotine–alcohol exposure. Microarray analysis was then used to profile the differential expression of both miRNAs and mRNAs from DA neurons of each treatment group to further explore the altered genes and related biological pathways modulated. Predicted and validated miRNA-gene target pairs were analyzed to further understand the roles of miRNAs within these networks following each treatment, along with their post transcription regulation points affecting gene expression throughout development. This study suggested that glutamatergic synapse and axon guidance pathways were specifically enriched and many miRNAs and genes were significantly altered following alcohol or nicotine–alcohol perinatal exposure when compared to saline control. These results provide more detailed insight into the cell proliferation, neuronal migration, neuronal axon guidance during the infancy in rats in response to perinatal alcohol/ or nicotine–alcohol exposure.

miRNAs and their target gene profiling. Using the list of our DEGs and DEmiRs following each treatment group, predicted and validated miRNA and mRNA targets were found and plotted using MultiMiR package. After alcohol treatment, 770 miRNA-gene target pairs with 455 nodes and 1,148 miRNA-gene target pairs with 581 nodes following nicotine-alcohol treatment were found based on conserved prediction sites. Figure 2 shows the predicted networks for miRNAs and their target genes after (a) perinatal alcohol compared to the saline control group, (b) perinatal nicotine-alcohol treatment compared to the saline control group, and (c)  www.nature.com/scientificreports/ perinatal nicotine-alcohol treatment compared to the alcohol group. Among the DEmiRs, mir-30b was found to have the greatest number of connections to DEGs in the perinatal alcohol vs. saline comparison as well as the perinatal nicotine-alcohol vs. saline comparison. Figure 2d,e illustrate the miR-30b DEGs after alcohol and nicotine-alcohol exposure, respectively. miR-30b was predicted to target Gnai2 and Cotl1 following alcohol exposure and Gnai2 and Bnip3l following nicotine-alcohol exposure ( Table 2). Among the DEmiRs, mir-26b was found to have the greatest number of connections to DEGs in the perinatal nicotine-alcohol group when compared against the alcohol treatment group as shown in Fig. 2f. miR-26b was predicted to target Nxpe3 following nicotine-alcohol exposure vs. alcohol exposure ( Table 2).
Enriched pathway analysis, and target predictions. The multiMiR 51 R package was used to find the miRNA and mRNA validated and predicted target pairings. This was done using inversely regulated DEmiRs and DEGs. Using DAVID v6.8 52,53 , we further analyzed our DEG and DEmiR lists and looked across the significantly enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) 54-56 pathways following each treatment to better understand the systemic effect of perinatal alcohol or nicotine-alcohol exposure not only on the neuronal level, but also throughout other biological pathways in the body. The results summarized in Table 3, revealed many enriched KEGG pathways associated with cancer (p < 0.001 for alcohol and p < 0.01 for nicotine-alcohol), Glutamatergic synapse (p < 0.001 for alcohol and p < 0.01 for nicotine-alcohol), axon guidance (p < 0.01 for alcohol and p < 0.01 for nicotine-alcohol), and Mitogen-Activated Protein Kinase (MAPK) signaling pathway (p < 0.001 for alcohol and p < 0.001 for nicotine-alcohol). The axonal guidance KEGG pathway was significantly enriched among the upregulated genes following both treatment groups. Additionally, the glutamatergic synapse KEGG pathway was enriched among the significantly differentially expressed up and downregulated genes following both treatment groups. This information was then used to obtain further insight regarding the enriched pathways, including possible underlying biological processes, and to identify potential regulatory points across treatment groups. KEGGgraph 57 was additionally used to generate gene networks. When perinatal nicotinealcohol treatment group was compared to the alcohol treatment group as shown in Supplementary Table S1, pathways in cancer, long term potentiation, Huntington's disease and Parkinson's disease were among the many biological pathways that were significantly altered. Figure 3a-c show the miRNA-gene network from the Glutamatergic synapse after perinatal alcohol vs. saline, perinatal nicotine-alcohol vs. saline, and perinatal nicotinealcohol vs. alcohol treatments, respectively (p < 0.001 for alcohol, p < 0.01 for nicotine-alcohol, and p > 0.05 for nicotine-alcohol vs. alcohol). Figure 3d-f illustrate the axon guidance KEGG pathways after perinatal alcohol vs. saline, perinatal nicotine-alcohol vs. saline, and perinatal nicotine-alcohol vs. alcohol treatments respectively (p < 0.01 for alcohol, p < 0.01 for nicotine-alcohol, and p > 0.05 for nicotine-alcohol vs. alcohol). Additionally, we used Metascape 58 enrichment network visualization to show the intra-cluster and intercluster similarities of enriched terms and to predict the interactions among biological pathways using the significant DEGs following each treatment 58 . Metascape uses human targets as a default and analyzes different model organisms, giving us a most comprehensive human-centric database 58 . Figure 4 shows our DEG lists analyzed as human species, allowing us to predict affected downstream pathways and protein complexes within the human genome. This analysis confirmed the axon guidance pathway was affected in each treatment group. After perinatal alcohol treatment within our upregulated DEG list, axon guidance pathway was interacting with signaling by interleukins and cell projection morphogenesis pathways (Fig. 4a). Following perinatal nicotine-alcohol, our upregulated DEG lists showed axon guidance pathway showed interaction with protein localization to membrane (Fig. 4b). After perinatal alcohol treatment, our downregulated DEG lists showed that axon guidance pathway was interacting with VEGFR2 mediated vascular permeability (Fig. 4c). After perinatal nicotine-alcohol treatment when compared to the alcohol treatment group, our upregulated DEG lists showed that axon guidance pathway interacts with membrane trafficking, establishment of protein localization to organelle, and cellular responses to stress (Fig. 4d).

Discussion
Recent statistics show that 7.2 and 11.5% of mothers use nicotine and alcohol during pregnancy, respectively 15,16 . The effects of drug use are lifelong and can be severe including birth complications and developmental disabilities 14 . In this study, we focused on identifying large scale miRNA and gene expression profiling in the DA neurons of the VTA following perinatal alcohol or nicotine-alcohol exposure in male rat pups. The time frame of exposure in rats was equivalent to the three trimesters of human pregnancy. In our recent study, we have found the expression of many miRNAs and mRNAs to be significantly altered following perinatal nicotine Table 1. Top 20 significantly differentially expressed miRNAs. (DEmiRs) in DA neurons of the VTA following perinatal (a) alcohol exposure compared to saline control, (b) nicotine-alcohol compared to saline control, and (c) nicotine-alcohol compared to alcohol exposure. Benjamini-Hochberg method was used for the statistical analysis (q value < 0.05, absolute log2 fold change > 0.5).   www.nature.com/scientificreports/ study encouraged us to investigate the influence of maternal alcohol intake and exposure and expand this study to include maternal alcohol and nicotine combined exposure. Subsequently, we used a similar protocol as in our previous study 31 which investigated perinatal nicotine exposure alone compared to saline. For the alcohol group, pregnant mothers were fed a Lieber-DeCarli ethanol diet. This widely practiced and established method provides high protein ethanol diet to the pregnant mother without inducing stress or compromising the mother's health 50,59 . Regarding the nicotine-alcohol group, we used our established nicotine animal model while feeding them the liquid ethanol diet throughout the 4 weeks of gestational exposure 24,30,35 . Following our differential expression analysis of the alcohol-treated group, we found 1,257 DEGs to be upregulated, and 330 DEGs downregulated. Our nicotine-alcohol treatment differential expression analysis identified 1,771 DEGs upregulated, and 269 DEGs downregulated. Among the microRNAs, there were 51 DEmiRs that were upregulated and 39 DEmiRs that were downregulated following perinatal alcohol exposure. Following perinatal nicotine-alcohol exposure, 51 DEmiRs were upregulated and 41 DEmiRs were downregulated. Validated and predicted correlation between our DEmiRs and DEGs target pairs was performed and analyzed following each treatment group to identify possible miRNA to gene pairs and putative miRNA targets resulting from perinatal alcohol or nicotine-alcohol exposure. To better interpret the function of these DEGs within biological processes during development, pathway enrichment analysis using DAVID was performed. Finally, we analyzed the DEG lists following perinatal alcohol and nicotine-alcohol exposure groups through Metascape to further study the interactions between enriched pathways, and apply our findings to genes within the human genome. After analyzing the predicted and validated miRNA-gene target pairs, we found the miRNA miR-30b had the greatest number of connections with the target genes in both alcohol and nicotine-alcohol treatment. This suggests the potential of miR-30b in controlling the expression of many genes involved in different biological processes during neurodevelopment following nicotine and/or alcohol exposure. In our study, following the alcohol treatment miR-30b was predicted to target the DEG GNAI2, a G protein subunit alpha i2. Recently, using whole-exome sequencing analysis, a de novo heterozygous missense mutation in the GNAI2 gene was found in an individual with periventricular nodular heterotopia and intellectual disability 60 . Additionally, in vivo studies conducted with GNAI2-knockdown mice determined a lack of social interaction, recognition and increased anxiety in these mice 60 . Together, these studies suggest an important role for GNAI2 in healthy brain development. Moreover, GNAI2 has been linked to long-term depression in the neurodevelopment 61 . Another predicted DEG that was targeted by miR-30b was COTL1, which is a coactosin-like F-actin binding protein 1. COTL1 competes with cofilin to bind to F-actin, and involves neuronal migration 62 . Early-and late-born cortical neurons display distinct migratory behaviors 63 . A study conducted by Li et al. revealed that COTL1 overexpression impaired migration of both early-and late-born neurons during mouse corticogenesis, which suggests COTL1 Figure 2. Integrated analysis of the predicted and validated miRNA-mRNA target network. Following perinatal (a) alcohol exposure compared to saline control, (b) nicotine-alcohol exposure compared to saline control, and (c) nicotine-alcohol exposure compared to alcohol exposure. For both (d), following alcohol exposure compared to saline control and (e), following nicotine-alcohol exposure compared to saline control, using negative correlation, rno-miR-30b-5p was predicted to target the greatest number of genes within our DEGs. For (f) following nicotine-alcohol exposure compared to alcohol exposure, using negative correlation, rno-miR-26b-5p was predicted to target the greatest number of genes within our DEGs.  −02   PHLPP1, CSF3, FGFR2, FGFR1,  YWHAZ, FGFR3, GRB2, PPP2R5C,  FOXO3, ITGB1, CDC37, SOS2,  GYS1, PRKAA1, MLST8, FGF1,  CHUK, IL7, PKN2, COL5A3,  YWHAE, CDK2, DDIT4, MAPK1,  YWHAH, CDKN1B, GNB2, GNB1,  VEGFA, YWHAQ, PDGFRA,   www.nature.com/scientificreports/ expression of BNIP3L promoted its localization to the mitochondria, triggered a loss of membrane potential, and increased reactive oxygen species production, which ultimately leads to cell death 68 . Therefore, the BCL2 family has been an important focus of neuroscientific interest due to its potential influence on neurodegenerative pathologies. Overall, these predicted targets suggest alterations during the neurodevelopmental processes at the cellular level following both alcohol and/or nicotine exposure through modulation of genes associated with neural migration, neurodevelopment and apoptosis. After analyzing the predicted and validated miRNA-gene target pairs in the nicotine-alcohol treatment group compared to the alcohol only group, we found the miRNA miR-26b had the greatest number of connections with the target genes. miR-26b has been previously shown to play a role in the basic mechanisms of brain neuroplasticity, stress response and in the pathogenetic mechanisms of several neuropsychiatric diseases 69 . Within the nicotine-alcohol vs. alcohol treatment group, miR-26b was predicted to target NXPE3, which is a neurexophilin family of neuropeptide-like glycoproteins promoting adhesion between dendrites and axons 70 . Defects in this may result in brain abnormalities as NXPE3 has been linked to epilepsy 71,72 . This suggests alterations in neural migration and connections during the neurodevelopmental processes following nicotine-alcohol perinatal exposure compared to the alcohol perinatal exposure group at the cellular level.
Following our differential expression analysis and miRNA-gene target pairs, we further analyzed the enriched KEGG pathways to understand putative processes regulated by perinatal alcohol and nicotine-alcohol exposures. We found that glutamatergic synapse and axonal guidance pathways following perinatal alcohol and perinatal nicotine-alcohol exposure in the DA neurons of the VTA were enriched. Glutamate is the major excitatory neurotransmitter in the mammalian brain accounting for approximately 70% of synaptic transmission within the central nervous system [73][74][75] . Glutamate pathways are linked to many other neurotransmitter pathways as glutamate receptors are found throughout the brain and spinal cord in neurons and glia 76 . A subset of DA neurons in the VTA co-release DA and glutamate to the NAc and are believed to play a role in behavioral activation following stimulants, illustrating a role in drug addiction [77][78][79] . This could partly explain the enriched glutamatergic synapse pathway as alcohol directly activates the DA, playing an essential role in neurodevelopment 80 . Glutamatergic Enriched KEGG pathways were shown. Glutamatergic synapse pathway was shown in (a-c) following perinatal alcohol compared to saline control, nicotine-alcohol compared to saline control, and nicotine-alcohol compared to alcohol exposure in DA neurons, respectively, (p < 0.001 for alcohol compared to saline control, p < 0.01 for nicotine-alcohol compared to saline control, and p > 0.05 for nicotine-alcohol compared to alcohol exposure). Axon guidance pathway was shown in (d-f) following perinatal alcohol compared to saline control, nicotine-alcohol compared to saline control, and nicotine-alcohol compared to alcohol exposure in DA neurons, (p < 0.01 for alcohol compared to saline, p < 0.01 for nicotine-alcohol compared to saline, and p > 0.05 for nicotine-alcohol compared to alcohol exposure).

Scientific RepoRtS
| (2020) 10:15016 | https://doi.org/10.1038/s41598-020-71875-1 www.nature.com/scientificreports/ synapse pathway was enriched following both treatments perinatally compared to control (p < 0.001 for alcohol, p < 0.01 for nicotine-alcohol). Among the DEmiRs within the glutamatergic synapse pathway following perinatal alcohol exposure, miR-410-3p (p value of < 0.05) and miR-298-5p (p value of < 0.01) had further connections. miR-410-3p was predicted to target and inhibit TRPC1, which encodes a membrane protein forming a non-selective channel permeable to calcium and other cations 81 . TRPC1 has been shown to affect the group I metabotropic glutamate receptors pathway and auditory signal processing at the midbrain level 82 . Additionally, Xiong et al. showed that low miR-410-3p expression was associated with the chemotherapy drug, gemcitabine resistance in human pancreatic cancer xenograft tumor tissues and pancreatic ductal adenocarcinoma (PDAC) cells as well as poor prognosis in PDAC patients. TRPC1 as one of the potential targets of miR-410-3p was also significantly affected by the miR-410-3p expression modifications 83 . This data may suggest that the glutamatergic system plays a role in non-neuronal tissues including tumor biology 84 . miR-298-5p was predicted to target and inhibit HOMER3 which encodes a postsynaptic density scaffolding protein 81 , which in part regulates signal transduction and maintains extracellular glutamate levels in corticolimbic brain regions 85 . A review study by Szumlinski et al. suggests members of the Homer protein family regulates drug-induced neuroplasticity through glutamate receptor trafficking 85 . The glutamatergic synapse pathway following perinatal nicotine-alcohol exposure shows a more complex pathway with several intracellular regulatory points. Among the many DEmiRs, miR-449a (p value of < 0.05) was found to significantly inhibit genes within different clusters. miR-449a showed a strong specificity for lung, testis, and adenocarcinoma tissues and to be involved in the development of carcinoma by being a potential inducer of cell death, cell-cycle arrest, and/or cell differentiation 86 . Our data showed that miR-449a was predicted to target PRKCB, which further connects and activates MAPK3. MAPK3 is known to be involved in the control of cell proliferation, cell differentiation, and cell survival 87 . Moreover, pathways and functional linkages in the large set of genes were associated with autism spectrum disorders (ASD). Based on the common ASD genes in the MAPK (MAPK3) and calcium signaling pathways (PRKCB), the overlapping function of these two pathways in ASDs were narrowed down to voltage-gated calcium channels and calcium activated PKC 88 . DEmiRs miR-290 (p value of < 0.05) and miR-7b (p value of < 0.05) also targeted the same gene cluster as miR-449a, but were specifically connected to the genes MAPK1 and Phospholipase D1(PLD1), respectively. PLD1 has been shown to negatively affect glutamate function under oxidative stress conditions, which highlighted the role of PLDs in glutamate transporter regulation in the synaptic endings exposed to oxidative injury 89 . miR-30b-5p was predicted www.nature.com/scientificreports/ to have the greatest number of connections downstream and was also differentially expressed in this pathway, affecting different gene clusters, predicted to connect and inhibit Shank3, which plays a role in the function of synapses, ensuring signals are received by the postsynaptic neuron in the brain 90 . The glutamatergic synapse pathway was looked at following perinatal nicotine-alcohol exposure compared to perinatal alcohol exposure group (p > 0.05) for comparison. During the development of the nervous system, neurons extend their axons to reach their targets, forming functional circuits. These circuits are the basis of neural function and their faulty assembly can result in disorders of the nervous system 91 . A study conducted by Lindsley et al. suggested that ethanol disrupted the way axons respond to guidance cues effecting axon growth and elongation 92 . Our results supported this data as we found the axonal guidance pathway to be enriched following both treatments perinatally compared to control (p < 0.01 for alcohol, p < 0.01 for nicotine-alcohol). The axonal guidance pathway following perinatal alcohol exposure had many DEmiRs, among which mir-15b-5p had a significant predicted interaction with a p value of < 0.001 targeting DEG GNAI1. Additionally, Lewohl et al. reported miR-15b to be up-regulated in the prefrontal cortex of human alcoholics 93 . This data may suggest possible similarities in the altered expression of genes involved in the neural networks of an addicted brain and a developing brain exposed to alcohol. miR-495 (p ≤ 0.02) was predicted to target and inhibit two gene from different gene clusters, KRAS, part of the RAS oncogene family, involved in cell growth, maturation, and death 94 , and Wnt4, which plays critical roles in many biological processes including embryonic development 95 . KRAS was further predicted to activate MAPK1/3, playing a role in regulating multiple physiological processes including mitosis, cell differentiation, and cell survival. The axon guidance pathway following perinatal nicotine-alcohol exposure compared to perinatal alcohol exposure group was not significant (p > 0.05).
Following perinatal nicotine-alcohol exposure, the axonal guidance pathway showed the significantly differentially downregulated mir-466b-5p (p value of ≤ 0.0001) inhibiting Integrin Subunit Beta 1 (ITGB1), which in turn activated UNC5B and PTK2. UNC5B encodes a protein that is part of the dependence receptors (DPRs) proteins, which are said to be involved in embryogenesis and cancer progression [96][97][98][99] . PTK2 gene activation is said to be a crucial step early on in the cell growth and intracellular signal transduction pathways 100,101 . miR-34c-5p (p < 0.05) was predicted to target and inhibit SLIT1, which was differentially upregulated and predicted to connect and activate SLIT2, SRC-like kinase FYN and ROBO1 genes. FYN has been implicated in the control of cell growth 102 and has been linked to cancer pathogenesis 103 . The Slit family of secreted glycoproteins were originally identified in the nervous system functioning as axon guidance cues and branching factors during development regulating neuronal axon guidance, neuronal migration, cell proliferation and cell motility through its binding to Robo receptors 104,105 . The overall results suggest neuronal development to be highly modulated at many putative points resulting in alterations within many biological processes crucial to development such as cellular growth, differentiation, signal transduction, synapses, and cell survival.
Pathways in cancer, Wnt signaling pathway, long-term potentiation, Huntington's disease, and Parkinson's disease were among the many significantly enriched KEGG pathways following the nicotine-alcohol perinatal exposure group compared to the alcohol perinatal exposure group. Prenatal alcohol exposure has been shown in previous studies to disrupt Wnt signaling pathway and has been a determinant of FASD 106,107 . In the central nervous system, Wnt signaling is known to modulate neuronal proliferation, migration, adhesion, differentiation, and axon outgrowth [108][109][110][111][112] . FASD have also been linked to abnormal neuronal plasticity responsible for normal wiring of the brain and involved in learning and memory 113 . Many studies have shown this by corresponding disruptions in long term potentiation with acute or chronic perinatal alcohol exposure 114-117 . We further looked at the regulation and connection of the enriched axon guidance pathway with other biological pathways and protein complexes by using Metascape, which converts the given gene lists from rat into human Entrez gene IDs. Following this analysis, axon guidance from our upregulated DEGs following perinatal alcohol exposure was shown to directly connect to cell projection morphogenesis. Cell morphogenesis is the process in which the neurons are generated, organized and targeted to a specific site in response to attractive or repulsive cues 118 . Axonogenesis refers to the morphogenesis of shape or form of the developing axon, which carries efferent action potentials from the cell body towards target cells. Axon guidance is one of the important parts of the axonogenesis during which the migration of an axon is directed to a specific target. This data also confirmed that an upregulation in these two processes closely linked them to each other. Axon guidance from our downregulated DEGs following perinatal alcohol exposure was connected to VEGFR2 mediated vascular permeability. Alcohol-induced premature retraction of the radial glia in the deep cortex and alcohol-induced retardation in neuronal migration have been already shown in the literature [119][120][121][122] . Considering VEGF is a chemoattractant for commissural axons in vitro and in vivo 123 and that growing axons are guided to their targets by attractive and repulsive cues, it is very likely that if there is a downregulation in the DEGs due to alcohol exposure, these two pathways will influence each other. Axon guidance from our upregulated DEGs following perinatal nicotine-alcohol exposure was directly connected to protein localization to membrane, which could suggest a broad modulation among many pathways affecting protein function and translation. Axon guidance was also connected to membrane trafficking, establishment of protein localization to organelle, and cellular responses to stress.
In conclusion, we have conducted a large-scale miRNome and transciptome study following perinatal alcohol and nicotine-alcohol exposure in the DA neurons of the VTA of male rat pups. We have investigated transcriptional and post transcriptional alterations and putative regulatory points within neurodevelopmental pathways in the postnatal brain and possible disruptions within biological pathways systemically. We identified enriched biological pathways following each treatment, and downstream gene network interactions between these significant pathways within the human genome. Our study suggested Glutamatergic synapse and axon guidance pathways to be significantly enriched and many miRNAs and genes were altered following nicotine or alcoholnicotine exposure perinatally. Cell growth, proliferation, neuronal migration, neuronal axon guidance, and cell www.nature.com/scientificreports/ survival cues were among pathways in which many genes and miRNAs were significantly altered in response to perinatal alcohol or nicotine-alcohol exposure. Our nicotine-alcohol exposure compared to saline group showed the nicotine addiction pathway was enriched, which was also seen in our previous nicotine only study, comparing perinatal nicotine exposure to saline 31 . Additionally, the glutamatergic synapse pathway was enriched in all groups (nicotine, alcohol, and combined nicotine-alcohol) when compared to saline 31 . Although, our results using rat model of 3-trimester gestational exposure to both alcohol and nicotine indicate that perinatal alcohol and nicotine exposure alters the expression of miRNAs and genes in infant rats, it should be noted that it is important to validate these models by translating them into human studies. Variable factors in this study include the alcohol intake, nicotine doses and body weight variations across dams which could potentially affect the results. Additionally, we limited our initial study to only include male pups. Such a limitation requires further studies to be conducted using female pups to explore gender differences following perinatal alcohol and/ or nicotine exposure. Further studies need to be conducted to better understand the systemic putative pathway regulation points at a molecular level from gene expression profiling to protein translation, as well as to investigate therapeutic approaches that target disorders associated with gestational addictive substance exposure such as FASDs.

Materials and methods
Animal treatment. All experimental protocols and surgical procedures were approved by the Institutional Animal Care and Use Committee (IACUC) and the University of Houston Animal Care Operations (ACO) and were performed in accordance with accepted guidelines and regulations. Sprague Dawley (SD) rats are one of the most commonly used outbred rat strains for biomedical research, providing 95% or greater accuracy on timed-pregnant gestation and ideal for safety and efficacy testing, aging, behavior, reproduction and surgical modifications 124,125 . Pregnant female SD rats purchased from Charles River (Charles River, Wilmington, MA, USA) were housed in the animal facility and were maintained at 22 ± 2 °C with 65% humidity, on a 12-h light/12-h dark cycle. The animal treatment method has been further detailed in Keller et al. 31,32 . Briefly, rats were acclimated to the animal facility upon arrival for 72 h before the subcutaneous insertion of an osmotic pump (Alzet, Cupertino, CA, USA) containing either nicotine hydrogen tartrate (Sigma-Aldrich, St. Louis, MO, USA), which released nicotine at a rate of 6 mg/kg/day in order to simulate the nicotine plasma level found in moderate smokers 31,32 , or an equal volume of saline for control 126 . Liquid diet containing 36% kcal from ethanol F1265SP or F1264SP control, purchased from Lieber-DeCarli (Bio-Serv, Flemington, NJ, USA) was placed in Richter feeding tubes (Bio-Serv), and was gradually introduced to the pregnant mothers based on the protocol provided from Bio-Serv. This liquid diet model has reliably produced blood alcohol concentrations (BACs) between 80 and 180 mg/dl in rats, which are accompanied by neurological deficits similar to what is observed in children with FASD 50,[127][128][129][130][131][132][133] . Pregnant dams on average consumed around 80-100 ml of liquid diet per day. A total of 12 dams were used in this study. 4 were used for saline, 4 for alcohol, and 4 for nicotine-alcohol. Pups would be exposed to either alcohol or nicotine and alcohol for four weeks starting from gestational day 6 to postnatal day 14, equivalent to the three trimesters in human gestation during which rapid brain growth and synaptogenesis occur 38,49 . The central nervous system development of a rat at postnatal age of 7-14 days is suggested to correspond approximately to the human brain at term 134 . The male pups collected from P10 to P14 were all considered as "infants", therefore, they were pooled from each litter for this study 135,136 . Supplementary Figure S1 summarizes the overall methodology. Perinatal drug exposure including nicotine, alcohol and combined exposure has been known to produce changes in a sex dependent manner [137][138][139][140] . Previous studies considering early-life exposure to these drugs have been done predominantly on male subjects or rodents 141 . Perinatal nicotine exposure has been shown to have more significant deleterious effects on the cholinergic and serotonergic markers in males than females 140,142 . Drug reward is also shown to be altered by prenatal nicotine exposure with increased preference for nicotine in males than females 140,143 . For this reason, males only were examined in this study. Male and female rat pups were distinguished by a larger genital papilla and longer ano-genital distance in male vs. female pups 137 . Postnatal 10 to postnatal 14-day old male pups were anesthetized using isoflurane gas before decapitation on a VT1200 semiautomatic vibrating blade microtome (Leica, Nussloch, Eisfeld, Germany). 1 mm thick horizontal brain slices containing the VTA were sliced and 1 mm biopsy punch (Integra Miltex, VWR, Radnor, PA, USA) was used to collect the VTA bilaterally. Brain punches from 4 to 7 pups from each litter were pooled and placed on ice in Hibernate A (Gibco, Thermo Fisher Scientific, USA) to preserve and maintain cell viability. Brain punches were pooled for each litter and a total of four samples, 16-28 pups/sample for each alcohol, nicotine-alcohol, and saline treated groups were collected, processed, and analyzed, for a total of 48-84 pups used in this study.
Brain slice preparations and FACS cell sorting. Brain tissue punches containing the VTA were collected and dissociated into a single cell suspension before sorting using FACS as previously reported by Guez-Barber et al. 144 . Briefly, collected tissue punches were dissociated in Accutase (Gibco, Thermo Fisher Scientific, Waltham, MA, USA) and shaken at 4 ℃ for 30 min. Cells were centrifuged and pelleted at 425×g and resuspended in Hiberante A medium (Gibco). Cell aggregates were then further dissociated through gentle pipetting with increasingly smaller pipette tips; the supernatant containing single cells were collected. Cellular debris was then removed using serial filtration. First, the cell suspension was run through a pre-wetted 100 µm cell strainer and then through a pre-wetted 40 µm cell strainer while on ice. Further removal of small cellular debris was done by density centrifugation. The cell suspension was added to the top of a three-density gradient which was made using Percoll (GEHealthcare, VWR, USA) and centrifuged for 3 min at 430×g. After centrifugation, the cloudy top layer, which contained cellular debris was removed. The remaining cell suspension was pelleted through centrifugation at 550×g for 5 min. www.nature.com/scientificreports/ To fix cells for immunolabeling, cells were resuspended in 1 ml of Hibernate A and 1 ml of cold absolute ethanol, gently vortexed and stored on ice for 15 min. Cells were incubated and labelled 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) and rotated for 30 min at 4 ℃. Cells were then washed with PBS and centrifuged at 950×g for 3 min before they were resuspended in PBS. Flow cytometry was performed on an (LSR II) FACS Aria (BD Biosciences) instrument and analyzed using FlowJo software at the Baylor College of Medicine Cytometry and Cell Sorting Core (One Baylor Plaza, Houston, TX, USA). DA neurons with intact nuclei were labelled and sorted based on their positive double staining of NeuN and TH, populations were distinguished by their forward and side scatter, and two-parameter density plots were measured with gating parameters set at around 10 3 for NeuN FITC and 10 4 for TH-PE expression.
RNA extraction. Following FACS, cells were pelleted by centrifugation at 2,650×g for 8 min at 18 ℃ and total RNA was extracted using miRNeasy Micro Kit (Qiagen, Hilden, Germany) including DNAse treatment following manufacturer's instructions. A NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific) was used to check the RNA purity and quantity according to the optical density (OD) of each sample at 260 nm and 280 nm. Only samples with a 260/280 ratio of 1.9 or greater were used in experiments.
Microarray preparation, labeling, and hybridization for gene and microRNA expression profiling. All kits that were used for both mRNA and miRNA gene expression analysis were purchased from Agilent (Santa Clara, CA, USA) unless stated otherwise. SurePrint G3 Rat Gene Expression v2 8 × 60 K microarray (ID: 074036) with 30,584 unique genes was used for the mRNA expression profiling using 25 ng of total RNA. Samples were prepared and labeled according to manufacturer's instructions for the One-Color Microarray-Based Gene Expression Analysis using the One-Color Low Input Quick Amp Labeling kit with RNA Spike-Ins. The labeled amplified complimentary RNA was then purified using the RNeasy Mini Kit (Qiagen) and quantified using a NanoDrop 2000 spectrophotometer. The cRNA yield and specificity was then calculated according to manufacturer's instructions. Microarrays were then hybridized for 17 h at 65 ℃ using the Gene Expression Hybridization kit according to manufacturer's instructions. miRNA expression profiling was performed using an 8 × 15 K Rat miRNA Microarray, Release 21.0 (ID: 070154) containing 758 mature miRNAs. A starting quantity of 100 ng total RNA containing miRNAs was used following manufacturer's instructions of the MicroRNA Microarray System with miRNA Complete Labeling and Hyb Kit. Purification of the labeled RNA was done on a Micro Bio-Spin P-6 gel column (Bio-Rad, Hercules, CA, USA) following manufacturer's instructions. Samples were then dried using a vacuum concentrator with heater at 50 ℃ and hybridized at 55 ℃ for 20 h. Slides were then washed following manufacturer's instructions using Gene Expression Wash Buffers containing Triton X-102. Both gene and miRNA expression slides were scanned following hybridization using a G4900DA SureScan Microarray Scanner. The raw microarray dataset was then collected from the resulting images using the Feature Extraction Software v12.0.1.

Data analysis.
Analysis and comprehension of the genomic data was performed using Bioconductor packages in R version 3.6.1 145 . The limma 146 package was used for mRNA expression data loading and following procedures including the removal of outliers, background correction and quantile normalization were performed. We defined low expression as any intensity less than 75% brighter than 90% intensity of negative controls. Replicates were removed and a total of 11,791 genes were collected for analysis. Fold change and standard errors were calculated using lmFit function, which fits multiple linear models by weighted least squares. Standard errors were moderated using eBayes function, which computes log-odds of differential expression by an empirical Bayes model. DEGs were identified using a series of p values (i.e. 0.05, 0.01, 0.001), with the minimum log2 fold change > 1 and adjusted using the BH method. DEmiRS were identified using a series of p values (i.e. 0.05, 0.01, 0.001), with the minimum log2 fold change > 0.5 and adjusted using the BH method.
The AgiMicroRna 147 package was used for the miRNA expression data loading and processing. Raw data was loaded using readMicroRnaAFE function, and preprocessed using rmaMicroRna function, which implements the robust multi-array average (RMA) algorithm. Data was then filtered using filterMicroRna function and only detected genes which were expressed in at least 50% of samples, with higher intensity than the mean value of negative control + 1.5 standard deviations, were collected for analysis. After preprocessing, 332 miRNAs remained. Similar to mRNA data processing, linear model was fitted to the miRNA expression data and moderated statistics were calculated using eBayes. Differential expression was identified using a series of p values (i.e. 0.05, 0.01, 0.001) and adjusted using the BH method.
Comprehensive analysis of differentially expressed gene and microRNAs. MultiMiR 51 package was used to identify the predicted and validated miRNA-mRNA pairs based on the inversely correlated regulation between miRNA and target genes. miRNA-mRNA pairs were identified using a series of p values (i.e. 0.05, 0.01, 0.001) and labeled accordingly. Network visualization was conducted using Gephi 148 , which is an open-source software for network visualization and analysis. Three colors were assigned to edges with different p values (i.e. green for p = 0.05, blue for p = 0.01, and red for p = 0.001). The color blue was then assigned to all the miRNAs and yellow for all the mRNAs. The size of the vertex was determined based on an out-degree, which represents the number of edges formed by an incident. A topology filter was then applied to filter out the vertex, which do not have any connections. Functional enrichment analysis was performed on our integrated network of DEGs and DEmiRS using DAVID version 6.8 52,53 . Enriched KEGG pathways were then identified and further analyzed using KEGGgraph 57 . Metascape 58 was used on our DEG list following alcohol and nicotine-alco-Scientific RepoRtS | (2020) 10:15016 | https://doi.org/10.1038/s41598-020-71875-1 www.nature.com/scientificreports/ hol treatment to show enriched pathways as a network, therefore further understand relationships among the enriched pathways and their correlation to each other through their downstream connections.