Stress resilience is associated with transcriptional remodeling in the VTA

Stress resilience is associated with transcriptional remodeling in the VTA | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results Stress resilience is associated with transcriptional remodeling in the VTA View ORCID Profile Adelaide R. Minerva , View ORCID Profile Brenna McMannon , View ORCID Profile Rixing Lin , View ORCID Profile Anna Zhukovskaya , View ORCID Profile Ilana B. Witten , View ORCID Profile Catherine Jensen Peña doi: https://doi.org/10.1101/2025.11.13.687606 Adelaide R. Minerva 1 Princeton Neuroscience Institute Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Adelaide R. Minerva Brenna McMannon 1 Princeton Neuroscience Institute Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Brenna McMannon Rixing Lin 1 Princeton Neuroscience Institute Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Rixing Lin Anna Zhukovskaya 1 Princeton Neuroscience Institute Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Anna Zhukovskaya Ilana B. Witten 1 Princeton Neuroscience Institute 2 Department of Psychology 3 Howard Hughes Medical Institute Princeton University , Princeton NJ 08544 USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Ilana B. Witten Catherine Jensen Peña 1 Princeton Neuroscience Institute Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Catherine Jensen Peña For correspondence: cpena{at}princeton.edu Abstract Full Text Info/History Metrics Supplementary material Preview PDF Summary Individual responses to chronic stress vary, with some individuals remaining resilient while others exhibit susceptibility. The ventral tegmental area (VTA), a region involved in reward learning, and the lateral habenula (LHb), a region involved in aversive learning, have been implicated in the pathophysiology of stress-related mood disorders. Here, we seek to understand the molecular adaptations in these regions at the level of single cells that mediate susceptibility versus resilience. In particular, it remains unclear whether, at the level of gene expression, different cell types within different brain regions mediate stress susceptibility versus resilience, or if these phenotypes are mediated by distinct trajectories within the same cell types. To address this gap, we performed single-nucleus RNA-sequencing of LHb and VTA of mice subjected to chronic social defeat stress. While we found minimal gene expression changes in the LHb after stress, the VTA exhibited widespread, cell type-specific transcriptional remodeling in resilient individuals and few gene expression changes in susceptible individuals. Across VTA cell types, resilience was associated with the coordinated upregulation of genes involved in intercellular signaling and neural communication, with maintenance of receptor-ligand interaction strength in resilience that was not present in susceptibility. Within VTA neurons, gene expression changes were most prominent in glutamatergic and dopaminergic clusters. Multivariate analyses of dopamine and glutamate subclusters showed that resilient neurons diverged more from control than susceptible neurons, but along a similar trajectory, supporting a model in which resilience reflects greater stress-related adaptations in these cell types. Together, these findings highlight the VTA as a key site of molecular plasticity in stress resilience and therefore a potential therapeutic target. Introduction Mood disorders are among the most pervasive and disruptive mental health disorders worldwide 1 , 2 . These illnesses are often triggered or exacerbated by chronic stress 3 – 7 . However, not all individuals who experience stress develop psychiatric symptoms 8 , and therefore an important scientific and public health goal is to identify the biological mechanisms that promote susceptibility or confer resilience. Chronic social defeat stress (CSDS) is a well-validated paradigm that captures individual differences in stress susceptibility in rodents 9 – 13 . In the CSDS model, a subset of animals exhibit social avoidance and anxiety-like behavior following repeated social defeat and are considered susceptible, while the remainder are characterized by a lack of overt behavioral changes compared to control and are considered resilient. Multiple brain regions likely mediate the heterogeneous effects of chronic stress. The lateral habenula (LHb) is an area that drives aversive learning 14 – 18 , is implicated in the pathophysiology of depression 19 – 29 , and shows stronger activity in susceptible animals during and following CSDS 30 . In contrast, the ventral tegmental area (VTA), a key node in the brain’s reward circuitry, is engaged in resilience 31 – 33 . During CSDS, resilient mice have higher dopamine (VTADA) activity during resilience-associated behavior, and stimulating VTADA neurons can increase the chance of resilience 32 . These changes underlying susceptibility or resilience could arise from shifts in gene regulation. Here, we seek to distinguish between multiple hypotheses regarding the structure of such shifts. The behavioral change in susceptible mice could be accompanied by underlying gene expression changes that support strong aversive learning (such as in the LHb), while resistance to behavioral change in resilient mice may reflect only minimal changes along the same molecular trajectory ( Figure 1A , Hypothesis 1). Alternatively, resilience may recruit more pronounced adaptations, but again in the same direction as susceptibility, (such as in the VTA) to counteract the effects of stress and maintain behavioral stability ( Figure 1A , Hypothesis 2). Additionally, within each region, susceptibility and resilience might represent divergent molecular trajectories away from control in different directions ( Figure 1A , Hypothesis 3). Under this hypothesis, the two phenotypes are not simply defined by the degree of transcriptional change in a given direction, but by fundamentally different gene expression programs. Download figure Open in new tab Figure 1. snRNA-sequencing in the VTA and LHb after chronic social defeat stress A. Hypotheses for how resilience and susceptibility may be instantiated in the brain at a transcription level. B. Schematic of chronic social defeat stress paradigm. C. During the social interaction test, time spent in the interaction zone (left) and corners (right) of the arena (1-way ANOVA for time in corners with behavior outcome as factors, F(2,69)=21.69, p=4.9e-08; Tukey HSD for control vs resilient p=1.0, control vs susceptible p=0.0, resilient vs susceptible p=0.0). Black outlines indicate which mice were used for sequencing. D. In the elevated plus maze, number of entries into the open arms (1-way ANOVA for open arm entries with behavior outcome as factors, F(2,69)=4.89, p=0.01; Tukey HSD for control vs resilient p=0.175, control vs susceptible p=0.008, resilient vs susceptible p=0.255). Black outlines indicate which mice were used for sequencing. E. Relationship between time in the interaction zone during the SI test and number of entries into the open arms of the EPM (Pearson correlation, r=0.247, p=0.036). Black outlines indicate which mice were used for sequencing. F. Schematic of tissue collection, library preparation, and snRNA-seq. See Table S1 for statistics. ∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001. Unless specified, data plotted as mean ± SEM. To test these hypotheses, we characterized transcriptomic adaptations in the LHb and VTA following chronic social stress. Cells in both regions are molecularly 34 – 43 and functionally diverse 15 , 17 , 21 , 40 , 44 – 82 and therefore it is likely that subpopulations may exhibit unique adaptations. Most prior work on gene expression changes underlying stress have used bulk tissue approaches 10 , 83 – 85 and thus would be unable to test these hypotheses directly. Therefore, unbiased study of all cell types in the VTA is necessary for a comprehensive understanding of stress-induced, phenotype-specific gene expression changes. To address this, we used single-nucleus RNA-sequencing (snRNA-seq) to profile transcriptional changes across cell types in both regions of mice resilient and susceptible to CSDS. More specifically, we test the idea that different regions and cell types primarily underlie stress susceptibility ( Figure 1A , Hypothesis 1), primarily underlie stress resilience ( Figure 1A , Hypothesis 2), or contribute to both processes by undergoing different types of changes ( Figure 1A , Hypothesis 3). Surprisingly, we found few gene expression changes in the LHb following stress. In contrast, in the VTA we found that (1) resilience exhibits more changes in gene expression than susceptibility, (2) these changes are most dramatic in glutamatergic and dopaminergic VTA neurons, and (3) the direction of transcriptional remodeling in VTA DA and glutamate neurons is similar in susceptibility and resilience. These findings support a model in which the VTA contributes to resilience through active, cell type-specific molecular adaptations in response to stress which serve to maintain behavioral stability in the face of adversity. Results CSDS produces susceptible and resilient groups To assess resilience versus susceptibility to stress, we subjected mice to CSDS. Mice underwent ten consecutive days of defeat from novel aggressors ( Figure 1B ) followed by 3-8 days of behavioral testing (Figure S1A). The social interaction test ( Figure 1C , S1B-E) was used to determine stress susceptibility and mice showed a spectrum of social interaction times with a novel aggressor following CSDS. In concordance with prior literature 30 , 32 , stressed mice were defined as susceptible if the time they spent in the interaction zone was less than one standard deviation below that of unstressed controls of the same sex; otherwise, they were considered resilient ( Figure 1C , S1C). Sex-specific thresholds were used because female mice were overall more social than males (two-way ANOVA with sex and stress as factors; main effect of sex, F=5.32, P=0.02). Unstressed male and female mice spent 60.23±8.3% (mean±s.d.) and 64.29±6.12% of time in the interaction zone near the novel aggressor, respectively. Using this measure, similar proportions of male and female mice were assigned to each group (Figure S1E). Susceptibility by this measure corresponded with anxiety-like behavior in non-social settings, such as the elevated plus maze (EPM; Figure 1D-E , S1J-M, S1O-P) and open field test (OFT; Figure S1F-I, S1N). However, it should be noted that behavior in these additional tests did not perfectly differentiate resilient from susceptible mice, likely due to behavioral complexity and individual variation in rodent models of stress. Transcriptional atlases of the mouse VTA and LHb using single-nucleus RNA-sequencing To examine how stress affects both the VTA and LHb, we generated snRNA-seq atlases of transcriptional diversity in both regions from ∼10 week old male and female C57BL/6J mice. Tissue was collected from control, susceptible, and resilient mice 10-15 minutes after a final social defeat re-exposure (or conspecific exposure for controls; Figure 1F ). We then isolated nuclei to generate snRNA-seq libraries and performed sequencing. Following quality control and preprocessing, we retained a total of 162,413 VTA nuclei across 16 samples and 228,187 habenula nuclei across 18 samples. VTA nuclei clustered into 3 neuronal and 6 non-neuronal populations ( Figure 2A , S2A-I; median unique molecular identifiers [UMIs]/nucleus=1,223; median genes/nucleus=821; Figure S2J). These populations, including well characterized glial cells (Figure S2), were recovered in proportions similar to those from previous literature 39 ( Figure 2B ) and were similarly represented across behavioral groups (Figure S3A-C). Download figure Open in new tab Figure 2. Across VTA cell types, resilient mice have greater gene expression changes and maintain receptor-ligand interaction strength. A. UMAP of all VTA nuclei with points colored by cell type classification. B. Proportion of VTA nuclei belonging to each cluster (n=162,413 total nuclei). C. Number of DEGs compared to control in susceptible (left; mixed neuron=5, DA neuron=0, GABA neuron=0, astro=17, micro=0, endo=0, OPC=0, premyelinating oligo=9, myelinating oligo=9) and resilient (right; mixed neuron=73, DA neuron=53, GABA neuron=0, astro=29, micro=0, endo=12, OPC=7, premyelinating oligo=15, myelinating oligo=18). D. DEGs counts in susceptible versus resilient groups compared to control, by cell type (Wilcoxon rank sum test with continuity correction, W=18, P=0.045; ∗p < 0.05). See Table S1 for statistics. E. Venn diagram of susceptible and resilient DEGs. F. Top GO terms associated with genes upregulated in resilience. G. Overall receptor-ligand interaction strength between all cell types. H-J. Circle plots showing differential receptor-ligand interaction strength between all cell types each pair of behavioral outcomes. Resilient nuclei maintain interaction strength similar to that of control (H; similar overall gray and gold edge width); susceptible nuclei have reduced interaction strength compared to control (I; greater overall gray than blue edge width); resilient nuclei have higher interaction strength than susceptible nuclei (J; greater overall gold than blue edge width). Consistent with prior reports 86 , our dataset showed that the overall cell type composition of the VTA is not significantly different between the sexes (Figure S4A-C). However, the VTA shows sex-dependent variation in response to stress 55 , 87 and drugs of abuse 88 – 91 , and receives input from brain regions that show substantial sex differences 92 . Therefore, we asked whether VTA cell types show sex-specific transcriptional signatures. Within each cluster, we pseudobulked data and performed pairwise comparisons of gene expression between males and females to reveal differentially expressed genes (DEGs) between the sexes (Figure S4D). The most dramatic sex differences were in known sex-linked genes (e.g. Tsix and Uty, Figure S4E-F). We also identified autosomal genes with strong sex bias across all cell types (e.g. Fkbp5 , which plays a key role in stress-related disorders 93 and has sex-dependent effects on stress susceptibility 94 ; Figure S4G) as well as those that exhibited sex-specific expression in certain cell types (e.g. Ube3a , Figure S4H). These differences may be associated with sex-dependent responses of the VTA in different physiological or pathological contexts. Habenula nuclei also clustered into distinct neuronal and non-neuronal populations (Figure S5A, left). We used known marker genes 35 , 36 , 43 to distinguish neurons from medial and lateral Hb (MHb and LHb, respectively) that were distinct from extraneously captured neurons from surrounding regions (Figure S5A). Since we could not be confident about which subregion non-neuronal cell types came from due to a lack of robust markers, we excluded them from further analysis. To examine neuronal heterogeneity within both MHb and LHb and ensure consistency with prior datasets 35 , 36 , 43 , we subclustered Hb neurons (Figure S5A, right). Subclustering revealed 7 transcriptionally distinct populations (4 LHb, 2 MHb, and one that could not be confidently assigned to either subregion) consisting of 108,270 total nuclei and 23,318 genes (UMIs/nucleus=2732; median genes/nucleus=1557; Figure S5F). All subclusters expressed high levels of Slc17a6 (Vglut2; Figure S5B,E) and low levels of GABAergic markers (Figure S5E), with LHb and MHb subclusters marked by expression of canonical markers Pcdh10 and Tac2 , respectively (Figure S5C-E). We also recapitulated discrete MHb subtypes expressing Chat or Calb2 (Figure S5E). Subpopulations were equally represented across sexes (Figure S5G) and behavioral groups (Figure S5H). LHb does not show strong patterns of gene expression change with stress Given significant research implicating hyperactivity of the LHb in aversive processing 14 – 18 , depression 17 , 20 – 22 , 61 , 64 , 95 , and stress susceptibility 23 , 30 , 96 , 97 , we asked whether stress there were significant gene expression changes in this region following stress. One study found divergent stress-associated gene expression in LHb, which hints at Hypothesis 3 23 . Our high resolution transcriptional atlas allowed us to also probe changes in MHb in addition to LHb neurons. Within LHb and MHb neurons separately, we performed pairwise comparisons of pseudobulked gene expression between control versus susceptible, control versus resilient, and resilient versus susceptible, including sex as a covariate. Unexpectedly, both regions showed few gene expression differences between conditions (adjusted pvalue<0.1; Figure S5I-N). We then calculated DEGs between LHb and MHb control samples and detected 5164 DEGs, indicating that our samples were sufficiently powered to detect large differences in habenula. Despite a lack of DEGs, LHb neurons of susceptible mice tended to be more active than those of resilient mice in response to a final aggressor exposure (measured by IEG expression; Figure S5O). Resilient mice show greater gene expression changes than susceptible across all VTA cell types while maintaining receptor-ligand interaction strength We then wondered if the VTA, a region associated with reward rather than aversion processing, may show greater transcriptional changes following stress. We first focused on both neuronal and non-neuronal clusters in order to understand how stress affects the entire VTA network. Within each cluster, we performed pairwise comparisons of pseudobulked gene expression between control versus susceptible and control versus resilient, including sex as a covariate. Resilience was consistently associated with a greater number of DEGs in all cell types where DEGs were detected ( Figure 2C-E ). To ensure that the higher number of DEGs in resilience was not due to one extra resilient sample of each sex (6 resilient vs 4 susceptible), we systematically removed every combination of one male and one female resilient sample and re-calculated DEGs. Across these held out combinations, the number of resilient DEGs remained significantly greater than susceptible DEGs (Figure S6). Resilient DEGs in non-neuronal cell types (Figure S7A-F) included several genes with putative roles in supporting network stability and function. For example, Grm5 (metabotropic glutamate receptor 5), which plays a role in modulating synaptic plasticity and glutamate homeostasis, was upregulated in astrocytes (Figure S7A,G). This finding is consistent with work implicating astrocytic glutamate regulation in stress response and psychiatric disease 98 – 100 . Ptgfrn is involved in cell adhesion and was upregulated in OPCs (Figure S7B,H). The glucocorticoid-responsive kinase Sgk1 , which was upregulated in myelinating oligodendrocytes (Figure S7D,I), has been implicated in processes important for maintaining circuit integrity such as myelin plasticity and stress hormone signaling 101 – 105 . Bmp6 was altered in endothelial nuclei (Figure S7F,J) and is a regulator of blood-brain barrier function, potentially suggesting vascular adaptations in resilience that promote behavioral stability. We performed gene ontology (GO) analysis on resilient DEGs to more systematically investigate the biological processes most altered in resilience. Across all clusters, genes upregulated in resilient animals were primarily associated with processes related to neural signaling and intercellular communication ( Figure 2F ), including chemical synaptic transmission (GO term: 0007268), cell junction assembly (GO term: 0034329), and cell adhesion (GO term: 0007155). The coordinated upregulation of genes involved in the molecular interfaces between cells and intercellular signaling in resilience may reflect the promotion of stable and efficient neural communication in order to encourage an adaptive stress response. In order to further probe cellular communication, we used CellChat 106 , 107 to compare inferred cell-cell interactions across the VTA network between each behavioral group (control, susceptible, resilient). CellChat revealed a decrease in overall interaction strength in susceptibility compared to control ( Figure 2G,I ), while resilience maintained interaction strength similar to that of control ( Figure 2G-H ). Comparing resilient and susceptibility directly highlighted stronger interactions in resilience ( Figure 2J ), supporting the hypothesis that the gene expression changes seen in resilience allow the VTA to maintain this interaction strength. Resilience is associated with greater gene expression changes in VTA neurons We found the greatest number of DEGs and the strongest cell-cell interaction strength within neuronal clusters and, given the diversity of neuronal subtypes in the VTA, we wanted to probe these effects with finer resolution. Subclustering of VTA neurons revealed a large dopaminergic (DA) cluster as well as glutamatergic, mixed, and multiple GABAergic populations ( Figure 3A-E ). Download figure Open in new tab Figure 3. Resilience is associated with greater transcriptional changes in VTA DA and glutamate neurons, while Sst-expressing GABA neurons show differences between susceptibility and resilience. A. UMAP of all VTA nuclei with neurons highlighted (left) and neuronal subclusters (right) with points colored by cell type classification. B. Feature plot showing expression of Th across VTA neuron UMAP space. C. Same as B, but for Gad2. D. Same as B-C, but for Slc17a6 (Vglut2). E. Number of nuclei per VTA neuronal subcluster (top) and expression of select marker genes across VTA neuronal subclusters (bottom). F. Number of DEGs compared to control in susceptible (left; DA=0, Mixed=0, Glut=2, GABA1=3, GABA2=1, GABA3=1, GABA4=0, GABA5=0, GABA6=0) and resilient (right; DA=61, Mixed=11, Glut=177, GABA1=4, GABA2=1, GABA3=1, GABA4=21, GABA5=3, GABA6=10) neurons. G. Summary of DEGs counts in susceptible versus control and resilient versus control, by cell type (Wilcoxon rank sum test with continuity correction, W=6.5, P=0.003; ∗∗p < 0.01). See Table S1 for statistics. H. Venn diagram of susceptible and resilient DEGs. I. Volcano plots showing susceptible DEGs (left) and resilient DEGs (right) in DA neurons. Nominal p-values are plotted, but points are colored based on adjusted p-values. Adjusted p-value (Bonferroni post hoc correction) cutoff is 0.1. J. Expression of the DA DEG Grik2 across behavioral groups. Top: each point represents the normalized expression of Grik2 across all DA neurons in each sample. Differential expression was assessed using DESeq2 on pseudobulked counts (log2FC=0.38, Wald statistic=4.07, p.adjusted=0.03). Bottom: each point represents the proportion of DA nuclei in each sample that express Grik2. K. Same as I, but for glutamate neurons. L. Same as J, but for the glutamate DEG Kcna1 (log2FC=0.94, Wald statistic=3.64, p.adjusted=0.04). M. Top GO terms associated with genes altered in DA neurons of resilient mice compared to control. N. Top GO terms associated with genes altered in glutamate neurons of resilient mice compared to control. O. Number of DEGs between susceptible and resilient (DA=8, Mixed=1, Glut=0, GABA1=1, GABA2=1, GABA3=0, GABA4=0, GABA5=5, GABA6=65) neurons. P. Volcano plot showing DEGs between susceptible and resilient in GABA6 neurons. Nominal p-values are plotted, but points are colored based on adjusted p-values. Adjusted p-value (Bonferroni post hoc correction) cutoff is 0.1. Top GO terms associated with genes differentially expressed in GABA6 neurons of resilient versus susceptible mice. The six GABAergic clusters, marked by high levels of Gad2 and Slc32a1 ( Figure 3C,E , S8A-H), reflected the transcriptomic diversity of VTA GABA neurons. These neurons modulate local DA activity 65 , 66 , 108 , 109 , send long-range projections 75 , 110 – 116 , play diverse roles in processing reward and aversion 53 , 65 , 66 , 70 – 72 , 74 , 76 , and exhibit transcriptomic heterogeneity 117 – 120 . Two clusters (GABA3 and GABA4) expressed high levels of Htr2c (Figure S8B,H) and likely represent populations previously identified in the VTA 39 , 121 . GABA5 and GABA6 expressed canonical interneuron markers such as Sst , Kit , Nos1 , Oprm1 , Vip, and Sox6 122 (Figure S8F-H). While these genes typically mark non-projecting GABA neurons in other brain regions, there is still debate surrounding whether canonical interneurons exist within the VTA and whether GABAergic VTA projection neurons also have local collaterals that can modulate the activity of neighboring neurons 123 . To further characterize these neuronal subtypes, we mapped our neuronal nuclei to the Allen Brain Cell (ABC) Atlas using the Allen Institute’s mapmycells function (see Methods). The mapping results are summarized in Table S6, which lists the top ABC cell classes accounting for up to 90% of the mapping for each of our assigned clusters. The majority of our clusters corresponded to midbrain neurons of the same neurotransmitter type, consistent with accurate anatomical and molecular annotation. Clusters GABA3 and GABA5 showed predominant mapping to pons GABA neurons, which may reflect minor dissection variability or shared transcriptional signatures among inhibitory neurons across neighboring regions. To assess cell type-specific transcriptional changes related to stress response, we performed pairwise pseudobulked comparisons of gene expression between control versus susceptible, control versus resilient, and susceptible versus resilient within each VTA neuronal subcluster. Compared to control, resilience was associated with greater or equal numbers of DEGs than susceptibility in all subclusters ( Figure 3F-H ). VTA DA and glutamate neurons had the greatest number of DEGs ( Figure 3F ). DA DEGs included potassium channels such as Kcnab2 , cell adhesion protein Cntnap2 that is known to cluster around potassium channels and play a crucial role in synaptic connections, amino acid transporter Slc1a4 , and both acetylcholine and glutamate receptors ( Chrm3 , Grik2 , Gria4 ; Figure 3I-J ). Glutamatergic DEGs included potassium channels ( Kcna1 , Kcng2 ), Pak6 (a kinase implicated in synaptic plasticity and neurotransmission), corticotropin releasing hormone receptor Crhr1 , and thyroid hormone transporter Slc16a2 ( Figure 3K-L ). Given that VTADA neurons are more heterogeneous than previously thought, we also looked for DA subcluster-specific DEGs. Subclustering of VTADA neurons (Figure S9A) revealed similar subpopulations to those found in previously collected sc- and snRNA-seq datasets generated from mouse and rat midbrain 37 – 39 , 124 . We observed ubiquitous expression of canonical DA neuron markers (e.g. Th, Slc6a3, Ddc, Slc18a2 ) as well as additional VTADA subpopulations previously identified to be anatomically and functionally-distinct (Figure S9B-H). While we saw the same pattern of greater change in resilience than susceptibility, we surprisingly found fewer DEGs overall among VTADA subclusters than in the higher level DA cluster (Figure S9I). Subcluster DA-2 had the highest number of DEGs (Figure S9I) and may be driving the resilience-associated changes in DA neurons, but the overall low number of DEGs at the subcluster level suggests that transcriptional changes associated with resilience may be distributed across dopaminergic subtypes rather than concentrated within specific ones. Due to the low number of DA subcluster DEGs, we conducted downstream analyses at the higher DA cluster level, but future studies may probe the role of individual DA subclusters in responding to stress. We performed GO analysis on resilient DEGs in DA and glutamate neurons ( Figure 3M-N ). In resilient animals, DA neuron DEGs were associated with synaptic transmission (GO terms: 0007268, 0098815, 0097106), neuronal excitability (GO terms: 0019228, 0099610), and structural connectivity (GO terms: 0008090, 0000902, 0007157), suggesting adaptations that may fine-tune how these neurons communicate in real-time. In contrast, glutamate neuron DEGs were associated with developmental and structural plasticity processes such as neuron migration (GO term: 0001764) and projection growth (GO terms: 0031175, 0048812) as well as MAP kinase signaling (GO term: 0043410). These terms suggest long-term remodeling that might support behavioral flexibility and learning associated with resilience. The one population that exhibited a large number of DEGs between susceptible and resilient animals was a specific population of GABAergic VTA neurons (GABA6; Figure 3O-P ). This cluster was characterized by high expression of Sst, Oprm1, Kit, Nos1, and other peptidergic and synaptic regulatory genes (Figure S8G-H). According to our ABC mapping, GABA6 most closely corresponded to ABC midbrain (MB) GABA neurons (Table S6). This mapping reinforces that the observed transcriptional differences between resilient and susceptible animals are intrinsic to a midbrain GABAergic population rather than driven by contamination from neighboring regions. GO analysis of susceptible versus resilient GABA6 DEGs highlighted differences in neuronal development, synapse assembly, and signaling pathways, suggesting altered circuit connectivity and plasticity between the two phenotypes in this population ( Figure 3Q ). Linear discriminant analysis reveals greater molecular divergence in resilience than susceptibility compared to control In VTA DA and glutamate neurons, we observed more DEGs between resilient and control than between susceptible and control. These findings challenge the hypothesis that susceptibility represents a greater molecular shift from the control state than resilience ( Figure 1A , Hypothesis 1). In addition, the small number of DEGs between susceptible and resilient in these neuronal clusters challenge the idea that susceptibility and resilience diverge in distinct directions away from control ( Figure 1A , Hypothesis 3). Instead, our DEG analysis supports the idea that resilience represents a greater shift in expression from control in these cell types ( Figure 1A , Hypothesis 2). Alternatively, one population of GABAergic VTA neurons characterized by high expression of peptidergic and synaptic regulatory genes (GABA6) had few DEGs in either resilient or susceptible compared to control but showed large differences between resilient and susceptible directly ( Figure 3O ). This pattern suggests that susceptibility and resilience diverge more from each other than from control ( Figure 1A , Hypothesis 3). To investigate the overall transcriptional structure between behavioral groups and test these hypotheses more directly, we applied linear discriminant analysis (LDA) to each of these three clusters (DA, glutamate, GABA6). This approach allowed us to identify the axes that best distinguish all three behavioral groups ( Figure 4A ). In both DA and glutamate neurons, group centroids laid approximately along a linear trajectory ( Figure 4B, D ), exhibiting a graded shift in transcriptional state from control to susceptible to resilient that supports Hypothesis 2. In DA neurons in particular, when nuclei from susceptible mice were misclassified, they were more often labeled as control than resilient ( Figure 4C ), reinforcing their transcriptional proximity to control. This also occurred to a lesser extent among glutamate neurons ( Figure 4E ). These results suggest that a one-dimensional axis could effectively capture the main group differences in these subclusters and that, in these cell types, gene expression in susceptible mice is most similar to control. Download figure Open in new tab Figure 4. LDA identifies cell type-specific transcriptional trajectories in resilient and susceptible VTA neurons. A. Schematic of the cross-validation (CV) workflow for 2D LDA. B. 2D LDA projections of held-out nuclei from DA neurons, colored by true group labels. C. Confusion matrix showing classification accuracy of held-out DA nuclei, aggregated across all folds. D. Same as B, but for glutamate neurons. E. Same as C, but for glutamate neurons. F. Same as B, but for GABA6 neurons. G. Same as C, but for GABA6 neurons. H. 1D LDA projection of DA nuclei onto LD1 trained only and control and resilient. I. Pie chart summarizing predicted classifications for susceptible nuclei projected on the control-resilient LD1 axis. J. Same as F, but for glutamate nuclei K. Same as G, but for glutamate nuclei. LDA supports hypothesis 2 for how resilience and susceptibility are instantiated in VTA dopamine and glutamate neurons and hypothesis 3 for a specific population of VTA GABA neurons (GABA6). To further test this idea, we performed one-dimensional LDA on DA and glutamate neurons using only control and resilient nuclei to derive the linear discriminant (LD1) that best separates these two groups. We then projected susceptible nuclei onto this axis to examine their relative positioning ( Figure 4H, J ). In both DA and glutamate neurons, susceptible nuclei indeed landed closer to control than to resilient along LD1, and were more often classified as control than resilient ( Figure 4H-K ). These results reinforce the idea that in these cell types, susceptibility maintains overall transcriptional structure much closer to control than resilience does ( Figure 4L , Hypothesis 2). In contrast, in GABA6 neurons all three group centroids occupied distinct regions in LDA space ( Figure 4F ) and were classified well above chance ( Figure 4G ). These results support the idea that resilience and susceptibility are instantiated by discrete gene expression structures away from control within this population ( Figure 4L , Hypothesis 3). Together, these analyses uncover cell type-specific transcriptional structures in VTA that distinguish susceptible and resilient phenotypes following chronic social stress. Discussion Chronic stress is a significant trigger for depression and other mood disorders, and therefore the study of mechanisms that confer resilience to stress are of great importance for public health. Here, we demonstrated that resilience to chronic stress involves active, cell type-specific transcriptional changes in the VTA and is not simply a passive lack of susceptibility. Surprisingly, despite its well-established role in aversive signaling, the LHb exhibited minimal gene expression changes. Limited transcriptional changes in the LHb suggest alternative mechanisms of plasticity While the LHb is widely implicated in aversive processing and stress-related behavior 14 – 30 , we observed few transcriptional changes in this region in either susceptible or resilient animals following CSDS. These unexpected findings suggest that the LHb may exert its influence through alternative mechanisms that do not involve sustained transcriptomic remodeling within neurons. One possibility is that stress induces transient gene expression changes that drive synaptic remodeling, which can maintain changes in neuronal activity with normal gene expression. A similar pattern of rapid transcriptional response to acute stress followed by long-term dampening of this response as stress becomes chronic has recently been reported in ventral hippocampus 125 . This underscores the importance of assessing multiple time points and complementing transcriptomic studies with circuit-level and physiological approaches to fully capture regional dynamics. In addition, studies have shown that stress-induced alterations in LHb astrocyte function can modulate LHb neuron firing and depressive-like behavior 57 – 59 , 126 . These effects may be reflected as astrocyte-specific transcriptional changes and warrant further investigation in future studies. Cell type specificity of stress adaptations in the VTA Within VTA neurons, we observed distinct gene expression patterns across dopaminergic, glutamatergic, and GABAergic subpopulations. These results reveal that susceptibility and resilience diverge in their transcriptional trajectories across neuron subtypes and highlight that individual cell types within the VTA adapt to stress heterogeneously. Specifically, within VTA DA and glutamate neurons, our analyses support the idea that resilience recruits more pronounced adaptations in a similar direction as susceptibility away from control ( Figure 4L , Hypothesis 2). On the other hand, in a specific population of Sst- expressing GABA neurons (GABA6), our data supports the hypothesis that susceptibility and resilience represent fundamentally different gene expression programs, moving away from control in divergent directions ( Figure 4L , Hypothesis 3). Together, this reinforces the importance of examining unique cell types in this region to understand the precise molecular mechanisms that differentiate susceptibility and resilience to chronic stress. Resilience as an active adaptive response in the VTA The VTA showed widespread transcriptional changes across both neuronal and non-neuronal cell types, with resilience consistently associated with more DEGs than susceptibility, supporting and extending earlier bulk RNA-seq findings 10 . These changes were in genes related to synaptic transmission, structural plasticity, and cell-cell communication, suggesting that resilient animals engage regulatory programs to stabilize neural function and preserve network connectivity. These changes were largest in VTA DA and glutamate neurons. Notably, susceptibility involved far fewer transcriptional changes despite producing strong behavioral phenotypes, suggesting that resilience requires a more coordinated molecular response. Further, multivariate LDA analysis of overall transcriptional state found that within both VTA DA and glutamate neurons, group centroids laid approximately along a linear trajectory from control to susceptible to resilient. This graded positioning reinforces the asymmetry we found in the number of DEGs: greater transcriptional distance between control and resilient nuclei than control and susceptible nuclei reflects more extensive transcriptional remodeling in resilience, whereas susceptibility maintains overall transcriptional structure closer to control. In contrast, the closer proximity of susceptible nuclei to control is represented by a more subtle shift in gene expression, represented by fewer DEGs. These findings support the hypothesis that within VTA DA and glutamate neurons resilience is not merely the absence of susceptibility, but is instead an active transcriptional state that extends further from control than susceptibility does and likely engages stronger regulatory programs. This view aligns with electrophysiological findings that resilient animals maintain stable VTADA neuron firing, likely through compensatory, homeostatic mechanisms that preserve neural function under stress 33 . These findings challenge the prevailing view that to ameliorate the changes in mood, behavior, and neural activity associated with stress susceptibility, we must reverse the molecular consequences of stress. By identifying resilience as an active and widespread transcriptional program in the VTA, our work provides a foundation for future efforts to target cell type-specific pathways that support adaptive responses to chronic stress. Whereas susceptibility may arise from a failure to activate sufficient protective mechanisms, resilience appears to involve the engagement of gene networks that promote neural stability and flexible behavioral coping. Similar forms of structural plasticity have been observed in other brain regions such as the hippocampus, where resilience is associated with compensatory postsynaptic remodeling in CA1 127 . Identifying upstream drivers of these cell type-specific transcriptional programs could inform new strategies to enhance natural resilience mechanisms. Such approaches may offer a compelling alternative to current treatments, which largely aim to suppress symptoms of mood disorders rather than promote resilience. This concept is supported by prior studies in other brain regions such as the nucleus accumbens, where similar efforts to harness endogenous resilience mechanisms have shown potential as therapeutic strategies for depression 85 , 128 , 129 . Caveats and Considerations Although we performed snRNA-seq on 16 VTA samples from 32 mice and 18 un-pooled Hb samples, we are likely underpowered to detect small changes in expression given the need to correct for multiple comparisons across the genome 130 . However, we performed several analyses to ensure that our findings are not simply due to small sample sizes. Iterative DEseq2 analysis of VTA data, in which we dropped every combination of male and female resilient samples, confirmed that an imbalance in sample size does not account for a greater number of resilient DEGs (Figure S6). The conclusion that resilience involves greater sustained transcriptional change compared to susceptibility is further supported by multivariate LDA modeling. Although we did not identify stress-associated DEGs in the Hb, comparison of gene expression between LHb and MHb control samples identified thousands of DEGs, indicating the current sample sizes are sufficient to detect the most robust differences. Nevertheless, given sample size considerations, it would be naïve to conclude that stress does not result in any gene expression differences in the VTA of susceptible mice or in the Hb of either group; rather, we can only conclude that the most robust changes are detectable among resilient mice. Although both males and females were used and sex was specifically considered in our study design and analysis, additional considerations remain. To examine both sexes, we used chronic nondiscriminatory social defeat stress (CNSDS) 11, 131 , a modified version of CSDS in which a male and a female mouse are placed together into the homecage of a novel aggressor to induce nondiscriminatory attacks, including towards the female. Although male mice tend to be attacked more than females in this paradigm, and therefore the severity of the stress may not be the same between the sexes, this approach nonetheless allowed us to use the same stressor to generate resilient and susceptible phenotypes from both sexes, which was key to the goals of this study. We controlled for the effects of sex by using sex as a factor in all DEseq2 models, but this study was not powered to examine sex-specific stress effects. Given prior evidence of sex-dependent variation in response to stress in the VTA 55, 87 Here, transcriptional profiling was performed at a single time point following CSDS, limiting the ability to capture dynamics or transient gene expression changes over time. In particular, additional time points immediately following the first stress experience, and later time points after recovery, may resolve additional signatures of susceptibility and resilience, respectively. While comprehensive longitudinal snRNA-seq remains prohibitively expensive at the scale necessary to fully capture these dynamics, reduced time-course designs or computational trajectory inference could better resolve the temporal progression of stress-induced molecular adaptations. Finally, the analyses here have focused on characterization of the molecular signatures of resilience and susceptibility after chronic adult social stress, and have not yet functionally established the causal role of any of the specific gene expression changes observed. Future studies utilizing gene perturbation or functional assays will be necessary to test whether identified DEGs causally influence resilience. Contributions A.R.M, I.B.W., and C.J.P. conceived of the project and designed the experiments. A.R.M., B.M., R.L., and A.Z. collected the data. A.R.M. analyzed the data. A.R.M., I.B.W., and C.J.P. wrote the manuscript. Declaration of Interests The authors declare no competing interests. Data and Code Availability Single-nucleus RNA-sequencing data have been deposited at GEO and are awaiting a GEO accession number. Code for analyzing behavior and snRNA-seq data is available on GitHub: https://github.com/aminerva/CSDS-SingleCell-VTA-Hb . Any additional information required to reanalyze the data reported in this paper is available from the Lead Contact upon request. Supplemental Figures Figure S1: Post hoc behavioral tests correlate with susceptible and resilient phenotypes from the social interaction test Figure S2. Characterization and quality control metrics from snRNA-sequencing of VTA Figure S3. Quality control metrics from snRNA-sequencing of VTA across behavioral groups Figure S4. Sex differences across VTA cell types Figure S5. Gene expression in the Hb does not uniquely define resilience or susceptibility Figure S6. Differences in resilient vs control and susceptible vs control DEG counts is not due to sample size differences Figure S7. Example DEGs in non-neuronal cell types of the VTA Figure S8. GABAergic neuron diversity in the VTA Figure S9. DA neuron diversity in the VTA Supplemental Tables Table S1: Stats for main figures Table S2: Stats for supp figures Table S3: DESeq2 results for all VTA cells Table S4: DESeq2 results for VTA neurons Table S5: DESeq2 results for Hb Table S6: Results from mapping neuronal clusters to the Allen Institute’s brain atlas cell classes Methods Mice All experiments were approved by the Princeton University Institutional Animal Care and Use Committee and were in accordance with National Institutes of Health standards. Prior to and throughout experimental assays, animals were housed under a 12H light-dark cycle. Food and water were given ad libitum . Experimental animals used in this study were comprised of 72 mice (36 males and 36 females) on a C57/BL6 WT background (Jackson Laboratory) between the ages of 10 and 16 weeks. Tissue from 32 of these mice was used for sequencing. Aggressor mice were Swiss Webster male retired breeders (Taconic). Chronic social defeat stress (CSDS) Both male and female mice were exposed to chronic non-discriminatory social defeat stress (CNSDS) 11 , 131 , a modified version of chronic social defeat stress 9 . Male mice were placed in the shoebox home cage (#5 Expanded Mouse Cage 22.2cm x 30.8cm x 16.2cm, Thoren Caging Systems, Inc.) of a novel aggressor mouse for 3 minutes of free interaction, followed by an additional 5 minutes during which a female was added. While aggressor mice do not reliably attack females on their own, the inclusion of both a male and a female induces nondiscriminatory attacks from the aggressor. Males were separated from the aggressor by a perforated acrylic barrier (Tap Plastics) within the aggressor cage, while females were housed in a new cage across an acrylic barrier from a different aggressor with a handful of soiled bedding from the aggressor they were defeated by that day. These housing conditions remained overnight. 24H later, mice were placed in the cage of a new aggressor and males and females were rotated in different directions so that trios were unique each day. This continued for a total of 10 days. Unstressed control mice were pair-housed with a perforated barrier separating the two mice. They were handled and their cages rotated each day for 10 days. Following the defeat on day 10, aggressors were removed and all mice were singly housed through the remaining stages of the experiment. Social interaction test To assess social avoidance behavior following the final day of CSDS, experimental mice were tested in a two-stage social interaction test during their light cycle from ∼9am-12pm. In the first 2.5 min stage, experimental mice were allowed to freely explore an arena (44 × 44 × 20 cm) containing a novel plexiglass and wire mesh enclosure (10 × 6 cm) centered against one wall of the arena. The mouse was then removed from the chamber, and a novel Swiss Webster aggressor mouse was placed within the plexiglass and mesh enclosure. The experimental mouse was then returned to the arena for an additional 2.5 min. Mouse location was tracked via Ethovision (Noldus) and used to measure time spent in the “interaction zone” (14 × 26 cm) surrounding the enclosure, “corner zones” (10 × 10 cm), and “distance traveled”. Time spent in the social interaction zone when the social target was a novel Swiss Webster aggressor was used to delineate resilient and susceptible mice. We defined susceptible as 1 standard deviation below the mean social interaction time of the unstressed control group as done in previous literature 30 , 32 . Elevated plus maze Mice were tested in the elevated plus maze (EPM) during their light cycle from ∼9am-12pm (85-101 lux in open arms). Each mouse was placed in the center of the EPM (2 enclosed arms and 2 open arms; each arm 76cm long and 6.5cm wide). The mouse explored the maze for 7 minutes, while its centroid location was tracked via Ethovision (Noldus) and used to measure the time spent in the open arms and center of the maze. Open field test Mice were tested in the open field test (OFT) during their light cycle from ∼9am-12pm (85-101 lux). Each mouse was placed into the center of an empty area (50cm x 50cm). Lamps were used to illuminate the arena on the left and right so that there was a shaded area along the left and right walls. The center was 42cm x 42cm centered at the center of the arena. The animal explored for 10 minutes, while its centroid location was tracked via Ethovision (Noldus) and used to measure the time spent in the center of the arena. Tissue collection and single-nuclei dissociation for sequencing We collected VTA and habenula tissue from male and female mice in adulthood (>P60). Following 10 days of CSDS, mice underwent 3-8 days of behavioral testing (described above). On day 19, mice were re-exposed to a novel aggressor (or novel conspecific for control mice) for 5 minutes. ∼10-15 minutes after re-exposure, animals from all conditions were cervically dislocated. This time delay is in line with prior reports indicating sufficient IEG exposure at this time following a social exposure 132 . Brains were extracted and sectioned into 1mm coronal sections on ice using a brain block. Bilateral punches of VTA and one punch encompassing both hemispheres of habenula were taken, flash frozen in Eppendorf tubes, and kept at -80C until nuclei isolation and library preparation. Nuclei were isolated into suspensions for snRNA-seq as previously described, 133 with minor modifications. For VTA, tissue punches from two animals of the same sex and condition were pooled together to ensure a high enough concentration of nuclei per sample for quality sequencing and were then Dounce homogenized. Pairs were selected randomly from within each sex and condition. The sample was filtered through a 40um strainer and nuclei were isolated via iodixanol gradient centrifugation. For habenula, a subset of the animals pooled for VTA samples were used for sequencing and each sample contained tissue from only one animal. Nuclei isolation continued as above, except that nuclei were pelleted by centrifugation rather than isolated via gradient. Nuclei were hand-counted using a hemocytometer on an EVOS M5000 microscope and the volume of suspension used for library preparation was adjusted in order to load a target of 20,000 nuclei and capture 10,000 nuclei per sample. snRNA-seq library prep and sequencing Nuclei were captured and RNA was reverse transcribed and barcoded for library preparation using the 10X Genomics Chromium v3 (VTA) or v4 (habenula) platform according to the manufacturer’s instructions. cDNA was amplified and adapters were added for sequencing on an Illumina NovaSeq SP 100nt Lane v1.5. Library preparation and sequencing was performed by the Princeton University Genomics Core. snRNA-seq read mapping The 10X Genomics package CellRanger (v6.1.1) was used to map transcripts to the mm10 mouse reference genome. In order to capture the high percentage of pre-spliced intronic mRNA present within the nucleus, the ‘–include-introns’ flag was used to map unspliced reads to corresponding genes. QUANTIFICATION AND STATISTICAL ANALYSIS Statistical analysis of behaviors All statistical tests of behavior (SI test, EPM, OFT) were analyzed using custom Python scripts. To assess differences in each behavioral metric across groups, a one-way analysis of variance (ANOVA) was performed with behavioral metric as the dependent variable and behavior group as the independent factor. When a significant main effect was observed, post hoc comparisons were conducted using Tukey’s Honestly Significant Difference (HSD) test to determine pairwise group differences. These statistical analyses were performed using the statsmodels module. To assess correlations between behavioral metrics, a Pearson correlation analysis was performed between the two metrics. Statistical values were computed using the scipy.stats module. snRNA-sequencing analysis snRNA-seq read mapping The 10X Genomics package CellRanger (v6.1.1) was used to map transcripts to the mm10 mouse reference genome. In order to capture the high percentage of pre-spliced intronic mRNA present within the nucleus, the ‘--include-introns’ flag was used to map unspliced reads to corresponding genes. After quality control and initial filtering, we recovered a total of 228,187 habenula nuclei and 162,413 VTA nuclei. Preprocessing CellRanger output count matrices were further analyzed using Seurat (v4.0.0) 134 , 135 in R (v4.0.3). During microfluidic encapsulation of nuclei on 10X devices, some droplets may contain more than one nucleus. During reverse transcription, transcripts from co-encapsulated nuclei are labeled with the same cell barcode, creating a ‘doublet’. To remove likely doublets in each sample, nuclei with >2500 genes were removed. In addition, nuclei with 5% of reads mapping to mitochondrial genes were removed. Mitochondrial and ribosomal genes were then removed in order to remove these genes as a confounding source of variation in downstream analysis. Dataset integration, UMAP embedding, and clustering VTA and Hb datasets were analyzed separately. To integrate all data from the same brain region and assign cell type identity to nuclei, each sample dataset was log normalized and the top 2000 highly variable features were identified for downstream sample integration using the FindVariableFeatures() function. These features were used as input to the function SelectIntegrationFeatures() . We then scaled and ran PCA on each dataset separately. The 16 separate VTA datasets–6 control (3 male and 3 female), 6 resilient (3 male and 3 female), and 4 susceptible (2 male and 2 female)–and 18 separate habenula datasets–6 control (3 male and 3 female), 6 resilient (3 male and 3 female), and 6 susceptible (3 male and 3 female)–were then integrated using the FindIntegrationAnchors() and IntegrateData() functions. Next, we scaled and ran PCA on both integrated datasets. The RunUMAP() function was used to embed the nuclei into UMAP space and the FindNeighbors() and FindClusters() functions were used to cluster nuclei based on overall gene expression similarity. Known marker genes were used to assign cell type identities to each cluster. Quality control metrics In both datasets, nuclei from male and female samples were distributed evenly across clusters (Figure S4A-C, S5G). Nuclei were similarly well distributed across samples from different behavioral groups (control, susceptible, resilient; Figure S3A-C, S5H). Calculation of DEGs between sexes across all VTA cell types To identify DEGs between males and females in each VTA cluster, we performed pseudobulk DESeq2 analysis similar to above, but with slight modifications. The design formula included behavioral condition as a covariate to adjust for stress-specific effects: design = ∼ sex + condition. Normalized counts were calculated and genes were filtered based on expression thresholds to improve statistical power: genes with normalized counts above the 5th percentile in at least four samples were retained. The same FDR adjusted p-value threshold was used to determine significance. Identification of Hb neuronal subtypes We further subclustered neurons in the Hb dataset to identify neuronal subtypes in both the LHb and MHb. Clusters were identified as neuronal based on expression of Syt1 and as belonging to LHb or MHb based on expression of Pcdh10/Gap43 or Tac2 , respectively. Neurons that did not express either of these markers were determined to be from surrounding areas outside of the Hb. LHb and MHb neuron clustered were filtered out for additional subclustering. The RunUMAP() function was used to embed neuronal nuclei into UMAP space and the FindNeighbors() and FindClusters() functions were used to cluster based on overall gene expression similarity. Known marker genes were used to assign cell type identities to each subcluster. Calculation of stress DEGs within Hb clusters To identify DEGs between control, susceptible, and resilient in LHb and MHb, we performed pseudobulk DESeq2 136 analysis. For the LHb and MHb cluster separately, raw RNA counts were aggregated by sample ID to create pseudobulk count matrices, representing summed gene counts per sample within each cell type. The known sex chromosome markers Xist and Tsix were removed to reduce noise. For each cell type and comparison (susceptible vs control, resilient vs control, susceptible vs resilient), we created DESeq2 datasets using the pseudobulk counts and metadata. The design formula included sex as a covariate to adjust for sex-specific effects: design = ∼ sex + condition . Condition factors were re-leveled such that the reference group was control for comparisons to control or resilient for susceptible vs resilient. Normalized counts were calculated and genes were filtered based on expression thresholds to improve statistical power: genes with normalized counts above the 5th percentile (or a minimum of 25 counts) in at least two samples were retained. DESeq2 Wald tests were performed to identify DEGs and an adjusted p-value threshold (FDR) of 0.1 was applied to determine significance. This approach was iteratively applied to both the LHb and MHb cluster and each behavioral comparison, yielding comprehensive DEG lists. Gene ontology (GO) analysis Gene Ontology enrichment analysis was performed using the Database for Annotation, Visualization and Integrated Discovery (DAVID) to identify significantly enriched biological processes among differentially expressed genes. DEG lists were submitted using official gene symbols and analyzed against the Mus musculus background. Enriched GO terms were considered significant at a Benjamini-adjusted p-value (FDR) < 0.05. CellChat for analysis of VTA receptor-ligand interactions We used CellChat 106 , 107 to compare inferred cell-cell interactions across the entire VTA network in each behavioral group (control, susceptible, resilient). For each group, a CellChat object was initialized using the assigned cell type identities. Overexpressed genes and ligand-receptor pairs were identified, followed by computation of cell-cell communication probabilities as previously described 106 , 107 . Interactions involving fewer than 10 nuclei were filtered out. The CellChat objects for each group were merged using the mergeCellChat() function for comparative visualization. Identification of VTA neuronal subtypes In the VTA dataset, clustering confirmed the presence of expected neuronal and non-neuronal populations, which we labeled based on expression of established marker genes 38 , 39 ( Figure 2A-B , S3A-I). These populations included neuronal and non-neuronal cell types. Clusters were identified as neuronal based on expression of Syt1 (Figure S3A-B) and filtered out for additional subclustering. The RunUMAP() function was used to embed the neuronal nuclei into UMAP space and the FindNeighbors() and FindClusters() functions were used to cluster neuronal nuclei based on overall gene expression similarity. Known marker genes were used to assign identities to each subcluster ( Figure 3A-E ). Mapping of VTA neuronal subtypes to Allen Brain Cell Atlas To assess correspondence between experimentally derived neuronal clusters and reference-defined brain cell types, we used the Allen Institute’s Allen Brain Cell (ABC) Atlas MapMyCells tool (RRID:SCR_024672). A raw count matrix was submitted to the mapmycells function, which performs transcriptome-to-reference matching using a cosine similarity–based classifier trained on the ABC taxonomy. The output included the top-ranked ABC cell classes and subclasses corresponding to each input cluster, along with their relative mapping proportions. Results are summarized in Table S6, where ABC cell classes accounting for up to 90% of the mapping for each cluster are listed. This analysis was used to ensure accurate regional and neurotransmitter-level annotation. Calculation of DEGs within VTA neuronal subtypes To identify DEGs within each VTA neuronal subcluster, we performed pseudobulk DESeq2 analysis as described above. LDA analysis of VTA neurons To assess transcriptional similarity across behavioral outcomes, we performed shrinkage discriminant analysis (SDA), a regularized form of linear discriminant analysis (LDA) suited for high-dimensional data. Analyses were conducted separately for VTA DA, glutamate, and GABA6 neuron populations using normalized gene expression matrices. All SDA analyses were performed using the top 2,000 most variable genes as input features. Nuclei were grouped by outcome (control, susceptible, or resilient). We used hold-one-sample-out cross-validation to evaluate classification performance while avoiding overfitting to individual animals. Each fold held out all nuclei from one sample, and the SDA model was trained on the remaining samples. We used the sda() function from the sda R package 137 , which implements shrinkage estimation of the within-class covariance matrix using the Ledoit-Wolf estimator. For 2D SDA, all three groups (control, susceptible, resilient) were included in training and the top two discriminant axes (LD1 and LD2) were computed to maximize separation between all three groups. Test set nuclei were then projected into this space and assigned predicted class labels based on the trained model. Classification performance was assessed across all folds by computing confusion matrices and overall accuracy. To examine the positioning of susceptible nuclei along the control-resilient axis, we conducted a separate 1D SDA analysis. The SDA model was trained only on control and resilient nuclei to derive a single discriminant axis (LD1). Susceptible nuclei were then projected onto this axis and their predicted labels were assigned using the trained SDA classifier. Acknowledgements We thank Annegret Falkner, Fenna Krienen, Chris Zimmerman, and Lindsay Willmore for feedback, as well as support from the rest of the Peña and Witten labs. We also thank Wei Wang and Jennifer Miller of the Princeton Genomics Core for conducting library preparation and sequencing. This research was funded by NSF GRFP DGE-2039656 (A.R.M.), Princeton University Elizabeth Procter Fellowship (A.R.M.), Princeton Neuroscience Institute Research Innovator Award (C.J.P.), New York Stem Cell Foundation (C.J.P. & I.B.W.), the Howard Hughes Medical Institute (I.B.W.), NIH DP1-MH136573 (I.B.W.). Funder Information Declared U.S. National Science Foundation, https://ror.org/021nxhr62 , DGE-2039656 Princeton University , Elizabeth Procter Fellowship , Princeton Neuroscience Institute Research Innovator Award New York Stem Cell Foundation , Robertson Neuroscience Award Howard Hughes Medical Institute National Institutes of Health, https://ror.org/01cwqze88 , DP1-MH136573 References 1. ↵ Kessler , R. C. et al. Lifetime prevalence and age-of-onset distributions of DSM-IV disorders in the National Comorbidity Survey Replication . Arch. Gen. Psychiatry 62 , 593 – 602 ( 2005 ). OpenUrl CrossRef PubMed Web of Science 2. ↵ Charney , D. S. & Manji , H. K . Life stress, genes, and depression: multiple pathways lead to increased risk and new opportunities for intervention . Sci. STKE 2004, re5 ( 2004 ). 3. ↵ Huhman , K. L . Social conflict models: can they inform us about human psychopathology? Horm. Behav . 50 , 640 – 646 ( 2006 ). OpenUrl CrossRef PubMed Web of Science 4. Breslau , N. & Davis, G. C. Chronic stress and major depression . Arch. Gen. Psychiatry 43 , 309 – 314 ( 1986 ). OpenUrl CrossRef PubMed Web of Science 5. Kendler , K. S. , Karkowski , L. M. & Prescott , C. A . Causal relationship between stressful life events and the onset of major depression . Am. J. Psychiatry 156 , 837 – 841 ( 1999 ). OpenUrl CrossRef PubMed Web of Science 6. McEwen , B. S. & Akil , H . Revisiting the stress concept: Implications for affective disorders . J. Neurosci . 40 , 12 – 21 ( 2020 ). OpenUrl Abstract / FREE Full Text 7. ↵ Daviu , N. , Bruchas , M. R. , Moghaddam , B. , Sandi , C. & Beyeler , A . Neurobiological links between stress and anxiety . Neurobiol. Stress 11 , 100191 ( 2019 ). 8. ↵ Hoge , E. A. , Austin , E. D. & Pollack , M. H . Resilience: research evidence and conceptual considerations for posttraumatic stress disorder . Depress. Anxiety 24 , 139 – 152 ( 2007 ). OpenUrl CrossRef PubMed Web of Science 9. ↵ Berton , O. et al. Essential role of BDNF in the mesolimbic dopamine pathway in social defeat stress . Science 311 , 864 – 868 ( 2006 ). OpenUrl Abstract / FREE Full Text 10. ↵ Krishnan , V. et al. Molecular adaptations underlying susceptibility and resistance to social defeat in brain reward regions . Cell 131 , 391 – 404 ( 2007 ). OpenUrl CrossRef PubMed Web of Science 11. ↵ Yohn , C. N. et al. Chronic non-discriminatory social defeat is an effective chronic stress paradigm for both male and female mice . Neuropsychopharmacology 44 , 2220 – 2229 ( 2019 ). OpenUrl CrossRef PubMed 12. Chaudhury , D. et al. Rapid regulation of depression-related behaviours by control of midbrain dopamine neurons . Nature 493 , 532 – 536 ( 2013 ). OpenUrl CrossRef PubMed Web of Science 13. ↵ Cao , J.-L. et al. Mesolimbic dopamine neurons in the brain reward circuit mediate susceptibility to social defeat and antidepressant action . J. Neurosci . 30 , 16453 – 16458 ( 2010 ). OpenUrl Abstract / FREE Full Text 14. ↵ Stamatakis , A. M. & Stuber , G. D . Activation of lateral habenula inputs to the ventral midbrain promotes behavioral avoidance . Nat. Neurosci . 15 , 1105 – 1107 ( 2012 ). OpenUrl CrossRef PubMed 15. ↵ Wang , D. et al. Learning shapes the aversion and reward responses of lateral habenula neurons . Elife 6 , e23045 ( 2017 ). OpenUrl CrossRef PubMed 16. Matsumoto , M. & Hikosaka , O . Lateral habenula as a source of negative reward signals in dopamine neurons . Nature 447 , 1111 – 1115 ( 2007 ). OpenUrl CrossRef PubMed Web of Science 17. ↵ Hu , H. , Cui , Y. & Yang , Y . Circuits and functions of the lateral habenula in health and in disease . Nat. Rev. Neurosci . 21 , 277 – 295 ( 2020 ). OpenUrl CrossRef PubMed 18. ↵ Matsumoto , M. & Hikosaka , O . Representation of negative motivational value in the primate lateral habenula . Nat. Neurosci . 12 , 77 – 84 ( 2009 ). OpenUrl CrossRef PubMed Web of Science 19. ↵ Cameron , S. , Weston-Green , K. & Newell , K. A . The disappointment centre of the brain gets exciting: a systematic review of habenula dysfunction in depression . Transl. Psychiatry 14 , 499 ( 2024 ). 20. ↵ Li , B. et al. Synaptic potentiation onto habenula neurons in the learned helplessness model of depression . Nature 470 , 535 – 539 ( 2011 ). OpenUrl CrossRef PubMed Web of Science 21. ↵ Li , K. et al. βCaMKII in lateral habenula mediates core symptoms of depression . Science 341 , 1016 – 1020 ( 2013 ). OpenUrl Abstract / FREE Full Text 22. ↵ Lecca , S. et al. Rescue of GABAB and GIRK function in the lateral habenula by protein phosphatase 2A inhibition ameliorates depression-like phenotypes in mice . Nat. Med . 22 , 254 – 261 ( 2016 ). OpenUrl CrossRef PubMed 23. ↵ Cerniauskas , I. et al. Chronic stress induces activity, synaptic, and transcriptional remodeling of the lateral habenula associated with deficits in motivated behaviors . Neuron 104 , 899 – 915 .e8 ( 2019 ). OpenUrl CrossRef PubMed 24. Morris , J. S. , Smith , K. A. , Cowen , P. J. , Friston , K. J. & Dolan , R. J . Covariation of activity in habenula and dorsal raphé nuclei following tryptophan depletion . Neuroimage 10 , 163 – 172 ( 1999 ). OpenUrl CrossRef PubMed Web of Science 25. Winter , C. , Vollmayr , B. , Djodari-Irani , A. , Klein , J. & Sartorius , A . Pharmacological inhibition of the lateral habenula improves depressive-like behavior in an animal model of treatment resistant depression . Behav. Brain Res . 216 , 463 – 465 ( 2011 ). OpenUrl CrossRef PubMed Web of Science 26. Sachs , B. D. , Ni , J. R. & Caron , M. G . Brain 5-HT deficiency increases stress vulnerability and impairs antidepressant responses following psychosocial stress . Proc. Natl. Acad. Sci. U. S. A . 112 , 2557 – 2562 ( 2015 ). OpenUrl Abstract / FREE Full Text 27. Sartorius , A. et al. Remission of major depression under deep brain stimulation of the lateral habenula in a therapy-refractory patient . Biol. Psychiatry 67 , e9 – e11 ( 2010 ). OpenUrl CrossRef PubMed Web of Science 28. Tchenio , A. , Lecca , S. , Valentinova , K. & Mameli , M . Limiting habenular hyperactivity ameliorates maternal separation-driven depressive-like symptoms . Nat. Commun . 8 , 1135 ( 2017 ). OpenUrl CrossRef PubMed 29. ↵ Yang , Y. et al. Ketamine blocks bursting in the lateral habenula to rapidly relieve depression . Nature 554 , 317 – 322 ( 2018 ). OpenUrl CrossRef PubMed 30. ↵ Zhukovskaya , A. et al. Heightened lateral habenula activity during stress produces brainwide and behavioral substrates of susceptibility . Neuron 112 , 3940 – 3956 .e10 ( 2024 ). OpenUrl PubMed 31. ↵ Tye , K. M. et al. Dopamine neurons modulate neural encoding and expression of depression-related behaviour . Nature 493 , 537 – 541 ( 2013 ). OpenUrl CrossRef PubMed Web of Science 32. ↵ Willmore , L. , Cameron , C. , Yang , J. , Witten , I. B. & Falkner , A. L . Behavioural and dopaminergic signatures of resilience . Nature 611 , 124 – 132 ( 2022 ). OpenUrl CrossRef PubMed 33. ↵ Friedman , A. K. et al. Enhancing depression mechanisms in midbrain dopamine neurons achieves homeostatic resilience . Science 344 , 313 – 319 ( 2014 ). OpenUrl Abstract / FREE Full Text 34. ↵ Geisler , S. , Andres , K. H. & Veh , R. W . Morphologic and cytochemical criteria for the identification and delineation of individual subnuclei within the lateral habenular complex of the rat . J. Comp. Neurol . 458 , 78 – 97 ( 2003 ). OpenUrl CrossRef PubMed Web of Science 35. ↵ Hashikawa , Y. et al. Transcriptional and spatial resolution of cell types in the mammalian habenula . Neuron 106 , 743 – 758 .e5 ( 2020 ). OpenUrl CrossRef PubMed 36. ↵ Sylwestrak , E. L. et al. Cell-type-specific population dynamics of diverse reward computations . Cell 185 , 3568 – 3587 .e27 ( 2022 ). OpenUrl CrossRef PubMed 37. ↵ Tiklová , K. et al. Single-cell RNA sequencing reveals midbrain dopamine neuron diversity emerging during mouse brain development . Nat. Commun . 10 , 581 ( 2019 ). 38. ↵ Saunders , A. et al. Molecular Diversity and Specializations among the Cells of the Adult Mouse Brain . Cell 174 , 1015 – 1030 .e16 ( 2018 ). OpenUrl CrossRef PubMed 39. ↵ Phillips , R. A ., 3rd et al. An atlas of transcriptionally defined cell populations in the rat ventral tegmental area . Cell Rep . 39 , 110616 ( 2022 ). 40. ↵ Zhang , L. et al. A GABAergic cell type in the lateral habenula links hypothalamic homeostatic and midbrain motivation circuits with sex steroid signaling . Transl. Psychiatry 8 , 50 ( 2018 ). 41. Aizawa , H. , Kobayashi , M. , Tanaka , S. , Fukai , T. & Okamoto , H . Molecular characterization of the subnuclei in rat habenula . J. Comp. Neurol . 520 , 4051 – 4066 ( 2012 ). OpenUrl CrossRef PubMed 42. van de Haar , L. L. et al. Molecular signatures and cellular diversity during mouse habenula development . Cell Rep . 40 , 111029 ( 2022 ). 43. ↵ Wallace , M. L. et al. Anatomical and single-cell transcriptional profiling of the murine habenular complex . Elife 9 , e51271 ( 2020 ). OpenUrl CrossRef PubMed 44. ↵ Margolis , E. B. , Lock , H. , Hjelmstad , G. O. & Fields , H. L . The ventral tegmental area revisited: is there an electrophysiological marker for dopaminergic neurons? J. Physiol . 577 , 907 – 924 ( 2006 ). OpenUrl CrossRef PubMed Web of Science 45. Lammel , S. , Ion , D. I. , Roeper , J. & Malenka , R. C . Projection-specific modulation of dopamine neuron synapses by aversive and rewarding stimuli . Neuron 70 , 855 – 862 ( 2011 ). OpenUrl CrossRef PubMed Web of Science 46. Matsumoto , H. , Tian , J. , Uchida , N. & Watabe-Uchida , M . Midbrain dopamine neurons signal aversion in a reward-context-dependent manner . Elife 5 , e17328 ( 2016 ). OpenUrl CrossRef PubMed 47. Engelhard , B. et al. Specialized coding of sensory, motor and cognitive variables in VTA dopamine neurons . Nature 570 , 509 – 513 ( 2019 ). OpenUrl CrossRef PubMed 48. de Jong , J. W. et al. A neural circuit mechanism for encoding aversive stimuli in the mesolimbic dopamine system . Neuron 101 , 133 – 151 .e7 ( 2019 ). OpenUrl CrossRef PubMed 49. Lammel , S. , Lim , B. K. & Malenka , R. C . Reward and aversion in a heterogeneous midbrain dopamine system . Neuropharmacology 76 Pt B , 351 – 359 ( 2014 ). OpenUrl CrossRef PubMed Web of Science 50. Lammel , S. et al. Input-specific control of reward and aversion in the ventral tegmental area . Nature 491 , 212 – 217 ( 2012 ). OpenUrl CrossRef PubMed Web of Science 51. Ford , C. P. , Mark , G. P. & Williams , J. T . Properties and opioid inhibition of mesolimbic dopamine neurons vary according to target location . J. Neurosci . 26 , 2788 – 2797 ( 2006 ). OpenUrl Abstract / FREE Full Text 52. Morales , M. & Root , D. H . Glutamate neurons within the midbrain dopamine regions . Neuroscience 282 , 60 – 68 ( 2014 ). OpenUrl CrossRef PubMed Web of Science 53. ↵ Bouarab , C. , Thompson , B. & Polter , A. M . VTA GABA neurons at the interface of stress and reward . Front. Neural Circuits 13 , 78 ( 2019 ). 54. Lowes , D. C. et al. Ventral tegmental area GABA neurons mediate stress-induced blunted reward-seeking in mice . Nat. Commun . 12 , 3539 ( 2021 ). OpenUrl CrossRef PubMed 55. ↵ Bouarab , C. et al. Sex-specific adaptations to VTA circuits following subchronic stress . Neuroscience ( 2023 ). 56. Cathomas , F. et al. Beyond the neuron: Role of non-neuronal cells in stress disorders . Neuron 110 , 1116 – 1138 ( 2022 ). OpenUrl CrossRef PubMed 57. ↵ Xin , Q. et al. Neuron-astrocyte coupling in lateral habenula mediates depressive-like behaviors . Cell 188 , 3291 – 3309 .e24 ( 2025 ). OpenUrl CrossRef PubMed 58. Cui , Y. et al. Astroglial Kir4.1 in the lateral habenula drives neuronal bursts in depression . Nature 554 , 323 – 327 ( 2018 ). OpenUrl CrossRef PubMed 59. ↵ Song , S. et al. Role of phospholipase Cη1 in lateral habenula astrocytes in depressive-like behavior in mice . Exp. Mol. Med . 57 , 872 – 887 ( 2025 ). OpenUrl PubMed 60. Yang , J. et al. Ventral tegmental area astrocytes modulate cocaine reward by tonically releasing GABA . Neuron 111 , 1104 – 1117 .e6 ( 2023 ). OpenUrl PubMed 61. ↵ Proulx , C. D. , Hikosaka , O. & Malinow , R . Reward processing by the lateral habenula in normal and depressive behaviors . Nat. Neurosci . 17 , 1146 – 1152 ( 2014 ). OpenUrl CrossRef PubMed 62. Lee , E. H. & Huang , S. L . Role of lateral habenula in the regulation of exploratory behavior and its relationship to stress in rats . Behav. Brain Res . 30 , 265 – 271 ( 1988 ). OpenUrl CrossRef PubMed Web of Science 63. Kim , H. J. et al. Postmortem gene expression profiles in the habenulae of suicides: implication of endothelial dysfunction in the neurovascular system . Mol. Brain 15 , 48 ( 2022 ). 64. ↵ Yang , Y. , Wang , H. , Hu , J. & Hu , H . Lateral habenula in the pathophysiology of depression . Curr. Opin. Neurobiol . 48 , 90 – 96 ( 2018 ). OpenUrl CrossRef PubMed 65. ↵ van Zessen , R. , Phillips , J. L. , Budygin , E. A. & Stuber , G. D . Activation of VTA GABA neurons disrupts reward consumption . Neuron 73 , 1184 – 1194 ( 2012 ). OpenUrl CrossRef PubMed Web of Science 66. ↵ Tan , K. R. et al. GABA neurons of the VTA drive conditioned place aversion . Neuron 73 , 1173 – 1183 ( 2012 ). OpenUrl CrossRef PubMed Web of Science 67. Stamatakis , A. M. et al. A unique population of ventral tegmental area neurons inhibits the lateral habenula to promote reward . Neuron 80 , 1039 – 1053 ( 2013 ). OpenUrl CrossRef PubMed Web of Science 68. Wang , H.-L. , Qi , J. , Zhang , S. , Wang , H. & Morales , M . Rewarding effects of optical stimulation of ventral tegmental area glutamatergic neurons . J. Neurosci . 35 , 15948 – 15954 ( 2015 ). OpenUrl Abstract / FREE Full Text 69. Al-Hasani , R. et al. Ventral tegmental area GABAergic inhibition of cholinergic interneurons in the ventral nucleus accumbens shell promotes reward reinforcement . Nat. Neurosci . 24 , 1414 – 1428 ( 2021 ). OpenUrl CrossRef PubMed 70. ↵ Bocklisch , C. et al. Cocaine disinhibits dopamine neurons by potentiation of GABA transmission in the ventral tegmental area . Science 341 , 1521 – 1525 ( 2013 ). OpenUrl Abstract / FREE Full Text 71. Cohen , J. Y. , Haesler , S. , Vong , L. , Lowell , B. B. & Uchida , N . Neuron-type-specific signals for reward and punishment in the ventral tegmental area . Nature 482 , 85 – 88 ( 2012 ). OpenUrl CrossRef PubMed Web of Science 72. ↵ Zhou , Z. et al. A VTA GABAergic neural circuit mediates visually evoked innate defensive responses . Neuron 103 , 473 – 488 .e6 ( 2019 ). OpenUrl CrossRef PubMed 73. Yoo , J. H. et al. Ventral tegmental area glutamate neurons co-release GABA and promote positive reinforcement . Nat. Commun . 7 , 13697 ( 2016 ). 74. ↵ Johnson , S. W. & North , R. A . Opioids excite dopamine neurons by hyperpolarization of local interneurons . J. Neurosci . 12 , 483 – 488 ( 1992 ). OpenUrl Abstract / FREE Full Text 75. ↵ Li , Y. et al. Rostral and caudal ventral tegmental area GABAergic inputs to different dorsal raphe neurons participate in opioid dependence . Neuron 101 , 748 – 761 .e5 ( 2019 ). OpenUrl CrossRef PubMed 76. ↵ Wakabayashi , K. T. et al. Chemogenetic activation of ventral tegmental area GABA neurons, but not mesoaccumbal GABA terminals, disrupts responding to reward-predictive cues . Neuropsychopharmacology 44 , 372 – 380 ( 2019 ). OpenUrl CrossRef PubMed 77. Root , D. H. , Mejias-Aponte , C. A. , Qi , J. & Morales , M . Role of glutamatergic projections from ventral tegmental area to lateral habenula in aversive conditioning . J. Neurosci . 34 , 13906 – 13910 ( 2014 ). OpenUrl Abstract / FREE Full Text 78. Li , J. et al. The convergence of aversion and reward signals in individual neurons of the mice lateral habenula . Exp. Neurol . 339 , 113637 ( 2021 ). 79. Congiu , M. , Trusel , M. , Pistis , M. , Mameli , M. & Lecca , S . Opposite responses to aversive stimuli in lateral habenula neurons . Eur. J. Neurosci . 50 , 2921 – 2930 ( 2019 ). OpenUrl CrossRef PubMed 80. Flanigan , M. E. et al. Orexin signaling in GABAergic lateral habenula neurons modulates aggressive behavior in male mice . Nat. Neurosci . 23 , 638 – 650 ( 2020 ). OpenUrl CrossRef PubMed 81. Groos , D. & Helmchen , F . The lateral habenula: A hub for value-guided behavior . Cell Rep . 43 , 113968 ( 2024 ). 82. ↵ Mondoloni , S. , Mameli , M. & Congiu , M . Reward and aversion encoding in the lateral habenula for innate and learned behaviours . Transl. Psychiatry 12 , 3 ( 2022 ). 83. ↵ Bagot , R. C. et al. Circuit-wide transcriptional profiling reveals brain region-specific gene networks regulating depression susceptibility . Neuron 90 , 969 – 983 ( 2016 ). OpenUrl CrossRef PubMed 84. Lorsch , Z. S. et al. Estrogen receptor α drives pro-resilient transcription in mouse models of depression . Nat. Commun . 9 , 1116 ( 2018 ). OpenUrl CrossRef PubMed 85. ↵ Bagot , R. C. et al. Ketamine and imipramine reverse transcriptional signatures of susceptibility and induce resilience-specific gene expression profiles . Biol. Psychiatry 81 , 285 – 295 ( 2017 ). OpenUrl CrossRef PubMed 86. ↵ Chung , A. S. , Miller , S. M. , Sun , Y. , Xu , X. & Zweifel , L. S . Sexual congruency in the connectome and translatome of VTA dopamine neurons . Sci. Rep . 7 , 11120 ( 2017 ). 87. ↵ Price , K. M. & Polter , A. M . Interactions of sex and stress in modulation of ventral tegmental area dopaminergic activity . Curr. Opin. Behav. Sci . 61 , 101477 ( 2025 ). 88. ↵ Kalamarides , D. J. , Singh , A. , Wolfman , S. L. & Dani , J. A . Sex differences in VTA GABA transmission and plasticity during opioid withdrawal . Sci. Rep . 13 , 8460 ( 2023 ). OpenUrl PubMed 89. Zhu , M. , Rogers , N. G. , Jahad , J. V. & Herman , M. A . Sex differences in the impact of electronic nicotine vapor on corticotropin-releasing factor receptor 1 neurons in the mouse ventral tegmental area . J. Neurosci . 43 , 3081 – 3093 ( 2023 ). OpenUrl Abstract / FREE Full Text 90. Calipari , E. S. et al. Dopaminergic dynamics underlying sex-specific cocaine reward . Nat. Commun . 8 , 13877 ( 2017 ). 91. ↵ Zhang , D. , Yang , S. , Yang , C. , Jin , G. & Zhen , X . Estrogen regulates responses of dopamine neurons in the ventral tegmental area to cocaine . Psychopharmacology (Berl .) 199 , 625 – 635 ( 2008 ). OpenUrl CrossRef PubMed 92. ↵ Simerly , R. B . Wired for reproduction: organization and development of sexually dimorphic circuits in the mammalian forebrain . Annu. Rev. Neurosci . 25 , 507 – 536 ( 2002 ). OpenUrl CrossRef PubMed Web of Science 93. ↵ Zannas , A. S. , Wiechmann , T. , Gassen , N. C. & Binder , E. B . Gene-stress-epigenetic regulation of FKBP5: Clinical and translational implications . Neuropsychopharmacology 41 , 261 – 274 ( 2016 ). OpenUrl CrossRef PubMed 94. ↵ van Doeselaar , L. et al. FKBP51 in glutamatergic forebrain neurons promotes early life stress inoculation in female mice . Nat. Commun . 16 , 2529 ( 2025 ). OpenUrl PubMed 95. ↵ Fan , Z. et al. Neural mechanism underlying depressive-like state associated with social status loss . Cell 186 , 560 – 576 .e17 ( 2023 ). OpenUrl CrossRef PubMed 96. ↵ Park , H. et al. Exposure to stressors facilitates long-term synaptic potentiation in the lateral habenula . J. Neurosci . 37 , 6021 – 6030 ( 2017 ). OpenUrl Abstract / FREE Full Text 97. ↵ Nuno-Perez , A. et al. Stress undermines reward-guided cognitive performance through synaptic depression in the lateral habenula . Neuron 109 , 947 – 956 .e5 ( 2021 ). OpenUrl PubMed 98. ↵ Kaul , D. et al. Astrocytic glutamate regulation is implicated in the development of stress-related psychiatric disorders . Neuroscience ( 2024 ). 99. Wu , J. et al. Integrating spatial and single-nucleus transcriptomic data elucidates microglial-specific responses in female cynomolgus macaques with depressive-like behaviors . Nat. Neurosci . 26 , 1352 – 1364 ( 2023 ). OpenUrl PubMed 100. ↵ Mayhew , J. , Beart , P. M. & Walker , F. R . Astrocyte and microglial control of glutamatergic signalling: a primer on understanding the disruptive role of chronic stress . J. Neuroendocrinol . 27 , 498 – 506 ( 2015 ). OpenUrl CrossRef PubMed 101. ↵ Anacker , C. et al. Role for the kinase SGK1 in stress, depression, and glucocorticoid effects on hippocampal neurogenesis . Proc. Natl. Acad. Sci. U. S. A . 110 , 8708 – 8713 ( 2013 ). OpenUrl Abstract / FREE Full Text 102. Dattilo , V. , Amato , R. , Perrotti , N. & Gennarelli , M . The emerging role of SGK1 (serum- and glucocorticoid-regulated kinase 1) in major depressive disorder: Hypothesis and mechanisms . Front. Genet . 11 , 826 ( 2020 ). 103. Miyata , S. et al. Plasma corticosterone activates SGK1 and induces morphological changes in oligodendrocytes in corpus callosum . PLoS One 6 , e19859 ( 2011 ). OpenUrl CrossRef PubMed 104. Miyata , S. et al. Sgk1 regulates desmoglein 1 expression levels in oligodendrocytes in the mouse corpus callosum after chronic stress exposure . Biochem. Biophys. Res. Commun . 464 , 76 – 82 ( 2015 ). OpenUrl PubMed 105. ↵ Hinds , L. R. et al. Dynamic glucocorticoid-dependent regulation of Sgk1 expression in oligodendrocytes of adult male rat brain by acute stress and time of day . PLoS One 12 , e0175075 ( 2017 ). OpenUrl CrossRef PubMed 106. ↵ Jin , S. et al. Inference and analysis of cell-cell communication using CellChat . Nat. Commun . 12 , 1088 ( 2021 ). OpenUrl CrossRef PubMed 107. ↵ Jin , S. , Plikus , M. V. & Nie , Q . CellChat for systematic analysis of cell-cell communication from single-cell transcriptomics . Nat. Protoc . 20 , 180 – 219 ( 2025 ). OpenUrl CrossRef PubMed 108. ↵ Omelchenko , N. & Sesack , S. R . Ultrastructural analysis of local collaterals of rat ventral tegmental area neurons: GABA phenotype and synapses onto dopamine and GABA cells . Synapse 63 , 895 – 906 ( 2009 ). OpenUrl CrossRef PubMed Web of Science 109. ↵ Theile , J. W. , Morikawa , H. , Gonzales , R. A. & Morrisett , R. A . GABAergic transmission modulates ethanol excitation of ventral tegmental area dopamine neurons . Neuroscience 172 , 94 – 103 ( 2011 ). OpenUrl CrossRef PubMed 110. ↵ Van Bockstaele , E. J. & Pickel , V. M . GABA-containing neurons in the ventral tegmental area project to the nucleus accumbens in rat brain . Brain Res . 682 , 215 – 221 ( 1995 ). OpenUrl CrossRef PubMed Web of Science 111. Taylor , S. R. et al. GABAergic and glutamatergic efferents of the mouse ventral tegmental area: Mouse VTA projections . J. Comp. Neurol . 522 , 3308 – 3334 ( 2014 ). OpenUrl CrossRef PubMed 112. Steffensen , S. C. , Svingos , A. L. , Pickel , V. M. & Henriksen , S. J . Electrophysiological characterization of GABAergic neurons in the ventral tegmental area . J. Neurosci . 18 , 8003 – 8015 ( 1998 ). OpenUrl Abstract / FREE Full Text 113. Morales , M. & Margolis , E. B . Ventral tegmental area: cellular heterogeneity, connectivity and behaviour . Nat. Rev. Neurosci . 18 , 73 – 85 ( 2017 ). OpenUrl CrossRef PubMed 114. Brown , M. T. C. et al. Ventral tegmental area GABA projections pause accumbal cholinergic interneurons to enhance associative learning . Nature 492 , 452 – 456 ( 2012 ). OpenUrl CrossRef PubMed 115. Breton , J. M. et al. Relative contributions and mapping of ventral tegmental area dopamine and GABA neurons by projection target in the rat . J. Comp. Neurol . 527 , 916 – 941 ( 2019 ). OpenUrl PubMed 116. ↵ Creed , M. C. , Ntamati , N. R. & Tan , K. R . VTA GABA neurons modulate specific learning behaviors through the control of dopamine and cholinergic systems . Front. Behav. Neurosci . 8 , 8 ( 2014 ). 117. ↵ Paul , E. J. , Tossell , K. & Ungless , M. A . Transcriptional profiling aligned with in situ expression image analysis reveals mosaically expressed molecular markers for GABA neuron sub-groups in the ventral tegmental area . Eur. J. Neurosci . 50 , 3732 – 3749 ( 2019 ). OpenUrl CrossRef PubMed 118. Margolis , E. B. , Toy , B. , Himmels , P. , Morales , M. & Fields , H. L . Identification of rat ventral tegmental area GABAergic neurons . PLoS One 7 , e42365 ( 2012 ). OpenUrl CrossRef PubMed 119. Lahti , L. et al. Differentiation and molecular heterogeneity of inhibitory and excitatory neurons associated with midbrain dopaminergic nuclei . Development 143 , 516 – 529 ( 2016 ). OpenUrl Abstract / FREE Full Text 120. ↵ Olson , V. G. & Nestler , E. J . Topographical organization of GABAergic neurons within the ventral tegmental area of the rat . Synapse 61 , 87 – 95 ( 2007 ). OpenUrl CrossRef PubMed Web of Science 121. ↵ Bubar , M. J. & Cunningham , K. A . Distribution of serotonin 5-HT2C receptors in the ventral tegmental area . Neuroscience 146 , 286 – 297 ( 2007 ). OpenUrl CrossRef PubMed Web of Science 122. ↵ Wamsley , B. & Fishell , G . Genetic and activity-dependent mechanisms underlying interneuron diversity . Nat. Rev. Neurosci . 18 , 299 – 309 ( 2017 ). OpenUrl CrossRef PubMed 123. ↵ Oriol , L. et al. Ventral tegmental area interneurons revisited: GABA and glutamate projection neurons make local synapses . Elife 13 , ( 2025 ). 124. ↵ La Manno , G. et al. Molecular Diversity of Midbrain Development in Mouse, Human, and Stem Cells . Cell 167 , 566 – 580 .e19 ( 2016 ). OpenUrl CrossRef PubMed 125. ↵ Waag , R. et al. Distinct molecular mechanisms of stress habituation in the mouse hippocampus . Neuroscience ( 2025 ). 126. ↵ Tan , W. et al. Anxiety control by astrocytes in the lateral habenula . Neurosci. Res . 205 , 1 – 15 ( 2024 ). OpenUrl PubMed 127. ↵ Bączyńska , E. et al. Stress resilience is an active and multifactorial process manifested by structural, functional, and molecular changes in synapses . Neurobiol. Stress 33 , 100683 ( 2024 ). 128. ↵ Wilkinson , M. B. et al. Imipramine treatment and resiliency exhibit similar chromatin regulation in the mouse nucleus accumbens in depression models . J. Neurosci . 29 , 7820 – 7832 ( 2009 ). OpenUrl Abstract / FREE Full Text 129. ↵ Dias , C. et al. β-catenin mediates stress resilience through Dicer1/microRNA regulation . Nature 516 , 51 – 55 ( 2014 ). OpenUrl CrossRef PubMed Web of Science 130. ↵ Minerva , A. R. et al. A practical guide to sequencing in neuropsychiatric research . NPP Digit. Psychiatry Neurosci . 3 , 21 ( 2025 ). 131. ↵ Dieterich , A. , Liu , T. & Samuels , B. A . Chronic non-discriminatory social defeat stress reduces effort-related motivated behaviors in male and female mice . Transl. Psychiatry 11 , 125 ( 2021 ). 132. ↵ Walker , H. & Frost , N. A . Protocol for the generation of single-nuclei RNA-seq libraries and quantification of heterogeneous cell types activated during social interaction . STAR Protoc . 5 , 103395 ( 2024 ). 133. ↵ Hrvatin , S. et al. Neurons that regulate mouse torpor . Nature 583 , 115 – 121 ( 2020 ). OpenUrl CrossRef PubMed 134. ↵ Stuart , T. et al. Comprehensive Integration of Single-Cell Data . Cell 177 , 1888 – 1902 .e21 ( 2019 ). OpenUrl CrossRef PubMed 135. ↵ Hao , Y. et al. Integrated analysis of multimodal single-cell data . Cell 184 , 3573 – 3587 .e29 ( 2021 ). OpenUrl CrossRef PubMed 136. ↵ Love , M. I. , Huber , W. & Anders , S . Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2 . Genome Biol . 15 , 550 ( 2014 ). 137. ↵ Ahdesmaki , M. , Zuber , V. , Gibb , S. & Strimmer , K. Sda: Shrinkage discriminant analysis and CAT score variable selection. CRAN: Contributed Packages The R Foundation doi: 10.32614/cran.package.sda ( 2008 ). View the discussion thread. Back to top Previous Next Posted November 13, 2025. Download PDF Supplementary Material Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Stress resilience is associated with transcriptional remodeling in the VTA Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Stress resilience is associated with transcriptional remodeling in the VTA Adelaide R. Minerva , Brenna McMannon , Rixing Lin , Anna Zhukovskaya , Ilana B. Witten , Catherine Jensen Peña bioRxiv 2025.11.13.687606; doi: https://doi.org/10.1101/2025.11.13.687606 Share This Article: Copy Citation Tools Stress resilience is associated with transcriptional remodeling in the VTA Adelaide R. Minerva , Brenna McMannon , Rixing Lin , Anna Zhukovskaya , Ilana B. Witten , Catherine Jensen Peña bioRxiv 2025.11.13.687606; doi: https://doi.org/10.1101/2025.11.13.687606 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Neuroscience Subject Areas All Articles Animal Behavior and Cognition (7635) Biochemistry (17690) Bioengineering (13892) Bioinformatics (41935) Biophysics (21451) Cancer Biology (18587) Cell Biology (25499) Clinical Trials (138) Developmental Biology (13377) Ecology (19899) Epidemiology (2067) Evolutionary Biology (24318) Genetics (15609) Genomics (22506) Immunology (17736) Microbiology (40394) Molecular Biology (17181) Neuroscience (88601) Paleontology (666) Pathology (2832) Pharmacology and Toxicology (4824) Physiology (7641) Plant Biology (15152) Scientific Communication and Education (2045) Synthetic Biology (4294) Systems Biology (9825) Zoology (2271)