Differential transcriptional profiles of vagal sensory neurons in female and male mice

Introduction Differences in metabolic homeostasis, diabetes, and obesity between males and females are evident in rodents and humans. Vagal sensory neurons in the vagus nerve ganglia innervate a variety of visceral organs and use specialized nerve endings to sense interoceptive signals. This visceral organ-brain axis plays a role in relaying interoceptive signals to higher brain centers, as well as in regulating the vago-vagal reflex. I hypothesized that molecularly distinct populations of vagal sensory neurons would play a role in causing differences in metabolic homeostasis between the sexes. Methods SnRNA-Seq was conducted on dissociated cells from the vagus nerve ganglia using the 10X Genomics Chromium platform. Results Single-nucleus RNA sequencing analysis of vagal sensory neurons from female and male mice revealed differences in the transcriptional profiles of cells in the vagus nerve ganglia. These differences are linked to the expression of sex-specific genes such as Xist, Tsix, and Ddx3y. Among the 13 neuronal clusters, one-fourth of the neurons in male mice were located in the Ddx3y-enriched VN1 and VN8 clusters, which displayed higher enrichment of Trpv1, Piezo2, Htr3a, and Vip genes. In contrast, 70% of the neurons in females were found in Xist-enriched clusters VN4, 6, 7, 10, 11, and 13, which showed enriched genes such as Fgfr1, Lpar1, Cpe, Esr1, Nrg1, Egfr, and Oprm1. Two clusters of satellite cells were identified, one of which contained oligodendrocyte precursor cells in male mice. A small population of cells expressed Ucp1 and Plin1, indicating that they are epineural adipocytes. Discussion Understanding the physiological implications of distinct transcriptomic profiles in vagal sensory neurons on energy balance and metabolic homeostasis would help develop sex-specific treatments for obesity and metabolic dysregulation.

It is widely recognized that there are significant differences in metabolic homeostasis, diabetes, and obesity between males and females, both in rodents and humans.Female rodents are less susceptible to dietinduced obesity, insulin resistance, hyperglycemia, and hypertriglyceridemia (Hevener et al., 2002;Hong et al., 2009;Shi and Clegg, 2009;Hwang et al., 2010;Pettersson et al., 2012;Mauvais-Jarvis, 2015;Varlamov et al., 2015).Similarly, premenopausal women are generally protected from metabolic diseases (Morselli et al., 2016), whereas men and postmenopausal women are more likely to accumulate fat in the intra-abdominal depot, which increases the risk of developing obesity-related metabolic complications (Shi and Clegg, 2009).Sex hormones, particularly estrogen, play a significant role in causing this difference in metabolic physiology between the sexes (Mauvais-Jarvis, 2015; Morselli et al., 2016;Krause et al., 2021).In addition to the contribution of sex hormones to metabolic homeostasis, it is plausible that molecularly distinct types of vagal sensory neurons may differentially integrate interoceptive signals at the peripheral level and transmit this filtered information to higher brain centers in a sex-dependent manner.In other words, the expression of genes encoding membrane receptors, ion channels, neurotransmitters, neuropeptides, and intracellular signaling proteins in vagal sensory neurons would be sex-and cellspecific.In this study, I sought to identify differences in the transcriptional profiles of cells in the vagus nerve ganglia of male and female mice.

Animals
All mouse care and experimental procedures were approved by the Institutional Animal Care Research Advisory Committee of Albert Einstein College of Medicine.All experiments were performed in accordance with relevant guidelines and regulations.This study is reported in accordance with ARRIVE guidelines.8-9 week-old Rosa26-floxed-STOP-Cas9-eGFP mice (stock# 024857) were purchased from Jackson Laboratory.Mice were housed in cages at a controlled temperature (22°C) with a 12: 12 h light-dark cycle and fed a standard chow diet with water provided ad libitum.

Single-nucleus RNA sequencing
Nuclei isolation and single-nucleus RNA sequencing (snRNA-Seq) were performed by the Singulomics Corporation (Singulomics.com,Bronx, NY, United States).Mice were euthanized with overdose of isoflurane (5% or greater).Thirty vagus nerve ganglia were collected from eight male and seven female Rosa26-eGFP f mice injected with AAVrg-Cre into the liver 4 weeks after viral injection (Hwang et al., 2024).The transcriptome data of these GFP-positive liver-projecting vagal sensory neurons were described in the separate study (Hwang et al., 2024).Immediately after tissue collection, 16 vagus nerve ganglia from males and 14 vagus nerve ganglia from females were flash-frozen.Individual tissue samples were divided into separate pools for males and females.These pooled samples were subsequently homogenized and lysed using Triton X-100 in RNase-free water to isolate the nuclei.For nuclear isolation, we followed the 10× genomics protocol, CG000124 (Rev.F).Briefly, isolated cells were centrifuged at 400 g for 5 min at 4°C and the supernatant were removed without disturbing the cell pellet.Cells were completely suspended in 200 μL lysis buffer and lysed on ice for 1 min.We added 800 μL nuclei wash and resuspension buffer containing 1% bovine serum albumin and 0.2 U/μl RNase inhibitor and centrifuged the nuclei at 500 g for 10 min at 4°C.The isolated nuclei were diluted to 700 nuclei/μl for standardized 10x capture.The 10x Genomics single cell protocol was immediately proceeded using 10x Genomics Chromium Next GEM 3' Single Cell Reagent kits v3.1 (10× Genomics, Pleasanton, CA).Libraries were sequenced using an Illumina NovaSeq 6000 (Illumina, San Diego, CA, United States).
The libraries were sequenced with ~200 million PE150 reads per sample on Illumina NovaSeq.The raw sequencing reads were analyzed with the mouse reference genome (mm10) with the addition of the human growth hormone poly(A) [hgGH-poly(A)] transgene sequence found in AAVrg-Cre using Cell Ranger v7.1.0.Results containing this information have been submitted elsewhere.Introns were included in the analyses.To further clean the data, Debris Identification using the Expectation Maximization (DIEM) program (Alvarez et al., 2020) was used to filter the barcodes (from raw_feature_bc_matrix of each sample), keeping barcodes with debris scores < 0.5, number of features (genes) > 250, and UMI count > 1,000.
Aggregation of the samples was performed with the cellranger aggr function, normalizing for the total number of confidently mapped reads across libraries.We then used the cellranger reanalyze function to analyze the samples with only the filtered, barcodes.The UMAP is generated by the cellranger reanalyze function, running at the default setting.Differential gene expression was determined using the 10x Genomics Loupe Browser (version 7).It should be noted that single-cell DE analysis carries a high risk of false-positive detection of differentially expressed genes.The significant gene test in the Loupe Cell Browser replicated the differential expression analysis in Cell Ranger.

Statistics
To determine differential gene expression, statistical tests in Cell Ranger generated p-values that were adjusted for multiple testing using the Benjamini-Hochberg procedure to control the false discovery rate (FDR).To identify and visualize enriched GO terms, the web-based application Gorilla was used (Eden et al., 2009).Fisher's exact test was used to determine if there was a significant association between the number of cells in the same cluster (GraphPad, Prism 10).Data were considered significantly different if the p-value was less than 0.05.
Sex differences in the transcriptional profiles of vagal sensory neurons I then sought to determine if there are sex differences in the transcriptional profiles of the cells in the vagus nerve ganglia.As expected, X-linked X-inactive-specific transcript (Xist) RNA, which plays a crucial role in regulating X-chromosome inactivation (Loda and Heard, 2019), was present in female mice (Figure 2A).On the other hand, the DEAD-box helicase 3 Y-linked (Ddx3y) transcript, which is located on the Y chromosome, was detected in male mice (Figure 2A).Unsupervised uniform manifold approximation and projection (UMAP) plots revealed that there were 20 clusters for both males and females, but there were noticeable differences in the number of cells in the same clusters between females and males (Figures 2B,C).In other words, the number of cells in the same clusters was not evenly distributed, with females having more cells in VN4, VN7, VN10, and VN13, whereas males had more cells in VN1, VN2, VN3, VN5, VN8, and VN12 (Figure 2D).Additionally, males had a higher number of cells in SC1, whereas females had a greater number of cells in SC2.There was a significant difference in the number of satellite cells (Fisher's exact test, p < 0.001) and Schwann cells (Fisher's exact test, p < 0.001) between the sexes.Differential gene expression analysis revealed that more than 100 genes differed significantly between females and males (Figure 2E; Supplementary Table 3).These included genes encoding voltagedependent channels, ligand-gated channels, membrane receptors, intracellular signaling proteins, and transcription factors (i.e., Kcnd2,Kcnq3,Grid2,Grik2,Rora,Gnas,Pde10a,Ndrg4,Nfia,and Tcf4; Supplementary Table 3).Upon conducting a selective gene ontology (GO) analysis using GOrilla (Eden et al., 2009) of the significantly upregulated genes in male and female mice, it was observed that among the top 20 enriched GO terms in biological process, they shared eight GO terms related to cell development and differentiation (Supplementary Figure 1).Interestingly, one-fourth of the enriched GO terms in females were related to synaptic signaling, whereas one-third of them in males were associated with intracellular component organization and transport (Supplementary Figure 1).
I further evaluated whether the uneven distribution of cell numbers in the same clusters is associated with sex differences.Gene distribution plots revealed that both Xist and X (inactive)-specific transcript (Tsix) that binds Xist during X chromosome inactivation (Lee et al., 1999) were expressed in VN4, VN6, VN7, VN10, VN11, and VN13 (Figure 3A).In contrast, Ddx3y was enriched only in VN1 and VN8 (Figure 3A).In males, the percentage of cells in VN1 and VN8 was 25% (n = 879 out of 3,569 neurons), whereas in females, it was 12% of the total number of neurons (n = 301 out of 2,591 neurons).In contrast, females had a significantly higher proportion of VN4, VN6, VN7, VN10, VN11, and VN13, accounting for 69% of the cells in the neuronal clusters, compared to only 18% in males.VN9 expressed both Xist and Ddx3y, while the expression of the three genes was very low or minimal in VN2, VN3, VN5, and VN12 (Figure 3A).The UMAP plots further demonstrated the uneven distribution of the number of cells in the same clusters (Figure 3B).Comparison of the transcriptional profiles in the neural clusters revealed that more than 1,000 genes were differentially expressed among these clusters (Supplementary Table 4).Figure 3C shows the top five enriched genes in VN1, VN4, VN6, VN7, VN8, VN10, VN11, and VN13.

Sex differences in the transcriptional profiles of satellite cells
As the number of satellite cells in SC1 and SC2 was significantly different (Fisher's exact test, p < 0.001; Figures 2D, 4A), I examined whether this was also related to sex differences.The results showed that Xist expression levels were higher in SC2 than in SC1 (Figure 4B).Moreover, there were differential gene expression between the two clusters (Figure 4C).Tubulin beta 3 (Tubb3) expression was significantly higher in SC1 than in SC2.Tubb3, a gene typically associated with neurons, has also been discovered in intermediate cells that lie between oligodendrocyte precursor cells (OPCs) and astrocytes (28).Furthermore, several other genes, including laminin subunit alpha (Lmna), CNPase (Cnp), and SRY-box transcription factor 10 (Sox10), which are enriched in OPC (Akay et al., 2021;Pruvost et al., 2023), were found in SC1 (Figure 4B).Specifically, as LMNA expression levels rise as OPCs differentiate (Pruvost et al., 2023), SC1 may consist of some progenitor cells that differentiate into myelinating Schwann cells.The number of Schwann cells was significantly lower in male mice than in female mice (Fisher's exact test, p < 0.001; Figure 2D).In contrast, EC1 and EC2 had comparable numbers of cells (Fisher's exact test, p > 0.05; Figure 4D), and both clusters expressed Xist (Figure 4E).

Discussion
This study demonstrated distinct transcriptomic profiles of cells, including sensory neurons and glial cells in the vagus nerve ganglia, between male and female mice.Prior research has established the molecular and neuroanatomical identity of target organ-specific vagal  sensory neurons in several studies (Chang et al., 2015;Williams et al., 2016;Bai et al., 2019;Makhmutova et al., 2021;Zhao et al., 2022).Although some of these studies utilized both male and female mice (Chang et al., 2015;Bai et al., 2019;Zhao et al., 2022), they reported no sex differences without specifying which results were obtained from male and female mice.Vagal sensory neurons with sophisticated nerve terminals dependent on target areas in organs are vital for detecting and relaying interoceptive signals from organs to the brain (Bai et al., 2019;Zhao et al., 2022).This implies that they may possess the capacity to regulate the energy balance.In fact, activating gastrointestinal mechanoreceptors has been shown to reduce food intake (Bai et al., 2019).Metabolic homeostasis, diabetes, and obesity exhibit notable differences between males and females in both rodents and humans (Hevener et al., 2002;Hong et al., 2009;Shi and Clegg, 2009;Hwang et al., 2010;Pettersson et al., 2012;Mauvais-Jarvis, 2015;Varlamov et al., 2015;Morselli et al., 2016).In this regard, differences in the quantity and gene expression profiles of cells in the vagus nerve ganglia between the sexes may respond distinctly to nutrients and environmental cues, leading to distinct metabolic regulation.Thirteen neuronal clusters with unique molecular characteristics were identified.Approximately half of these neuronal clusters exhibited high expression levels of both Xist and Tsix genes, whereas only two showed elevated levels of the Ddx3y gene.One of these clusters (VN9) expressed both Xist and Ddx3y and showed high levels of Prdm12, a marker gene for neurons in the jugular ganglia (Kupari et al., 2019).Another cluster (VN1) also expressed Prdm12 and had both the jugular marker Tac1 and the nodose marker Htr3a (Kim et al., 2020), suggesting that a small subset of VN1 neurons were jugular neurons.A comparison of genes that have been extensively studied in vagal sensory neurons revealed that VN1 and VN8 displayed a higher enrichment of Trpv1, Piezo2, Htr3a, Vip, Calca, P2ry1, Npy2r, and Tac1 genes (Chang et al., 2015;Williams et al., 2016;Bai et al., 2019;Makhmutova et al., 2021;Zhao et al., 2022) compared to the clusters with high expression of the Xist and Tsix genes.As they expressed chemosensitive, temperature-sensitive and mechanosensitive nociceptors, they are polymodal sensory neurons and seem to be involved in various physiological processes.PIEZO2positive neurons with mechanosensory endings played a role in the baroreceptor reflex (Min et al., 2019), while activation of P2RY1-and NPY2R-positive neurons produced opposing effects on breathing; stimulation of P2RY1-positive cells silenced respiration, while activation of NPY2R neurons induced rapid, shallow breathing (Chang et al., 2015).TRPV1-positive vagal sensory neurons also innervated the lung (Kim et al., 2020), while CGRP-positive cells innervated the stomach (Bai et al., 2019) and pancreatic islets (Makhmutova et al., 2021).Moreover, serotonin released from pancreatic β-cells activates 5-HT3R-expressing vagal sensory neurons (Makhmutova et al., 2021), and VIP-positive neurons respond to intestinal stretch (Zhao et al., 2022).
There are similarities and differences between my findings and those in others (Bai et al., 2019;Kupari et al., 2019;Prescott et al., 2020).First, highly enriched genes in nodose (Ntrk2, Eya1, and Negr1) and jugular (Gm13425, Prrxl1, and Samsn1) sensory neurons were similar between my findings and those described in the study by Kupari et al. (2019).The significance of these similarities lies in the fact that these neurons seem to have distinct functions (e.g., polymodal mechanosensor vs. noxious touch and pain).
Additionally, the study of Bai et al. (2019) identified 12 genetically distinct vagal sensory neurons that innervate both the stomach and intestine.Of particular interest is that t05 (Calca+), t06 (Vip/Uts2b), and t12 (Dbh/Ebf3) fall into the VN1 and VN8 clusters, suggesting that VN1 and VN8 neurons are neurons that innervate the gut.Moreover, VN8 neurons express Uts2b, suggesting that they may project to the intestine rather than the stomach.Additionally, the neuron cluster t09 (Caln+) lies within both VN7 and VN11, suggesting that these neurons also innervate the gut.Interestingly, VN1 and VN8 express the Ddx3y gene, while VN7 and VN11 express the Xist gene.These results imply that males and females might have distinct neurons that innervate the gut at the molecular level.
The primary distinction between the present study and earlier ones (Bai et al., 2019;Kupari et al., 2019;Prescott et al., 2020) is the presence of cells that exhibit adipocyte marker genes, including Plin1, Ucp1, and Pparg in the current study, which implies that they are epineural adipocytes.This discrepancy could be attributable to the different methodologies employed, such as single-cell RNA-sequencing vs. single-nucleus RNA-sequencing.The observed sex differences could be attributable to technical or biological factors.If the observed differences are indeed a result of technical issues, then it would be expected that there would be some overlap in the Xist and Ddx3y genes between the male and female groups.However, there is a noticeable distinction in the Xist and Ddx3y gene expression between males and females.Based on these findings, it is very unlikely that that the differences in sex observed are the result of any issues with the sample preparation, cDNA library, sequencing, or clustering.
It has been shown that sex-dependent regulation of energy balance was significantly influenced by the quantity and activity of hypothalamic neurons involved in feeding regulation.Female mice exhibit a greater number and higher activity of anorexigenic proopiomelanocortin (POMC) neurons, which helps limit the development of obesity in females (Wang et al., 2018).In this regard, the uneven distribution of neurons in certain clusters between males and females could lead to sexual dimorphism in the energy balance and metabolic homeostasis.Interestingly, VN4 appeared to have a higher expression of the fibroblast growth factor receptor 1 (Fgfr1) gene.Fibroblast growth factor 21 (FGF21), a peptide regulated by nutrient availability (Kaur et al., 2022), played a role in controlling hepatic gluconeogenesis and insulin sensitivity through FGF21 and FGFR1 signaling (Fisher et al., 2011).In addition, VN4 and VN6 showed elevated levels of lysophosphatidic acid (LPA) receptor 1 (Lpar1).Circulating LPA levels were influenced by feeding and were higher under obesogenic conditions in animal models (D 'Souza et al., 2018).The administration of LPA has been shown to have a detrimental effect on glucose homeostasis in both normal and obese mice (Rancoule et al., 2013).VN7 was found to have elevated levels of carboxypeptidase E (CPE), an enzyme responsible for processing prohormones.Hyperglycemia was observed in CPE knockout mice, with female mice exhibiting more pronounced glucose intolerance than their male counterparts (Cawley et al., 2004).VN11, on the other hand, exhibited enhanced expression of estrogen receptor 1 (ESR1), which plays a role in regulating sex-dependent metabolic balance (Mauvais-Jarvis, 2015; Morselli et al., 2016;Krause et al., 2021).While the involvement of vagal sensory neurons in these previous findings is yet unknown, it is possible that differences in gene expression patterns in vagal sensory neurons between the sexes may contribute to the sexual dimorphism in metabolic homeostasis.Hence, further studies are necessary to explore these crucial issues and to examine the physiological consequences of differences in transcriptomic profiles of vagal sensory neurons on energy balance and metabolic homeostasis.

FIGURE 1
FIGURE 1 Molecular heterogeneity of cells in the vagus nerve ganglia.(A) Schematic illustration of snRNA-Seq workflow.(B) t-SNE projection showing 20 unique clusters, comprised of 13 neuronal clusters, three glial cell clusters, one microglia cluster, and one epineural adipocytes group.(C) Heatmap plot showing differential gene expression among 20 clusters.(D) Violin plots displaying the enriched genes specific to each cell type within each specific group.

FIGURE 2
FIGURE 2 Uneven distribution of the number of cells within the same clusters.(A) t-SNE projections highlighting the expression of Xist and Ddx3y in females and males, respectively.(B,C) UMAP projections showing 20 molecularly distinct clusters in females and males.Panel (C) illustrates a clear difference in the number of cells in the same clusters between females and males.(D) Heatmap plots showing the number and percentage of cells in each cluster between females and males.(E) Heatmap plot displaying the differential expression of genes between males and females.

FIGURE 4
FIGURE 4 Existence of sex-specific glial cell populations in the vagus nerve ganglia.(A) UMAP projections showing SC1 and SC2 clusters in females and males.(B) Violin plots displaying differential expression of Xist and OPC-related genes in SC1 and SC2 clusters.(C) Heatmap plot showing the differential gene expression in SC1 and SC2 clusters.(D and E) There were no differences in the number of EC and Xist expression in females and males.