Introduction

The ventral tegmental area (VTA) contains dopamine cells (VTA^DA) that play key roles in several behavioral processes, including reward, aversion, and motor control^1–3. VTA^DA cells are heterogeneous, with recent studies identifying subpopulations of VTA^DA cells that differentially contribute to adaptive and pathological behaviors^4–7. Numerous studies have mapped the input-output architecture of VTA cells generally, as well as subpopulations of constituent DAergic, GABAergic, and glutamatergic cells^7–12. These studies have helped us to understand the biased input and discrete output logic of the VTA, as well as generate hypotheses regarding the function of circuits including the VTA.

In addition to normal adaptive behaviors, VTA^DA cells critically contribute to a variety of behaviors related to substance use disorder (SUD), including reward, sensitization, and withdrawal elicited by drugs of abuse^13–18. For a variety of drugs of abuse, even a single exposure triggers long-lasting changes in synaptic plasticity onto VTA^DA cells^19,20. These changes include both a long-lasting enhancement of excitatory synaptic strength onto VTA^DA cells²⁰, and development of inhibitory synaptic plasticity^21,22. These changes are essential for a variety of drug-induced behaviors, linking plasticity at VTA^DA cells with downstream behavioral changes^22–24. Furthermore, the changes in excitatory synaptic plasticity specifically occur onto some cell types and not others⁵, indicating that these adaptations are likely cell-type and pathway-specific. However, it remains unknown which sets of inputs are modified, and if these changes are consistent across abused drugs. Furthermore, synaptic plasticity is not the only type of long-lasting change that can be induced by drugs of abuse. Unfortunately, even less is known about how other forms of experience-dependent changes, including elevations in input activity, may contribute to addiction-related behaviors.

We recently showed that a single dose of cocaine causes long-lasting changes in cellular activity in inputs to VTA^DA cells^15,16. These changes occur from the globus pallidus external segment (GPe) and bed nucleus of the stria terminalis (BNST) projections to the VTA that contribute to drug-induced reward and sensitization, or withdrawal anxiety, respectively^15,25. Notably, while these studies focused on inputs from the GPe and BNST, the viral mapping results also demonstrated the existence of other input changes to the VTA that were not pursued in those studies. These unexplored connections may provide critical information about how drugs alter connectivity onto VTA^DA cells, which in turn contributes to drug-induced behavioral changes^15,16.

In this study, we explore the brain-wide changes in rabies-labeled inputs to VTA cells elicited by a variety of different drugs, including drugs of abuse, an anesthetic mixture of ketamine/xylazine (K/X), and fluoxetine, a psychoactive but non-abused control substance. Our study provides a first-of-its-kind exploration into how even a single drug exposure can change brain-wide connectivity patterns, in addition to highlighting some concerns about the potential effects of K/X anesthesia on motivated behavior.

Materials and methods

Experimental Procedures

Animals

Mice were housed on a 12-hour light–dark cycle with food and water ad libitum. Males and females from a C57/BL6 background were used for all experiments in approximately equal proportions. Mice were approximately 3-4 months of age at the time of experiments. All surgeries were done under isoflurane or ketamine/xylazine anesthesia. All procedures complied with the animal care standards set forth by the National Institute of Health and were approved by the University of California, Irvine’s Institutional Animal Care and Use Committee (IACUC).

Stereotaxic surgery

Mice were anesthetized with 3-4% isoflurane and maintained during surgery at 1-1.5% isoflurane, or with a mixture of ketamine (100 mg/mL) and xylazine (5 mg/mL). Mice were secured in a stereotaxic apparatus (Stoelting). Under aseptic conditions, guide holes were drilled, and viruses were infused into the target sites using a glass capillary attached to a microinjection pump (WPI, UMP3T). 100 nL of AAVs or 500 nL of EnvA-pseudotyped RABV were infused at a rate of 1.6 nL/s. The glass capillary remained in place for 5 min following the infusion to allow for virus diffusion. After infusion, the surgical incision sites were closed with either sutures or Vetbond tissue adhesive (Patterson Veterinary). Vetameg (flunixin, Patterson veterinary) was administered for pain management and topical bacitracin was applied to prevent infection at the incision site.

For viral injections mapping experiments from a defined cell type, 100 nL of a 1:1 volume mixture of AAV₅-CAG-FLEx^loxP-TCB and AAV₈-CAG-FLEx^loxP-RABV-G was injected into the VTA of DAT-Cre²⁶, GAD2-Cre²⁷, or vGluT2-Cre mice²⁸. Thirteen days later, animals were given a single injection of a drug of abuse, fluoxetine, or saline. The following day, 500 nL of EnvA-pseudotyped RABV was injected into the VTA. Animals were sacrificed five days later. For cTRIO experiments shown in Figure 4, during the first surgery AAV₅-CAG-FLEx^FRT-TCB and AAV₈-CAG-FLEx^FRT-RABV-G were injected into the VTA, and CAV2-FLEx^loxP-Flp was injected into an output site (NAcLat, amygdala, NAcMed, mPFC) during the same surgery. All other procedures were the same as above. RABV-labeled cells throughout the brain were manually quantified, as done previously^7,10,15.

The injection coordinates and volumes used were as follows:

VTA: 100 nL AAV and 500 nL RABV, AP –3.2, LM 0.4, DV –4.2

NAcMed: 250 nL, AP +1.78, L0.4, DV –4.1

NAcLat: 250 nL, AP +1.45, LM 1.75, DV –4.0

Amygdala: 500 nL, AP –1.43, LM 2.5, DV –4.5

mPFC: 1 μL into mPFC (two injections of 500 nL, one at AP +2.15, LM 0.27, DV –2.1 and another at AP +2.15, L0.27, DV –1.6

The titers of the viruses used, based on quantitative PCR analysis, were as follows: AAV₅-CAG-FLEx^loxP-TC, 2.4 x 10¹³ genome copies (gc)/mL;

AAV₈-CAG-FLEx^loxP-G, 1.0 x 10¹² gc/mL;

AAV₅-CAG-FLEx^FRT-TC, 2.6 x 10¹² gc/mL;

AAV₈-CAG-FLEx^FRT-G, 1.3 x 10¹² gc/mL;

CAV-FLEx^loxP-Flp, 5.0 x 10¹² gc/mL.

The titer of RABVΛ1G was estimated to be 5.0 x 10⁸ colony forming units (cfu)/mL based on serial dilutions of the virus stock followed by infection of the 293-TVA800 cell line.

Drug administration

Cocaine was administered (intraperitoneal injections) at a dose of 15 mg/kg, morphine at 10 mg/kg, methamphetamine at 2 mg/kg, amphetamine at 10 mg/kg, nicotine at 0.5 mg/kg, and fluoxetine at 10 mg/kg.

Quantification of RABV data

RABV mapping data were analyzed as described previously^7,10. All RABV mapping data included were previously published^7,10,15, except for input tracing from mice treated with a single dose of methamphetamine. A total of 118 samples were included in this study (Table 1). Exclusion criteria included mice in which fewer than 200 input cells were labeled, or in cases where injections were mistargeted (e.g., starter cells in the substantia nigra pars compacta). Experimenters were blinded to the experimental condition when data quantification was being performed. For starter cell center of mass calculations, we mapped the distribution of starter neurons, as defined by TVA-mCherry expression along the A-P axis, and in the sections at the center point along this axis, defined the center of labeling in along the M-L axis.

CRISPRi experiments

Plasmid construction

The CRISPRi system was adapted from Kumar’s method²⁹. pJEP317 pAAV-U6SaCas9gRNA(SapI)-EFS-GFP-KASH-pA was obtained from Addgene and was used as a template to clone pJEP317-pAAV-U6-Cacna1eAC-EFS-GFP-KASH-pA using the Q5 Site-directed mutagenesis kit (primers: 5’-acctccgtaaGTTTTAGTACTCTGGAAACAG-3’ and 5’-atgaacccttGGTGTTTCGTCCTTTCCAC-3’) to insert the gRNA sequence (AAGGGTTCATACCTCCGTAA) for Cacna1e. The gRNA sequence was selected from the IDT Predesigned Alt-R™ CRISPR-Cas9 guide RNA resource (https://www.idtdna.com/site/order/designtool/index/CRISPR_PREDESIGN). The plasmid for CRISPRi (pJEP309-pAAV-EFS-dSaCas9-KRAB-Dio-pA), pAAV2-8, and pAAV2-5 were obtained from Addgene. pHelper plasmid was a gift from Matthew Banghart at the University of California, San Diego.

AAV production, purification, and titering

Generation of AAVs from plasmids and AAV Purification by Iodixanol Gradient Ultracentrifugation were carried out as described on Addgene’s website with minor modifications. Briefly, pHelper, pAAV2-5/pAAV2-8 and pJEP309/pJEP317-pAAV-U6-Cacna1eAC-EFS-GFP-KASH-pA plasmids were transfected into HEK293T cells using PEIMAX (Polysciences, cat# 24765-100). The molar ratio for pHelper:pAAV2-5:transfer plasmid was 1:1:1, and the ratio of total plasmids to PEIMAX was 1:3. 72 hrs post-transfection, cell lysate and medium precipitated with PEG were purified by iodixanol gradient. Purified AAVs were titered using primers against viral ITRs (5’-GGAACCCCTAGTGATGGAGTT-3’ and 5’-CGGCCTCAGTGAGCGA-3’).

Viral injection

Prior to surgery, DAT-Cre mice were anesthetized with isoflurane (4% induction, 1-1.5% maintenance) and secured on a stereotaxic frame for precise localization of injection sites. For assessment of in vivo genome editing, two conditions were generated. For the experimental condition, 600 nL of a 1:1:1 volume mix of AAV₈-EFS-dSaCas9-KRAB-Dio, AAV₉-CamKII-0.4-Cre-SV40 (Addgene, cat# 05558-AAV9) and AAV₅-U6-Cacna1eAC-EFS-GFP-KASH was injected into the NAcLat. For a no gRNA control condition, 400 nL of a 1:1 volume mix of AAV₈-EFS-dSaCas9-KRAB-Dio and AAV₉-CamKII-0.4-Cre-SV40 was infused into the NAcLat.

For mapping inputs to VTA cells in mice with Cacna1e downregulation in NAcLat, during a single surgery, 300 nL of a of a 1:1 volume mix of AAV₅-FLEx^loxP-TC and AAV₅-FLEx^loxP-RABV-G was injected into the VTA bilaterally, 600 nL of a 1:1:1 volume mix of AAV₈-EFS-dSaCas9-KRAB-Dio, AAV₉-CamKII-0.4-Cre-SV40, and AAV₅-U6-Cacna1eAC-EFS-GFP-KASH was injected into the left NAcLat, and 400 nL of a 1:1 volume mix of AAV₈-EFS-dSaCas9-KRAB-Dio and AAV₉-CamKII-0.4-Cre-SV40 was injected into the right NAcLat. This study design enabled controls within the same mice to prevent systematic biases due to differences in treatment of the mice. After 14 days, RABVι1G was injected into the VTA. Mice were sacrificed five days later, and brains processed as described below.

Brain section preparation, imaging, and cell counting

Five days after RABV injection, mice underwent transcardial perfusion with 1x PBS followed by 4% formaldehyde in 1x PBS. Brains were extracted and post-fixed in 4% formaldehyde for 24 hours, and then dehydrated with 30% sucrose in 1x PBS. Brains were cut into 60 µm slices using a cryostat and mounted on Superfrost Plus^TM microscope slides (Fisher, cat # 1255015). Once the brain sections were fully dried on the slide, they were rehydrated in 1xPBS containing DAPI at a dilution of 1:1000 (ThermoFisher, cat# D1306). After 10 minutes, the liquid was removed, and Fluoromount-G (SouthernBiotech, cat# 0100-01) was applied to the sections. Coverslips (Thermo Scientific, cat# 152460) were used to cover the brain sections on the slides. Images of the entire slides were captured using an Olympus IX83 microscope, and manual cell counting was performed on the acquired images.

Determination of CRISPRi-mediated knockdown efficiency

To determine the efficacy of CRISPRi knockdown of Cacna1e expression levels, 10 days following AAV injections into the NAcLat, NAcLat brain tissue was dissected out from the brains of injected mice, and tissue was pooled (5 mice for Cacna1e gRNA and 3 mice for no gRNA control). Total RNA was extracted with TRIzol, and 50 ng of total RNA was used for RT-qPCR using the Luna Universal One-Step RT-qPCR Kit (New England Biolabs, cat# E3005L) following the manufacturer’s instructions. Phosphoglycerate Kinase 1 (PGK1) was used as the housekeeping gene for normalization of Cacna1e mRNA levels.

Dimensionality reduction of RABV input data

For RABV input mapping experiments, quantified inputs cells were binned into 22 regions: Cortex (anterior cortex), NAcMed, NAcLat, NAcCore, DStr, VP, PO, Septum, BNST, EAM, EP, GPe, PVH, LHb, MHb, CeA, LH, ZI, DR, LDT, PBN, and DCN. Altogether, this dataset consists of the input cell counts for these 22 regions across 118 mouse brains, with n=4 or 5 for each condition.

RABV input tracing data were reported as the percentage of total quantified labeled RABV inputs, as reported previously^7,10. These data were further scaled using Z-score normalization. Principal Component Analysis (PCA) was then implemented on these normalized data using the scikit learn library in Python with default parameters³⁰. Weights for PCA models were visualized as heatmaps, reflecting the contribution of each feature to each principal component. To test for differences between groups visualized in PC space, we assessed the value for each brain for each PC and if comparing 2 conditions, used an unpaired t-test, and if >2 conditions, a one-way ANOVA and if p < 0.05, followed by pairwise t-tests with multiple comparisons corrections (Dunnett or Tukey, as appropriate). Input values were the first 3 principal components, and significant differences between conditions for each principal component (drugs versus controls, for example) were denoted as having a p-value of less than 0.05. Results are reported where relevant in the text or in Table S2.

Uniform Manifold Approximation and Projection (UMAP), implemented using the UMAP library in Python, was used to generate alternative visualizations of the RABV input data³¹. To choose and standardize our UMAP parameters across analyses, we investigated which methods placed mouse brains of the same type and treatment close together most consistently. For the distance metric parameter, the metrics we found to perform best were Euclidean and cosine, but we chose to use the Euclidean metric for consistency with previous analyses³². For the number of neighbors parameter, we chose a large number, ⅓ the size of each of our datasets, for a more global view of the data. Other parameters were kept at the default values. To account for the stochasticity of individual embeddings, we applied UMAP 20 times to each data set and found the results to be very similar, as found in previous analyses³². The results were then presented as correlograms. Clustered heatmaps and plots of these data were made using a combination of the Python packages seaborn³³ and matplotlib³⁴. The seaborn clustermap function was implemented with default values, including Euclidean distance as the distance metric.

To better visualize the distribution of a given condition in PCA and UMAP space, we overlaid the points with ellipsoids. Ellipsoids were centered at the average coordinate of a condition and stretched one standard deviation along the primary and secondary axes.

Brain Region Selection in the Allen Mouse Brain Atlas

To compare our data with connectivity data from the Allen Institute, we converted our 22 brain regions into regions or groups of regions defined by the Allen Atlas. This was necessary as our region annotations defined while doing cell counts of the RABV data were at times different from the Allen Mouse Brain Atlas-defined regions. 22 brain regions were used for RABV analyses, but these were reduced to 20 brain regions (NAcMed, NAcLat, and NAcCore were merged into the NAc) for comparison with Allen Atlas data. In several cases, multiple Allen-defined regions were used to constitute a single region that we defined, as in the case of the cortex. A summary of these conversions is as follows:

Allen Mouse Brain Connectivity Atlas projection portrait construction

Projection portraits of the VTA were made using data from the Allen Mouse Brain Connectivity Atlas³⁵, as described previously³². Briefly, data from 20 input brain regions, as defined in the above table, were downloaded and used for analyses. Cre-based and non-Cre-based experiments representative of each brain region were chosen for visualization. Labeled cell bodies at the site of injection often spread into adjacent brain regions or only covered a part of a region. Several experiments were chosen to account for this variation and give a better estimate for each region’s projections. Experiment IDs are recorded in Table S1 in our previous work³². Connectivity of brain regions of interest to the VTA was visualized by showing a representative slice of the VTA and coloring the image based on high or low connectivity according to averaged AAV tracing data. AAV tracing data were averaged across experiments and across brain regions that grouped together according to weights assigned by PCA, as shown in Figures 1, 3, and 4.

Figure 1:

Figure 2:

Figure 3:

Figure 4:

Data for projection portraits and for subsequent network construction was accessed using the python Allen SDK package, developed and maintained by the Allen Institute for Brain Science.

Gene Expression Data Access and Processing

Gene expression data was downloaded from the Allen Institute’s Anatomic Gene Expression Atlas (AGEA)³⁶, an in situ hybridization dataset which has annotated expression of 19,932 genes to every brain region defined in the Allen Mouse Brain Atlas. We used code written by Eli Cornblath (https://github.com/ejcorn/mouse_abi_tool) to access and download gene expression data for our regions of interest, as done in a previous study³⁷. To ensure only high-quality data were included in our analyses, we filtered experiments based on metrics used in previous analyses of this dataset³⁸. We retained genes for which there were data from coronal sections, or where the Pearson correlation coefficient across experiments was over 0.5 for a particular gene in our regions of interest. This filtering left us with 4,283 unique genes for the 36 Allen-defined brain regions to be aggregated to compare with our 22 quantified input sites from the RABV mapping data. We averaged data among brain regions that were defined differently in the atlas versus our RABV quantifications, which resulted in a total of 20 brain regions for comparative analyses, as detailed above. Data were normalized using a scaled sigmoid function, a method that is robust to outliers and was used in previous analyses of these data³⁹.

Gene expression data analysis and comparison to RABV data

To investigate how basal gene expression might be correlated to connectivity changes after drug exposure, we used processed data from the AGEA. For each experimental condition, the RABV data for the drug-treated animals were compared to those of the saline-treated controls. To quantify this comparison, we computed a ratio, as used in Figures 6 and S7-9, by dividing the average RABV counts in each brain region for the experimental condition by the average RABV counts in each brain region for the saline condition. This resulted in a vector of 22 ratios for each experimental condition, one ratio for each brain region. We also computed a difference, as used in Figures 5 and S4-S6, by subtracting the average RABV counts in each brain region for the experimental condition by the average RABV counts in each brain region for the saline control condition. This resulted in a vector of 22 differences for each drug comparison, one difference for each brain region.

Figure 5:

Figure 6:

We calculated the correlation of the condition-control differences with gene expression data for each region. The Spearman correlation coefficient between the gene expression across regions and our RABV differences was used because the correlative relationship between these values cannot be assumed to be linear. A similar method was used in previous analyses³⁸. A high positive correlation between a gene and an experimental condition ratio was interpreted as a greater likelihood of a brain region to exhibit RABV labeling changes after drug exposure if this gene is highly expressed before drug exposure. A negative correlation between a gene and an experimental condition ratio was taken to mean this region would be more likely to exhibit RABV labeling changes after drug exposure if the gene expression level is low.

We performed GO analysis on the top 50 most highly correlated genes and on the top 50 most negatively correlated genes. GO analysis was performed using the EnrichR package⁴⁰. Lists of genes with highly positive or negative correlations to our condition differences were compared to GO annotations using the coding language R. The 2021 version of two databases were queried: GO Biological Process and GO Molecular Function. GO enrichment analysis results were visualized in dot plots using the package ggplot2, as shown in Figures 5 and S4-S6.

We next compared the drug/isoflurane ratios for a group of addictive drugs to gene expression of subsets of genes of interest (Figure 6, S7). Similar analyses were done for ketamine and fluoxetine (Figure S8, S9). The subgroups of genes used were synapse-related genes⁴¹, ion-channel related gene⁴², neurotransmitter receptor-related genes⁴³, differentially expressed genes related to substance dependence GWAS⁴⁴, and exo- and endocytosis-related genes⁴⁵. Other subsets of genes of interest were defined using the Guide to Pharmacology database⁴⁶. Correlation coefficients (Spearman’s ρ) and p-values are reported in the figures.

Results

Dimensionality reduction methods identify differences in inputs of VTA cell populations

We and others previously showed that DAergic, GABAergic, and glutamatergic cells in the VTA receive inputs from similar brain regions, with some quantitative differences that varied depending on the study^10,12. In our work we also showed that connectivity differences in the VTA are largely a consequence of spatial organization of cells within the VTA^7,10,32. The fact that the distributions of these three VTA cell populations are largely spatially overlapping may explain why the global inputs to each population are similar. However, one limitation in interpreting these studies is that they focused on individual comparisons of input counts on a per region basis. While this is a common approach in analyzing viral tracing data, and useful for both hypothesis generation and testing, it fails to account for the higher-level organization of input patterns. Our goal here is to go beyond considering single brain regions and better understand the holistic connectivity patterns of different cell populations or drug treatment conditions.

Our first goal was to explore the potential differences in input patterns generated using RABV tracing to each of the DAergic, GABAergic, and glutamatergic neuron populations in the VTA using principal component analysis, PCA. A schematic of the experiment is shown in Figure 1A, with representative images of brain regions of interest in Figure 1B. A standard representation of these data, a bar plot, is shown in Figure 1C. PCA embedding resulted in the three cell populations clustering into largely distinct groups (Figure 1E). We considered the first two principal components (PCs), which explain the most variance in the data. In PC space, DAT-Cre brains were characterized by low PC1 and low PC2 values, GAD2-Cre brains by high PC1 values, and vGluT2-Cre brains by high PC2 values. The differences we noted were statistically significant: PC1 values of GAD2-Cre compared to DAT-Cre brains were significantly different (One-way ANOVA p = 0.010, DAT-Cre vs. GAD2-Cre multiple comparisons adjusted p = 0.0087; Figure 1F) and PC2 values of DAT-Cre compared to vGluT2-Cre brains were significantly different (One-way

ANOVA p = 0.016, DAT-Cre vs. vGluT2-Cre multiple comparisons adjusted p = 0.016; Figure 1G). An advantage of PCA over other dimensionality reduction approaches is that, because it is a linear approach, the weights of each feature’s contribution to each PC can be clearly interpreted. In this case, we can infer which brain regions contribute most to the positive and negative dimensions of each PC. For example, in our data, DAT-Cre and GAD2-Cre brains were separated along PC1. Brain regions that contributed to the positive dimension of PC1 were the BNST, entopeduncular nucleus (EP), GPe, central amygdala (CeA), and zona incerta (ZI) (Figure 1D). These inputs were enriched onto GABAergic cells relative to DA cells (Figure 1C). To identify where in the VTA these input populations projected, we assembled projection portraits of these inputs to the VTA as we did previously³². This approach uses axonal projection data available from the Allen Mouse Brain Connectivity Atlas to generate spatial heatmaps that show “hotspots” of innervation in the midbrain³⁵. We first selected experiments where GFP was expressed in the input regions of interest. Pixel values representing projection density were pulled from a representative coronal slice containing the VTA for each region. For a given set of inputs of interest, we averaged together their projections to the VTA to see if the group had any common spatial trends (Figure 1H). The projection portrait of the positive PC1 inputs indicated that they preferentially project dorsally and laterally to areas surrounding the VTA (Figure 1I). In contrast, the regions that contributed most strongly to the negative dimension of PC1 were the anterior cortex, nucleus accumbens medial shell (NAcMed), ventral pallidum (VP), preoptic area (PO), septum, lateral habenula (LHb), medial habenula (MHb), and lateral hypothalamus (LH). These inputs were more highly labeled in DAT-Cre mice. The projection portrait of these inputs highlights the whole VTA and not surrounding regions (Figure 1J). While the distribution of starter cells in the VTA in DAT-Cre and GAD2-Cre mice is spatially overlapping⁷, this analysis indicates that inputs with biases onto DA cells target the VTA itself while those biased onto GABAergic cells preferentially target areas dorsal and lateral to the VTA.

While DAT-Cre and GAD2-Cre conditions were separated along PC1, DAT-Cre and vGluT2-Cre conditions were separated along PC2. Regions that contributed to the positive dimension of PC2 included the LHb, MHb, EP, ZI, and deep cerebellar nuclei (DCN) as well as the cortex, while those contributing to the negative dimension included the NAcMed, nucleus accumbens lateral shell (NAcLat), nucleus accumbens core (NAcCore), dorsal striatum (DStr), and GPe (Figure 1K-L). The MHb and LHb both project medially in the VTA, whereas the NAcMed, NAcLat, NAcCore, DStr, and GPe all project laterally³². The projection portrait corresponding to the positive dimension of PC2 highlights the area dorsal of the VTA, as well as the medial aspect of the VTA (Figure 1K), while those corresponding to the negative dimension highlight the ventrolateral aspect of the ventral midbrain, including the substantia nigra pars reticulata (SNr) (Figure 1L). These results are consistent with the medial-lateral projection gradient that principally divides the projections into the VTA, with inputs from the basal ganglia projecting to the lateral portion of the VTA and the habenula projecting to the medial portion³². As inputs from the MHb and LHb contribute positive weights to the PC2 axis and vGluT2-Cre brains have the most positive PC2 values, it suggests that most of the vGluT2-expressing cells in the VTA are located within the medial portion of the VTA. This is consistent with previous reports using immunohistochemistry and in situ hybridization^47,48. In addition, the EP, ZI, and DCN all project predominantly to the red nucleus, and the cortex also innervates the red nucleus, which contains glutamatergic neurons and is located dorsal to the VTA (Figure 1K). While injections for our viral mapping studies targeted the VTA and most starter neurons were within the VTA, there were some starter neurons in vGluT2-Cre brains located in the red nucleus. In contrast, this was not the case in DAT-Cre brains¹⁰. As these starter populations were also largely spatially overlapping, the data indicate that inputs with biases onto glutamatergic cells target the dorsomedial aspect of the VTA, while those with a relative bias onto DA cells preferentially target the ventrolateral aspect of the VTA.

While PCA is particularly useful due to the interpretability of the weights of the PCs, it cannot capture non-linear relationships in high-dimensional data, which are common in complex systems like in neuroscience. To examine the overall relationship of input patterns between DAT-Cre, GAD2-Cre, and vGluT2-Cre cell populations, we performed a second analysis on inputs from each population using Uniform Manifold Approximation and Projection, UMAP, which can capture more complex relationships than PCA. We found that DAT-Cre and GAD2-Cre brains were the most separated in UMAP space, with vGluT2-Cre brains located between them (Figure 1M). This is consistent with the overlap in neurotransmitters released by these cells: cells expressing DA markers such as tyrosine hydroxylase (TH) and traditional GABAergic markers such as GAD2 are nearly mutually exclusive, while there is considerable overlap between cells expressing DA markers and glutamatergic markers, as well as cells expressing GABAergic markers and glutamatergic markers^{10,12,49–51}. Note that while DA neurons have been found to co-transmit DA and GABA^52,53, this occurs through a non-canonical GABA synthesis pathway, and therefore DA cells do not express typical markers of GABAergic cells, such as GAD2, used to define the Cre line in our study.

Since individual UMAP embeddings are stochastic depending on initial seeding conditions, we examined the average Euclidian distance between pairs of brains in 20 UMAP embeddings and plotted the averaged results as a correlogram. This showed that DAT-Cre and vGluT2-Cre brains were more like one another, and the GAD2-Cre brains were the most distinct from the other conditions (Figure 1N). These results are thus consistent with the known neurochemical overlap between these cells and indicate that DAT-Cre and vGluT2-Cre label more overlapping cells with each other than with GAD2-Cre. In sum, these results serve to demonstrate that we can recapitulate the locations of different neurochemically-defined cells in the VTA based on PCA and projection portraits based on information from RABV input mapping data. They also provide a validation of our approach for dissecting the relationships between input patterns and different cell types within the VTA, as the conclusions about the relationship between input regions and cell types agree with other studies using different, complementary methodologies.

A single injection of a drug of abuse changes brain-wide input patterns to VTA^DA cells

With our dimensionality reduction approaches validated on our cell type-specific RABV data, we next applied them to examining how drugs of abuse may alter VTA input patterns. We first assessed whether the psychoactive but non-addictive fluoxetine altered brain-wide input patterns relative to saline-injected mice. We found no clear differences in inputs between these conditions, and they overlapped completely in PC space – with no significant differences between PC1, PC2, or PC3 between the saline and the fluoxetine brains (p > 0.05 for PC1-3) - indicating fluoxetine did not alter brain-wide input patterns to VTA^DA cells (Figure 2A-F). We therefore combined these groups for additional statistical power and tested whether a single administration of one of several drugs of abuse caused changes in brain-wide input patterns to VTA^DA cells, first considered as a single group. When comparing combined controls to the drugs group, we observed a statistically significant difference in input patterns, including a reduction in inputs from the DStr and elevation in inputs from the GPe following drug injection (Figure 2G). In order to more fully assess which input sites may differentiate the control-versus drug-treated groups, we again used PCA. We observed separation between the control and drug groups principally along PC2, with no separation along PC1 and minimal separation along PC3 (Figure 2H-I). These differences in PC space were statistically significant: as we observed visually, the control group’s PC2 and PC3 values were significantly different from the drug group’s PC2 and PC3 values (p = 0.0009 and p = 0.024 for PC2 and PC3, respectively; Figure 2K-L). There was no significant difference between drugs and controls in PC1 (p = 0.78; Figure 2J). The positive axis of PC2, which was enriched in drug-treated brains, was driven by inputs from EAM, EP, GPe, LH, and PBN, while the negative axis was driven by the NAcMed, NAcCore, DStr, septum, and MHb (Figure 2M). Notably, the most negative weights for PC2, namely the NAcMed, NAcCore, DStr, septum, and MHb are all reduced in drug-treated brains (Figure 2G). Conversely, the most positive weights, namely the EAM, EP, GPe, LH, and PBN, are all increased in drug-treated brains (Figure 2G). Therefore, while PCA allows us to observe trends consistent with traditional methods, we can also observe the relationship between the percentage of inputs among different brains, and also identify relationships that are not immediately apparent in the data when presented as bar graphs. For example, while the BNST, CeA, and ZI do not show clear differences in control vs. drug-treated brains in the bar graphs, they are all positive weights in PC2 and thus, likely change in concert with regions such as the GPe and PBN (Figure 2G-M).

Importantly, unlike traditional bar graph approaches, PCA allows us to isolate technical variation from the experiments, enabling us to focus more clearly on changes induced by the relevant variable, here treatment with a control substance or drug of abuse. For example, PC1 does not distinguish control vs. drug-treated brains (Figure 2H, J), indicating this variable was inherent to the way the experiment was performed and not related to drug or control treatment. The strongest negative weights for PC1 were the NAcLat, NAcCore, DStr, and GPe, while the strongest positive weights for PC1 were the VP, PO, LHb, and DR (Figure 2M). This is consistent with connectivity differences along the medial-lateral axis of the VTA, which underlies many of the connectivity differences to different VTA^DA subpopulations^7,10,32. To assess whether PC1 may reflect differences the location of starter cells in the VTA, we first assessed whether there was a systematic difference in starter cell location between conditions. We noted a high level of overlap in starter cell distribution between conditions (Figure 2N), indicating that this alone likely did not explain differences between control and drug-treated groups. We then assessed whether brain to brain variation in starter cell distribution could explain the effects observed along PC1. We therefore colored each point in PC space according to the assessed center of injection for each experiment. We noted a clear medial-lateral gradient along PC1, but not PC2 or PC3 (Figure 2O-P). To test the significance of this relationship, we performed linear regression analysis, comparing the x-coordinate of the injection center relative to location along PC1, PC2, or PC3 for each brain. We observed a highly significant positive correlation with PC1, a slightly positive but non-significant correlation with PC2, and no correlation for PC3 (Figure 2Q-S). These results show that PC1 reflects variation in the injection site within the VTA, while PC2 reflects changes induced by drugs of abuse, allowing us to isolate the changes in input connectivity induced by administration of a single drug of abuse.

With our integrative methodology validated on cell type specific RABV data and demonstrating that drugs of abuse alter brain-wide input patterns to VTA^DA cells, we applied our approaches to explore the effects of each drug specifically, as well as potential effects of the 2x use of the K/X anesthetic mixture. We first asked two questions: 1) How do injections of different drugs of abuse change the connectivity patterns onto VTA cells? 2) How do these potential drug-induced differences relate to differences between unique cell types in the VTA? To answer these, we first assessed experiments with ten separate experimental conditions: 2 conditions of GAD2-Cre mice, and 8 conditions of DAT-Cre mice. The GAD2-Cre mice were all anesthetized with isoflurane and administered either saline or cocaine one day prior to RABV injection. The DAT-Cre mice were given either isoflurane or K/X-based anesthesia. The K/X-anesthetized mice were administered saline. The isoflurane anesthetized mice were administered either saline, fluoxetine, amphetamine, cocaine, methamphetamine, morphine, or nicotine. RABV mapping was performed from targeted VTA cells in all cases, and the brain-wide inputs were quantified. We then performed PCA on the combined RABV tracing datasets, which showed that brains from GAD2-Cre mice were in the positive dimension of PC1, while those from DAT-Cre mice were in the negative dimension (Figure 3A), regardless of the substance administered. This difference in PC1 was highly statistically significant (One-way ANOVA p < 0.0001; all conditions pairwise comparisons with GAD2-Cre, multiple comparisons adjusted p < 0.0001; Figure 3B). The positive dimension of PC1 received contributions from largely the same sets of regions as in the previous analysis of DAT-Cre, GAD2-Cre, and vGluT2-Cre brains, which included the BNST, EP, ZI, CeA, as well as the DCN (Figure 3H, 1D). Given that the same set of brain regions differentiated DAT-Cre and GAD2-Cre brains along PC1 in both datasets, this suggests that the intrinsic differences in inputs to different cell populations in the VTA are more pronounced than differences induced by drugs of abuse, fluoxetine, or the method of anesthesia. Notably, these patterns were unlikely to be obtained by chance; scrambling the sample identity and percentage of RABV-labeled cells in each input site resulted in groups that overlapped in both PCA and UMAP space (Figure S1).

To specifically explore potential drug-induced differences in input labeling onto VTA^DA cells, we removed GAD2-Cre brains from the dataset. Results and representative images are shown in Figure 3J-K. In PCA analyses of these brains only, PC1 and PC3 largely failed to differentiate the conditions, with no significant differences observed (Figure 3C-G), while PC2 significantly differentiated the conditions (One-way ANOVA p = 0.0093; saline vs. drugs multiple comparisons adjusted p = 0.034; Figure 3F). The most negative weights to PC1 were contributed by the basal ganglia inputs that project to the SNr, located ventrolaterally to the VTA (Figure 3I, L) while positive contributors came from brain regions that projected more medially (LHb, MHb, dorsal raphe (DR)) and centrally to the VTA (ventral pallidum (VP), PO, extended amygdala area (EAM), paraventricular nucleus of the hypothalamus (PVH), and LH) (Figure 3I, M). This is consistent with a gradient along the medial-lateral axis³². As it did not differentiate the conditions, this axis likely represents experimental variations in the location of starter cell populations within the VTA^7,10,32. The regions with the largest positive contribution to PC2 included the EAM, parabrachial nucleus (PBN), DR, EP, and LH, with lesser contributions from the GPe, DCN, ZI, CeA, and BNST (Figure 3I, N). These regions were over-represented in animals treated with a single administration of a drug of abuse. Notably, this group of inputs contains brain regions that are involved in the brain’s stress response such as the PBN, EAM, CeA, and BNST, and have been implicated in drug-induced behaviors^54,55. For example, we recently showed the importance of the GPe in conditioned place preference and locomotor sensitization and the BNST in drug-induced withdrawal^15,25. The regions with the largest negative contribution to PC2 included the NAcMed, NAcCore, DStr, Septum, and MHb (Figure 3I, O). These regions were underrepresented in animals given a single administration of a drug of abuse. These results together demonstrate that drugs of abuse cause long-lasting changes in input populations to DA cells, and that we can identify which sets of inputs distinguish drug-treated vs. control-treated brains. These changes may be functionally important, as we previously demonstrated the functional role of two inputs that we identified as differing between drug-treated and control conditions.

Lastly, we wanted to explore the overall relationship between input changes induced by different drugs. As before, we used UMAP to explore the higher-level organization of these conditions to compare the relationship of holistic input patterns between groups. We found that the mice anesthetized with isoflurane and treated with saline overlapped with the isoflurane-anesthetized, fluoxetine-treated group, but not the K/X-anesthetized and saline-treated group (Figure 3P). These conditions were segregated from animals treated with the drugs of abuse - psychostimulants, morphine, and nicotine. Our analysis of the relationships between individual conditions, generated using the Euclidian distances between points in 20 UMAP embeddings as before, demonstrated that the expected control conditions formed a cluster (saline, fluoxetine) distinct from the other groups. Interestingly, the K/X-anesthetized animals were more closely associated with the mice treated with addictive drugs, in particular morphine, than with controls (Figure 3Q). This suggests that drugs of abuse cause similar changes in brain-wide input patterns to one another in DAT-Cre mice, and that input patterns in K/X-anesthetized animals more closely resemble morphine-treated mice than saline-treated mice. This indicates that the method of anesthesia, either using isoflurane or K/X, results in different patterns of input labeling. Given the similarity of input patterns in K/X-anesthetized mice to isoflurane-anesthetized mice treated with a drug of abuse but not saline, we suspect that two doses of K/X-anesthesia as is required for RABV circuit mapping experiments triggers changes in input patterns in rodent brains that last at least the five-day duration of these experiments.

Anesthetic doses of K/X cause long-lasting changes in inputs to VTA^DA®NAcLat and VTA^DA®Amygdala cells

To examine how differences in anesthesia methods may impact brain-wide input patterns, we performed some more targeted comparisons. When comparing only two groups using PCA – isoflurane-anesthetized animals to K/X-anesthetized animals – these conditions completely segregated (Figure 4A). While the difference in PC1 was not quite significant (p = 0.086), the difference in PC2 was significant (p = 0.037). This separation was driven by positive weights from the VP, PO, EAM, PVH, LH, and PBN along PC1, which were more strongly represented in the K/X-anesthetized mice and whose combined projection portrait broadly targets the VTA (Figure 4A, C, P, U), and the NAcLat, DStr, GPe, and PVH along PC2, which were also more strongly represented in the K/X-anesthetized mice whose combined projection portrait targets ventrolateral to the VTA (Figure 4A, C, P, V). A similar separation of conditions was observed using UMAP (Figure 4B).

Next, we wanted to assess which VTA^DA cells may be impacted by whether the animals received isoflurane or K/X-based anesthesia. To test this, we performed cell type-specific Tracing of the Relationship between Inputs and Outputs (cTRIO) for each of four different VTA^DA subpopulations: those projecting to the NAcMed, NAcLat, amygdala, or mPFC, and receiving K/X anesthesia and saline, isoflurane and saline, or isoflurane and cocaine. The differences in input patterns between conditions were most pronounced for VTA^DA®NAcLat cells, where the isoflurane- and K/X-treated groups segregated, both in PCA and UMAP space (Figure 4D-F). Groups appeared to differ along both PC1 and PC2, though differences in PC2 were not statistically significantly different, with K/X-treated brains being characterized by lower PC1 values which received contributions principally from the anterior cortex, NAcMed, NAcLat, NAcCore, septum, PVH, and MHb, as well as higher PC2 values which received contributions from the VP, PO, EAM, LHb, and LH (PC1, One-way ANOVA p = 0.0003, saline vs. cocaine multiple comparisons adjusted p = 0.011, k/x vs. cocaine multiple comparisons adjusted p = 0.0002; Figure 4D, Q, W, X). The following comparisons in PC1 were statistically significant: cocaine versus isoflurane controls (p = 0.011) and cocaine versus K/X anesthesia (p = 0.001). Notably, K/X, saline-treated mice and isoflurane, cocaine-treated mice segregate from one another in both PCA and UMAP space, indicating that the K/X-induced changes are distinct from those induced by cocaine (Figure 4D-F). The similarity in projection portraits for weights that contribute more to the K/X-treated group in VTA^DA®NAcLat (Figure 4U, V) and all VTA^DA neurons (Figure 4W, X) suggests that largely overlapping DA cells may be targeted for each set of experiments and is consistent with the observation that VTA^DA®NAcLat neurons comprise the majority of VTA^DA neurons, as noted in previous studies^7,10. These results thus serve as another proof of principle of the reliability of our approach to pull out similar results from different datasets.

In addition to VTA^DA®NAcLat cells, we observed some segregation of isoflurane and K/X-treated conditions in VTA^DA®Amygdala cells (Figure 4G-I, R). The lack of clear differentiation in PCA space, with only PC3 showing significant differences between the isoflurane control group and K/X anesthesia (One-way ANOVA p = 0.014, saline vs. K/X multiple comparisons adjusted p = 0.0078; Figure S2) makes it difficult to interpret which input sites may be contributing to the differentiation in conditions that we observe. In the VTA^DA®NAcMed and VTA^DA®mPFC subpopulations, the separation between isoflurane-saline, isoflurane-cocaine, and K/X-saline conditions was much less pronounced (Figure 4J-O, S-T), and not to be significant; indeed, similar separations were observed by chance (Figure S3). Notably, when considering inputs to all VTA^DA cells, the projection portraits that indicate the inputs most highly elevated or depressed in K/X-treated mice highlight either the entire VTA, or the ventrolateral aspect of the VTA (Figure 4U, V). These projections are consistent with the locations of the VTA^DA®Amygdala and VTA^DA®NAcLat cells, respectively, within the VTA. In contrast, the portraits do not highlight the medial or ventromedial aspects of the VTA, which is the location of VTA^DA®mPFC and VTA^DA®NAcMed cells^32,56. These results together indicate that K/X anesthesia impacts inputs that preferentially target VTA^DA®NAcLat or VTA^DA®Amygdala cells.

Exploring gene expression patterns that predict changes in RABV input labeling

We next wanted to explore if drug-induced changes in input connectivity might be related to gene expression profiles in cells located in the input sites. The AGEA data are from all cells in the input regions and are not limited to those projecting to the VTA and show only baseline gene expression patterns and not those following drug administration. Nonetheless, these analyses can be a valuable starting point in associating spatial gene expression data with patterns of RABV-labeled inputs from VTA^DA cells. To do this, we calculated the correlation of gene expression data from the Allen Gene Expression Atlas (AGEA) with the difference (Figure 5) or the ratio (Figure 6) between inputs labeled in each region in experimental conditions and isoflurane/saline-treated controls. After filtering AGEA data for only high-quality experiments, we were left with 4,283 genes to which we could calculate correlations within our 22 input brain regions of interest.

We first identified genes or pathways that may be linked to elevations or depressions in RABV labeling for each condition. We selected genes that were highly positively or negatively correlated to differences between 1) isoflurane-saline controls and 2) isoflurane-psychostimulant, morphine, nicotine, or fluoxetine, or 3) K/X-saline mice. We used the top 50 genes in the AGEA that were either positively or negatively correlated (Spearman’s ρ) to differences in RABV labelling and performed gene ontology (GO) analyses. Starting with drugs of abuse in aggregate, the 50 most positively correlated genes (Spearman’s ρ > 0.47, p < 0.05) were most significantly associated with inhibitory ligand-gated chloride channels (Figure 5A), suggesting that higher basal expression of these genes may be related to larger changes in connectivity induced by drugs of abuse. Conversely, there were many more genes that were strongly negatively correlated to drug-induced changes in RABV labeling. The 50 most negatively correlated genes (Spearman’s ρ < - 0.68, p < 0.001) were most closely associated with negative regulation of long-term synaptic potentiation (Figure 5B), suggesting that regulation of synaptic plasticity in input cell populations may be a key indicator of the likelihood of a drug of abuse altering RABV input labeling. We repeated this analysis of the top 50 genes with each individual drug: cocaine, methamphetamine, amphetamine, morphine, and nicotine. In addition to findings in aggregate, other interesting classes of genes that were negatively associated (Spearman’s ρ < -0.50, p < 0.05) with individual drug-induced changes in RABV labeling include regulation of synapse structural plasticity (methamphetamine), regulation of calcium ion transport (nicotine), cellular response to Ca²⁺ ions (cocaine), and regulation of potassium ion transporters (nicotine) (Figure S4), which relate to synaptic and cellular processes of plasticity and intrinsic excitability. There were fewer than 50 genes that were significantly positively correlated, but the top 50 genes (Spearman’s ρ > 0.36, p < 0.12) were associated with terms including anion channels (cocaine, amphetamine, methamphetamine) and syntaxin binding (cocaine, amphetamine, nicotine) (Figure S5). We also explored which gene expression patterns may be related to input changes driven by K/X anesthesia. Anion channel activity, particularly GABA receptor expression, again shows up as one of the most significant associations to the top 50 most positively correlated genes (Spearman’s ρ > 0.37, p < 0.10) (Figure 5C), while negatively correlated genes were associated with a variety of metabolic processes (Figure 5D). Notably, these gene expression patterns were not correlated with input changes induced by fluoxetine (Figure S6).

To further assess which classes of genes were the most closely related to the differences in drug evoked RABV labeling relative to controls, we performed linear regressions of the ratios of RABV labeling in experimental vs. control conditions relative to the mean of the log normalized expression of AGEA RNA in situ hybridization data for a particular class of genes. We observed an overall significant and slightly negative correlation between the RABV labeling ratio and the gene expression mean (r = -0.30, p = 2.1 x 10^-3) for all 4,283 genes (Figure S7A). This genome-wide negative correlation is a useful reference to compare to correlation coefficients of subsets of genes. Among the gene expression classes that we examined, the most significant association was between the change in RABV-labeled inputs and ion channel gene expression (r = -0.56, p = 1.3 x 10^-9; Figure 6A), consistent with our results from GO analyses. The next most significant classes were neurotransmitter receptors (r = -0.46, p = 1.9 x 10^-6) and synapses (r = -0.42, p = 1.6 x 10^-5), followed by genes that have been linked to SUDs through Genome-Wide Association Studies (GWAS; r = -0.33, p = 7.7 x 10^-4) as well as endo/exocytosis genes (r = -0.30, p = 2.3 x 10^-3). Each of the above relationships showed a negative correlation of RABV ratios with gene expression (Figure 6B-C, S7B-C). Except for endo/exocytosis genes, each class had a more negative, significant correlation than our genome-wide reference. As a further comparison, we observed that the relationships between RABV labeling ratios and gene classes that are likely irrelevant to RABV labeling, such as mitochondrial genes (r = -0.24, p = 1.4 x 10^-2), sugar transporters (r = -0.26, p = 1.0 x 10^-2), and tumor-associated genes (r = -0.18, p = 7.0 x 10^-2) were weaker and less significant than our classes of interest (Figure S7D-F).

As the class of ion channel genes was the most strongly correlated to drug-induced differences in RABV-labeled input patterns, we assessed which ion channels are most closely linked to RABV labeling differences between conditions. The most significant relationship was between Ca²⁺ channel genes (r = -0.57, p = 4.4 x 10^-10), followed by Cl^- (r = -0.39, p = 6.1 x 10^-5), K⁺ (r = -0.30, p = 2.9 x 10^-3), and Na⁺ channels (r = -0.11, p = 0.30; Figure 6D-G). The expression levels of these channels, in particular Ca²⁺ channels, are linked to regulating a variety of cellular processes, including immediate early gene expression networks as well as cellular excitability. Ligand-gated channels (r = -0.40, p = 4.1 x 10^-5) were more strongly related with RABV labeling than voltage-gated ion channels (r = -0.33, p = 7.4 x 10^-4) or other types of ion channels (r = -0.13, p = 0.21; Figure S7G-I). Among the ligand-gated channels, the glutamate receptors were the most strongly correlated (r = -0.49, p = 3.1 x 10^-7; Figure 6H). The next most strongly correlated gene classes were acetylcholine receptors (r = -0.36, p = 2.4 x 10^-4), GABA receptors (r = -0.31, p = 1.8 x 10^-3), and glycine receptors (r = 0.32, p = 1.4 x 10^-3) (Figure 6I-K). Notably, glycine receptors were the only class of genes that we analyzed that was positively correlated with RABV labeling. Similar patterns were observed for K/X anesthetized mice (Figure S8), but not with isoflurane-anesthetized, fluoxetine-treated mice (Figure S9). These results together indicate that the basal expression levels of a variety of ion channels, particularly Ca²⁺ and ligand-gated ion channels, may be related to driving elevations or depressions in RABV labeling following a single injection of a variety of drugs of abuse, or K/X anesthesia.

As the correlation between RABV labeling ratios and Ca²⁺ channel gene expression is negative, lower expression of Ca²⁺ genes is linked to an elevation in drug-induced RABV input labeling. This relationship could be explained by several mechanisms. One, expression of Ca²⁺ channel genes could be inhibitory towards the drug-induced change in RABV transmission. Alternatively, lower Ca²⁺ channel gene expression levels could indicate more potential for a drug-induced rise in gene expression level. To differentiate between these possibilities, we performed CRISPR-based gene expression modulation of the voltage-gated calcium ion channel Cacna1e using CRISPRi, which reduces gene expression through blocking DNA transcription without cutting DNA⁵⁷ (Figure 7A-B). Three different adeno-associated viruses (AAV) were co-injected into the NAcLat, an input site to the VTA, in DAT-Cre mice. These viruses expressed 1) Cre, 2) Cre-dependent dCas9, 3) Cre-dependent small guide RNA (gRNA) targeting Cacna1e. A two AAV mixture lacking the AAV expressing the gRNA was used as a control. The three AAV experimental mix was administered into the NAcLat in one hemisphere, while the two AAV control mix was injected into the other hemisphere. During the same surgery, AAVs expressing Cre-dependent TVA-mCherry and rabies virus glycoprotein were bilaterally injected into the VTA, so that RABV mapping could be done in both hemispheres (Figure 7C-D). Because the NAcLat, NAcCore, and GPe unilaterally project to the VTA in the same hemisphere, and not the contralateral VTA^7,10,15,25, the ratio of labeled inputs in these sites on a single hemisphere can be used to assess the relative labeling in the NAcLat in experimental and control conditions within each mouse.

Figure 7:

We first validated that our CRISPRi strategy reduced Cacna1e gene expression. We designed our gRNA to cause only a mild reduction in Cacna1e expression to more naturally mimic potential changes in gene expression. We found that our CRISPRi perturbation reduced gene expression by ∼20% (Figure 7E). This knockdown of Cacna1e in the NAcLat via CRISPRi resulted in decreased RABV-labeled inputs from the NAcLat compared to control mice where the shRNA was not included (Figure 7F). These results indicate that reduced levels of Cacna1e lower the number of RABV-labeled inputs from the NAcLat, suggesting that expression levels of Cacna1e are directly related to RABV-labeled inputs. This means that expression of this gene is not inhibitory towards RABV transmission. Therefore, a more likely hypothesis is that a low basal gene expression of Ca²⁺ channels allows for greater changes in the expression of these genes following drug exposure.

Discussion

Here we explored the long-term consequences of a single injection of a variety of abused drugs on the connectivity of VTA^DA cells. We first confirmed that dimensionality reduction approaches could adequately differentiate the connectivity between different cell populations in the VTA, and showed which inputs differentiate DAergic, GABAergic, and glutamatergic cells in the VTA (Figure 1). We then showed that a variety of drugs of abuse as well as K/X-based anesthesia trigger long-lasting changes in the brain-wide input patterns to VTA cells, which were clearly observable though less substantial than the differences onto different cell types or natural variation in starter cell location (Figures 2-3). We used our cTRIO method to identify DA cell type-specific changes induced by cocaine as well as K/X anesthesia, implicating the VTA^DA®NAcLat and VTA^DA®Amygdala cells as those exhibiting the most substantial input changes (Figure 4). We then explored intrinsic gene expression patterns throughout the brain that relate to the elevation or depression in RABV labeling by different drugs, finding that expression levels of ion channels, neurotransmitter receptors, and synapse-related genes were all significantly negatively correlated with RABV input labeling ratios, with Ca²⁺ channels and glutamate receptors being the most strongly correlated (Figure 6). We then demonstrated the causative role of one calcium channel gene, Cacna1e, in influencing RABV input mapping efficiency (Figure 7).

Drugs of abuse induce changes in connectivity onto VTA^DA neurons from common and distributed input populations throughout the brain

Using dimensionality reduction approaches we differentiated the input patterns to different VTA cells (Figure 1) and showed that these differences were larger than the differences in inputs from the same cell populations in animals given a drug of abuse vs. saline (Figure 3A). Notably, these differences were unique from those triggered by the psychoactive but non-addictive substances fluoxetine and K/X (Figure 3C, F). Differences were apparent on the 2^nd PC axis, which were driven by positive weights from the GPe, EP, EAM, LH, ZI, DR, and PBN, regions that showed elevated connectivity in drug-treated mice. Many of these regions have been implicated in cocaine addiction^58–61. While these drug-induced differences in input labeling were less substantial than intrinsic cell type connectivity differences (Figure 3A) or differences in starter cell location (Figure 2N-S), these differences are more apparent when using dimensionality reduction methods such as PCA and UMAP. That different cell types exhibit larger differences in input patterns than the same cell type with and without drug administration is perhaps not surprising, as different populations of cells may receive different sets of inputs from one another, and these are likely more stable across experiments than the differences in connectivity triggered following a single exposure to a drug of abuse. The same is likely true for spatial gradients, which we determined to be a major source of variation in RABV inputs to the VTA³². Notably, the minor differences in starter cell location are likely not a technical concern specific to our study, but rather inherent to all RABV mapping experiments that rely on intracranial injections, as a small amount of biological or technical variation is to be expected. The PCA approach we employ here, however, can better isolate these effects than traditional bar graph approaches and thus enable clearer observation of the effect of the relevant variable being tested, here effects of a single drug administration (or K/X anesthesia) on the input connectivity to VTA^DA cells.

The results in this study further confirm our previous observations that a single exposure to a variety of drugs of abuse triggers long-lasting changes in inputs to VTA^DA neurons. Our interpretation of these data is that the changes are largely being driven by a change in the activity of input cells, as we have shown using fiber photometry and slice electrophysiology in three separate studies^15,16,62. However, we have not confirmed this to be true for all cases. As a single exposure to a variety of drugs of abuse alters the strength of both excitatory and inhibitory inputs in the VTA^5,20–22, we cannot rule out that changes in synapse number or strength may also influence the observed connectivity results.

Ketamine/xylazine anesthesia alters the input landscape of VTA^DA cells

Animals anesthetized with a K/X mixture showed a different pattern of input labeling to VTA^DA cells than animals anesthetized with isoflurane (Figure 4A-B). While we cannot perform RABV mapping on unanesthetized mice, making it impossible to definitively define which mode of anesthesia, or both, alters RABV input patterns we suspect that isoflurane anesthesia has minimal effect on the brain-wide inputs to VTA^DA cells. This is largely because animals treated with a drug of abuse with isoflurane anesthesia look more like those receiving K/X-saline than isoflurane-saline animals (Figure 3). In addition, ketamine is a drug with known abuse potential^63–66. While the anesthetic doses used here (100 mg/kg) are well above the doses used recreationally (∼1 mg/kg), the changes that we observed with K/X anesthesia may be related to the abuse potential of ketamine. Interestingly, the K/X-induced input changes were most similar to those induced by 10 mg/kg of morphine (Figure 3P-Q). Antidepressant doses of ketamine may work in part via the opioid system; indeed, it has been known for some time that ketamine binds the mu opioid receptor^67–70, and it was more recently shown that the opioid antagonist naltrexone blocks the antidepressant effects of ketamine^71,72. Given that anesthetic doses of ketamine (100 mg/kg) are roughly 10 times higher than antidepressant doses of ketamine (∼10 mg/kg), K/X anesthesia almost certainly engages the opioid system.

Relationship of gene expression with the extent of drug-induced RABV input labeling

Gene expression patterns in the mouse brain are correlated with the extent of elevations and depressions of RABV labeling induced by one of several drugs of abuse. We observed a slightly negative correlation of RABV labeling changes with the overall level of gene expression across the brain. We are not sure of the implications of this negative correlation, though it could suggest that if initial gene expression levels are lower, it leaves a higher potential for experience-induced elevations. We observed significant negative correlations of RABV labeling changes induced by drugs of abuse with the levels of expression of ion channel genes, neurotransmitter receptor-related genes, and synapse-related genes. These results and our experiments demonstrating the functional impact of downregulating Cacna1e for RABV input labeling (Figure 7) indicate that the expression levels of these genes could be used as predictive biomarkers of which regions are likely to undergo drug-induced changes in connectivity to VTA^DA cells. They also suggest that the expression levels of these gene classes are related to the drug-induced changes in RABV labeling.

While the AGEA is a useful resource for relating spatial gene expression information to our RABV mapping datasets, there are notable limitations. First, the AGEA dataset is not as quantitatively rigorous as modern single nucleus RNA sequencing or spatial transcriptomic datasets. Nonetheless, we believe that the overall patterns of high and low expression, especially when considered among classes of genes as was performed here, can still be informative. Furthermore, the AGEA dataset provides a static picture of gene expression levels throughout the brain, and thus does not account for potential changes in gene expression levels triggered by administration of one of the drugs used in this study. Third, the AGEA dataset considers expression levels in all cell types in a brain region, irrespective of their individual transcriptomic/molecular identities, while only a subset of cells in a region project to VTA^DA cells. Nonetheless, our CRISPRi experiment showing that a mild reduction of Cacna1e reduces RABV labeling indicates the potential for causative associations between gene expression patterns and drug-induced RABV input labeling changes (Figure 7). This result also directly links cellular excitability and the extent to which RABV labels these cells as monosynaptic inputs to defined cell populations.

Our RABV screen identified several input sites that show elevated or depressed labeling following a single administration of a drug of abuse, suggesting potential changes in excitability of these cells. Ca²⁺ channel genes were also highly negatively correlated with the extent of drug-induced changes in RABV input labeling. Ca²⁺ channels have been implicated in substance abuse through Genome-Wide Association Studies⁷³, and molecules that target these channels have shown some promise as addiction therapeutics^74–79. It is of particular interest that Ca²⁺ and not Cl^-, Na⁺, or K⁺ channels are highly correlated to RABV input labeling changes. Ca²⁺ is critically involved as a secondary messenger in neurons that affects a variety of processes including long-lasting changes in transcription and synaptic transmission, whereas the other ions are mostly involved only in facilitating or inhibiting formation and propagation of the action potential. Understanding how Ca²⁺ channel gene expression levels relate to drug-induced changes in the brain will help to illuminate the relationship between basal gene expression levels, drug-induced changes in RABV labeling, and substance abuse, providing mechanistic insights into how inter-individual differences influence susceptibility to substance abuse.

Acknowledgements

We would like to thank Hemmings Wu and Zhoule Zhu for their thoughtful critique of our manuscript.

Author contributions

KTB performed the viral tracing experiments, KB performed the majority of computational analyses of the data, PD assisted with network constructions and projection portraits, MH assisted with figure preparation, GT, JJV, and GA performed CRISPRi experiments, and KB and KTB wrote the manuscript.

Funding

This work was supported by the NIH (R00 DA041445, DP2 AG067666, R01 DA054374, R01 DA056599, R01 NS130044), Tobacco Related Disease Research Program (T31KT1437, T31P1426), American Parkinson Disease Association (APDA-5589562), Alzheimer’s Association (AARG-NTF-20-685694), New Vision Research (CCAD2020-002), Brightfocus Foundation (A2022031S), and the Brain and Behavior Research Foundation (NARSAD 26845) to KTB, NIH T32GM008620 and NIH F30DA056215 to MH, T32 GM136624 to KB and PD, and NSF GRFP DGE-1839285 to PD.

Competing interests

The authors have nothing to disclose.

Note

This reviewed preprint has been updated to use the correct text; previously, a version prior to the one reviewed was presented here.