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.

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.

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.

SnRNA-Seq was conducted on dissociated cells from the vagus nerve ganglia using the 10X Genomics Chromium platform (Figure 1A). Low-quality cells and possible doublets were removed using the Cell Ranger and DIEM programs. I analyzed a total of 8,714 cells (3,791 cells from females and 4,923 cells from males). The t-distributed stochastic neighbor embedding (t-SNE) projection revealed 20 individual clusters (Figures 1B,C; Supplementary Table 1). Of the 20 clusters, 13 clusters expressed neuronal marker genes, including microtubule associated protein 2 (Map2), synapsin 1 and 2 (Syn1 and Syn2), and solute carrier family 17 member 6 (Slc17a6), indicating that they were vagal sensory neurons (VN; n = 6,160 out of 8,714 cells). These neurons were further subdivided into 13 clusters (VN1-VN13). Among the 20 clusters, three expressed glial cell marker genes: periaxin (Prx), non-compact myelin-associated protein (Ncmap), and myelin protein zero (Mpz; Figure 1D). One of these clusters was identified as myelinating Schwann cells based on the high expression of peripheral myelin protein 2 (Pmp2), ADAMTS-like 1 (Adamtsl1), and Erb-B2 receptor tyrosine kinase 4 (Erbb4) compared to the other two clusters (Supplementary Table 2). These genes are considered myelinating Schwann cell markers (Chen et al., 2006; Yim et al., 2022). Only one cluster expressed mannose receptor c-type 1 (Mrc1), a microglia (MG) marker gene (n = 267 cells). Additionally, two clusters expressed endothelial cell (EC) marker genes: platelet-derived growth factor receptor alpha (Pdgfra) and endomucin (Emcn; Figure 1D). Interestingly, a small number of cells expressed adipocyte marker genes, including perilipin 1 (Plin1), uncoupling protein 1 (Ucp1), and peroxisome proliferator-activated receptor gamma (Pparg; Figures 1C,D), suggesting that they were epineural adipocytes (EA; n = 132 cells).

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 voltage-dependent 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.

I also examined the expression of genes that are commonly present in vagal sensory neurons, specifically chemosensitive transient receptor potential cation channel subfamily V member 1 (Trpv1), temperature-sensitive TRP subfamily M member 3 (Trpm3), and mechanosensitive piezo-type ion channel component 2 (Piezo2; Figure 3D). Intriguingly, Trpv1 was mainly found in VN1 and VN8 (Figure 3D). Although the degree of Trpm3 and Piezo2 expression was different, both genes were expressed across the clusters (Figure 3D). Molecularly distinct subsets of vagal sensory neurons have their own target organs. For example, vagal sensory neurons that control the digestive system include vasoactive intestinal peptide (VIP)-, glucagon-like peptide 1 receptor (GLP1R), oxytocin receptor (OXTR), G protein-coupled receptor 65 (GPR65), somatostatin (SST), 5-hydroxytryptamine receptor 3A (HTR3A)-, and calcitonin gene-related peptide alpha (CALCA)-expressing neurons (Williams et al., 2016; Bai et al., 2019; Zhao et al., 2022). The purinergic receptor P2Y1 (P2RY1)-and neuropeptide Y receptor Y2 (NPY2R)-expressing neurons innervate the lungs (Chang et al., 2015). The pancreatic islets are innervated by vagal sensory neurons expressing tachykinin 1 (TAC1), CALCA, and HTR3 (Makhmutova et al., 2021). I thus examined if there are differences in their gene distribution across the Xist- and Ddx3y-positive VN clusters. Intriguingly, VN1 and VN8 expressed genes, including Htr3a, Vip, Calca, P2ry1, Npy2r, and Tac1, although the extent of gene distribution was different (Figure 3D). On the other hand, VN4, VN6, VN7, VN10, VN11, and VN13 displayed low or minimal expression levels of these genes (Figure 3D). As the expression levels of Glp1r, Oxtr, Gpr65, and Sst genes were minimal, the difference could not be observed.

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).