A unified rodent atlas reveals the cellular complexity and evolutionary divergence of the dorsal vagal complex
Figures

Two layers of cell identity for the murine DVC cells from snRNA-seq data.
(A) Adult C57BL/6J mice were subjected to 10X Genomics-based snRNA barcoding and sequencing after dissection of their DVC (n = 30). (B) Labeled UMAP plot of the mouse DVC (cells n = 99,740) at the lowest resolution (layer-1) including 4 cell identities, and (C) the cell class (layer-2) including 18. (D) Heatmap of the main marker genes (average log2FC expression >80 percentiles of upregulated genes) for each neuronal class of the mouse DVC (MAST algorithm; neurons n = 48,131; adj. p < 0.05). Each column is one cell and each row is one marker gene. (E) Stacked bar plot of log-normalized counts of genes for 10 transcripts (i.e. Gad2, Gad1, Slc6a5, Slc17a6, Slc18a2, Addc, Sst, Npy, Cck, and Slc5a7) related to transmission/release of 8 neurotransmitters and neuropeptides, by neuronal cell class. Each bar represents one neuron (n = 48,131). snRNA-seq = single-nuclei RNA-sequencing; DVC = dorsal vagal complex; FACS = fluorescence-activated cell sorting; UMAP = uniform manifold approximation and projection; FC = fold-change; adj. = adjusted; GABA = gamma-aminobutyric acid.

Mapping of Gad2 and VGLUT2 in neurons.
UMAP plot of the murine neurons (n = 48,131) showing scaled log-normalized expression of (A) Gad2, a gene related to GABA synthesis and transmission and (B) VGLUT2, a glutamate transporter. Cells with higher expression of Gad2 belong to the neuronal magnaclass 1, while cells with higher expression of VGLUT2 (Slc17a6) belong to the magnaclass 2 (both highlighted in yellow in adjacent plots). However, the expression is not mutually exclusive among both classes, nor does every neuron belonging to each class express the respective transcript. UMAP = uniform manifold approximation and projection; GABA = gamma-aminobutyric acid; VGLUT2 = vesicular glutamate transporter-2.

Anatomical mapping of the dorsal vagal complex (DVC) mouse layer-3 identities.
Schematic diagram of the spatial localization for the 35 neuronal identities at layer-3 of resolution in the mouse DVC. Coronal sections from the mouse brain in the Allen Brain Atlas were used to manually map the top marker genes of each (MAST algorithm on log-normalized counts; adj. p < 0.05) to the DVC. AP = area postrema; NTS = nucleus of the solitary tract; DMV = dorsal motor nucleus of the vagus; CC = central canal; 12N = hypoglossal nucleus; Cu = cuneate nucleus.

Mapping of DVC populations in a murine hypothalamic dataset.
UMAP plots of a randomly selected subset of the HypoMap database (cells n = 76,260; samples n = 32) labeled (A) by sample and (B) by original curated cell class. (C) Balloon plot showing average expression of the top 5 upregulated gene markers of each neuronal DVC class (MAST algorithm; adj. p < 0.05) in the neurons of the HypoMap subset (n = 45,679). We selected the neuronal curated class cells from our initial HypoMap subset and reprocessed them, including clustering through k-nearest neighbors algorithm using the first 30 principal components (PCs) and Harmony reduction embeddings. Resulting clusters at resolution 0.5 are shown. Average expression in the HypoMap neurons was calculated using log-normalized counts. (D) Set of UMAP plots showing scaled expression of three genes (i.e. Gfap, Gria2, and Kcnj3), which have a specific pattern of expression in murine DVC Ca+-permeable astrocytes. The astrocytic population is circled in purple (astrocytes n = 7629), and while Gfap and Gria2 expression is in agreement to what we observe in Ca+-permeable astrocytes, in the hypothalamic astrocytes, Kcnj3 is not expressed. (E) Set of UMAP plots of expression of Aldh1a1 and Stk39 which are markers of p-fibroblasts and s-fibroblasts in the DVC. The hypothalamic fibroblasts are circled in magenta (fibroblasts n = 300). DVC = dorsal vagal complex; UMAP = uniform manifold approximation and projection; AMPAr = α-amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid receptor; GIRK1 = G-protein-coupled inwardly rectifying potassium channel 1; OPC = oligodendrocyte precursor cell; NG = Ng2 glial cell.

Dorsal vagal complex (DVC)-centric dissection protocol.
(A) Following euthanasia, the brain is removed from the mouse and briefly rinsed with PBS. To better visualize the brainstem, the cerebellum is removed with fine forceps. The brain is then placed within a mouse brain matrix with 1 mm spacing between razor guides. In this example, blades will be placed in the 12th (cut 1) and 14th (cut 2) razor guide (these sites are shown). Using fresh razors, make two clean cuts in a single, smooth motion. The bottom shows the reference images of the mouse brain (from the Allen Brain Institute reference 3D mouse atlas) visualizing the location of brainstem and the AP/NTS. In these reference images, the cerebellum is not visualized. (B) Resultant brainstem tissue from the first two cuts is then further cut twice more (cuts 3 & 4) to keep the middle of the tissue where the DVC is located. This tissue is cut once more (cut 5) to generate a DVC-centric dissection. Reference images of the mouse brain (coronal sections from the Allen Brain Institute reference atlas) displaying the rostral and caudal sides of tissue collected from cuts 1 & 2 and the location of cuts 3–5 are shown.

Third layer of non-neuronal murine DVC cell identities.
(A) Labeled UMAP plot of the mouse DVC with the 15 non-neuronal cell identities shown (total cells n = 99,740; non-neuronal cells n = 51,609). (B) UMAP plot showing scaled expression of Gfap, Gria2, and Kcnj3 genes in Ca+-permeable astrocytes (circled; total cells n = 99,740; Ca+-permeable astrocytes n = 1719). (C) Barplot of log-normalized expression of six genes in oligodendrocytes (premyelinating n = 267; myelinating intermediate n = 25,692; myelinating n = 6170). Only premyelinating oligodendrocytes lack expression of Mobp, involved in myelin formation. Reduction of Synpr and Opalin expression is observed as Anln and Aff3 increase in myelinating oligodendrocytes. (D) Balloon plot of the main cell types of non-neuronal cells and neurons showing average log-normalized counts of their marker genes (MAST algorithm; adj. p < 0.05). Markers shown are upregulated genes in >80% of cells per group with average log2FC >4, or upregulated in >70% of cells with average log2FC >8. DVC = dorsal vagal complex; UMAP = uniform manifold approximation and projection; APNTS = area postrema and nucleus of the solitary tract; OL = oligodendrocyte; OPC = oligodendrocyte precursor cell; Smc = smooth muscle cells.

Mapping of DVC Ca+-permeable astrocytes in other murine brain sites.
UMAP plots of a cortical/hippocampal astrocytes dataset by (A) sample, (B) brain area, and (C) original astrocyte type (cell n = 1811). Set of UMAP plots showing scaled expression of three genes (i.e. Gfap, Gria2, and Kcnj3) which have a specific pattern of expression in murine DVC Ca+-permeable astrocytes in (D) cortical/hippocampal, (E) forebrain, (F) pons, and (G) spinal cord astrocytes. There is some overlap in Gfap and Gria2 expression in cortex, hippocampus, and forebrain, with almost no expression of Kcnj3. This is different from what we observe in DVC Ca+-permeable astrocytes. In the pons and the spinal cord, there is some expression of Kcnj3, but the Gfap and Kcnj3 co-expression is rare. DVC = dorsal vagal complex; UMAP = uniform manifold approximation and projection.

Trajectory inference of the oligodendrocyte lineage.
Dot plot of the two principal components obtained through dimensionality reduction using the main 50 principal components analysis on log-normalized gene expression data in oligodendrocytes (cells n = 32,129) (A) merged and (B) by sample, and on Harmony embeddings in oligodendrocytes (C) merged and (D) by sample. Cellular trajectory was inferred using SCORPIUS also on normalized gene expression data (A, B) and on the Harmony embeddings (C, D), and drawn across the data points (black line). Each data point represents one cell. OL = oligodendrocyte.

Murine DVC clusters.
(A) UMAP plot of the mouse DVC clusters at resolution 1 (cells n = 99,740). Cells from mice exposed to five food and aversion-related treatments were integrated using the Harmony algorithm in R, and the resulting embeddings were used to calculate the clusters and UMAP projections. We subset and re-clustered three clusters containing oligodendrocytes, vascular cells, and fibroblasts, respectively (highlighted). (B) Violin plots of the two types of fibroblasts in cluster 26, distinguished by Skt39 (stress-response ‘s-fibroblasts’) and Aldh1a1 (perivascular ‘p-fibroblasts’) expression. Each data point represents one cell (s-fibroblasts n = 244; p-fibroblasts n = 559). (C) UMAP plot of subclusters from clusters 23 (oligodendrocytes) and (D) 27 (vascular cells) using the Harmony embeddings. Gene expression markers for cell populations were mapped (Supplementary file 2), and the smallest resolution that allowed to separating them was further used to label the cells. Cells in subclusters 1 and 2 from cluster 23 were labeled ‘premyelinating_OL’ while subclusters 0 and 3 were labeled ‘myelinating_intermediate_OL’. In cluster 27, subcluster 7 was labeled ‘pericyte’, subcluster 4 was labeled ‘smooth_muscle’ and the rest of the clusters were considered to be endothelial cells. DVC = dorsal vagal complex; Res. = resolution; UMAP = uniform manifold approximation and projection; OL = oligodendrocyte; smc = smooth muscle cell.

Mouse microglial activation states.
(A) UMAP plot showing the mapping in the mouse DVC microglia (n = 2787) of monocyte (Emilin2 and Gda) and microglial (P2ry12 and Fcrls) known differential markers, corroborating the absence of monocytes in our dataset. (B) UMAP plot of the microglia from the mouse DVC showing scaled log-normalized expression of Cd5, Zbp1, Cxcl9, and Cxcl10. Microglia was subset and reprocessed using the harmony integration embeddings to obtain the UMAP coordinates. (C) Subclusters at resolution 0.2 of these microglial cells were performed and the (D) top marker genes are shown in a balloon plot with average log-normalized counts (MAST algorithm; adj. p < 0.05) along with astrocytes and oligodendrocytes (microglia n = 2787, oligodendrocytes n = 32,129, astrocytes n = 11,693). Cluster C3 shares high expression of markers with oligodendrocytes (in purple) while cluster C5 shares markers with astrocytes (in red). Microglial markers are highlighted in orange. UMAP = uniform manifold approximation and projection; DVC = dorsal vagal complex.

Neuronal populations of the murine DVC.
(A) Labeled UMAP plot with the neuronal layer-3 cell identities (neurons n = 48,131). (B) Monoamine class balloon plot showing average log-normalized expression of marker genes for each cell identity (MAST algorithm; M0 cells n = 770, M1 cells n = 551, M2 cells n = 666, and M3 cells n = 2593; adj. p < 0.05), and (C) heatmap of expression of monoamine-related genes. In the heatmap, each column represents one cell and each row, one gene. (D) Barplot of Gcg average log-normalized expression by neuronal cell identity. (E) UMAP plot showing the expression of leptin receptor and prolactin receptor genes in neurons. The GLP1 cluster is highlighted (GLP1 cells n = 153). (F) Stacked bar plot of the proportion of cells expressing one or more GABA-related genes (i.e. Slc32a1, Gad1, Gad2) and glutamate-related genes (i.e. Slc17a6, Slc17a7). The proportion of cells per cell group co-expressing GABA and glutamate associated genes is shown in purple. Only cells with log-normalized counts >0 were considered to express the genes. DVC = dorsal vagal complex; UMAP = uniform manifold approximation and projection; Res. = resolution; GABA = gamma-aminobutyric acid.

Mixed neurons subclustering.
UMAP plots with the mixed neurons of the mouse DVC (mixed neurons n = 15,190) showing their (A) original layer-3 identity and (B) final subcluster. (C) Heatmap of the top 5 marker genes (MAST algorithm; adj. p < 0.05) expression for each subcluster in all neurons (i.e. mixed neurons and all other neuronal identities, neurons n = 48,131). Each column is one cell and each row is one marker gene. Some marker genes were duplicated for more than one subcluster (e.g. Nrxn3 is a marker of subclusters 2, 3, 8, and 10; Pcdh9 is a marker of sub-clusters 6 and 8); therefore, these genes are only shown once in the plot and highlighted in red. UMAP = uniform manifold approximation and projection; DVC = dorsal vagal complex.

Th and Cck co-expression in the DVC.
(A) UMAP plot of neuronal Th and Cck scaled expression and cells with co-expression of both genes (neurons n = 48,131; Th-expressing neurons n = 764; Cck-expressing neurons n = 3821; Th/Cck co-expressing neurons n = 80). On the right plot, cells highlighted in magenta are the co-expressing neurons. The three neuronal cell identities in which the majority of nuclear co-expression of Th and Cck mRNA is found, as well as the percentage of total co-expressing neurons, are shown. (B) In situ hybridization of Cck and Th mRNA in coronal mouse DVC sections corresponding to –7.2 and –7.56 mm relative to bregma showing overlap of Th/Cck. Some of the cells with overlapping signal are highlighted with an arrow in the enhanced merged image. DVC = dorsal vagal complex; UMAP = uniform manifold approximation and projection.

The murine dorsal vagal complex (DVC) cell hierarchy.
(A) Dendrogram of the harmonized hierarchy which incorporates cells from Ludwig and our datasets. Orange cell identities are a fourth layer of cell identity resolution obtained as some of the Ludwig dataset identities are subgroups of our original cell identities at their highest resolution. Magenta labels represent layer-3 of cellular granularity. These two layers are considered high resolution. (B) UMAP plot of the integration between our and Ludwig datasets using treeArches (cells n = 171,868). (C) Pairwise heatmap showing the proportion of cells originally labeled in this study and in Ludwig dataset (y-axis), predicted to belong to each identity group (x-axis) by treeArches using our learned harmonized hierarchy. In blue are the labels considered non-specific for a high-resolution cell identity (n = 2324; 1.3%). In pink are the cell labels for rejected cells, therefore not assigned any identity by treeArches (n = 3108; 1.8%). (D) UMAP plot of the Dowsett dataset labeled through treeArches using the learned harmonized hierarchy representation from our murine DVC atlas. The ‘unlabeled’ cells include those rejected by the algorithm and thus without an assigned identity, and those with an unspecific layer-3/4 label. For example, some were labeled ‘neuron’ or ‘monoamine class’ but without a high-resolution cell identity. UMAP = uniform manifold approximation and projection; APNTS = area postrema and nucleus of the solitary tract; Lnc = long non-coding; OL = oligodendrocyte; OPC = oligodendrocyte precursor cell; Neur = neuronal; COE = collier/Olf1/EBF transcription factor.

Initial murine cell hierarchy using our samples.
(A) UMAP plot of the integration of our samples (based on the five mouse treatments administered to mice prior to euthanasia; see Methods) using treeArches. There were two cohorts of refed mice. (B) UMAP plot of the fifty layer-3 cell identities obtained using the treeArches pipeline. (C) Dendrogram showing the initial hierarchy of the three layers of cell identity established by us in this study. UMAP = uniform manifold approximation and projection; inj. = injected; APNTS = area postrema and nucleus of the solitary tract; Lnc = long non-coding; OL = oligodendrocyte; OPC = oligodendrocyte precursor cell; COE = collier/Olf1/EBF transcription factors.

Incorporation of the Ludwig dataset to our initial hierarchy using treeArches.
(A) Barplot of log-normalized expression of OPC markers (i.e. Pdgfra, Cspg4, C1ql1) and premyelinating OL markers (i.e. Bmp4, Fyn) in the OPC-labeled cells from the Ludwig dataset (total cells n = 1,770; cells expressing >1 marker gene n = 1104). Cells labeled as OPCs in the Ludwig dataset are a mix of OPC and premyelinating OL. Some cells lack expression of the five marker genes used in this plot; they are represented on the left side with an empty bar, and those expressing Fyn/Bmp4 but not expressing the OPC markers are likely premyelinating OLs. Cells were ordered using hierarchical clustering based on these five gene expression data through the Ward’s method in R. We only indicate the ‘likely cell identity’ for cells with expression of the five marker genes. (B) UMAP showing the manual mapping of Ludwig and Dowsett neuronal clusters and identities in neurons from this study. We used the gene expression markers for all clusters from Ludwig and Dowsett to identify which of our neurons shared identity with those described in both publications. Only groups and clusters successfully mapped manually are shown. The Ludwig identities are in black, Dowsett identities are in blue, and our identities are in purple. Most of the neuronal clusters in both publicly available datasets are subgroups of our cell identities. OPC = oligodendrocyte precursor cell; OL = oligodendrocyte; UMAP = uniform manifold approximation and projection; APNTS = area postrema and nucleus of the solitary tract; Lnc = long non-coding; COE = collier/Olf1/EBF transcription factors.

The snRNA-seq-derived cell identities for the rat DVC.
(A) Schematic of the pipeline for snRNA-seq of the rat DVC. (B) Labeled UMAP plots of the low resolution (i.e. layers 1 and 2 of granularity) and (C) high resolution (i.e. layer-3) cell identities in our rat DVC dataset (cells n = 12,167). We labeled 4 layer-1, 19 layer-2, and 52 layer-3 cell identities. Those labels novel to this dataset and not present in the murine data are highlighted with a pink ♦ symbol. We found a small cluster (cells n = 35) that we could not corroborate its identity at any layer of cellular granularity that we named ‘unspecific’. snRNA-seq = single-nuclei RNA-sequencing; DVC = dorsal vagal complex; FACS = fluorescence-activated cell sorting; UMAP = uniform manifold approximation and projection; APNTS = area postrema and nucleus of the solitary tract; Lnc = long non-coding; OL = oligodendrocyte; OPC = cursor cell; Neur = neuronal; COE = collier/Olf1/EBF transcription factors.

Two novel neuronal classes specific to the rat DVC.
(A) Balloon plots comparing the expression of the marker genes of the two novel neuronal classes found in rats (cells n = 12,167): immunity-akin and the ortus-akin classes (framed in gray). Expression is shown across the layer-2 rat cell identities. Based on the expression of two DVC neuronal markers Mtus2 and Syt1, we corroborated those cell classes contain neurons, and not microglial or vascular/endothelial cells (In orange and red, respectively). The immunity-akin class (n = 48 neurons) shares high expression of Cst3, Inpp5d, Csf1r, Slco2b1, and Lyn with microglial cells. Meanwhile, the ortus-akin class (n = 34 neurons) shows higher overlap of gene expression markers (i.e. Arhgap31, Pdgfra, Bcas1, Ppfibp1, and Itga9) with OPCs (highlighted with a ♦ symbol). Mouse cell identities (cells n = 99,740) do not show these expression patterns (framed in black). Average expression was calculated using log-normalized counts. (B) Violin plots of scaled expression of 10 genes (i.e. Glp1r, Gfral, Calcr, Ramp3, Gipr, Lepr, Cckar, Cckbr, Mc4r, and Prlr) coding for metabolism-associated receptors in the Pdgfra and immunity-related rat neurons (Pdgfra neurons n = 34; immunity-related neurons n = 48). Pdgfra neurons are the only layer-3 identity within the ortus-akin class. An overlapping dot plot shows one dot per cell. (C) Immunofluorescent detection of PDGFRA in mouse and rat DVC. Pink arrowheads point to cells with OPC morphology, white arrowheads point to cells with fibroblast morphology, and the yellow arrowhead points to cells with neural morphology. (D) Co-staining for PDGFRA and HuC/D in rat DVC with high magnification (right images) of area postrema PDGFRA-expressing cells. Yellow arrowheads indicate colocalization of HuC/D and PDGFRA. DVC = dorsal vagal complex; OPC = oligodendrocyte precursor cell; AP = area postrema; NTS = nucleus of the solitary tract; PDGFRA = platelet-derived growth factor receptor A; CC = central canal.

Rat DVC neurons.
(A) Labeled UMAP plot with the neuronal layer 3 cell identities in our rat dataset (neurons n = 4721). Neurons were subset and re-processed to obtain the UMAP coordinates. (B) Violin plots of scaled expression of 12 genes (i.e. Gad2, Gad1, Slc6a5, Slc17a7, Slc17a6, Th, Ddc, Nos1, Sst, Npy, Cck, and Slc5a7) related to transmission/release of GABA, glutamate, monoamines, nitric oxide, somatostatin, neuropeptide-Y, cholecystokinin, and acetylcholine in the Pdgfra and immunity-related rat neurons (Pdgfra neurons n = 34; immunity-related neurons n = 48). An overlapping dot plot shows one data point representing each cell. DVC = dorsal vagal complex; UMAP = uniform manifold approximation and projection; COE = collier/Olf1/EBF transcription factors; Lnc = long non-coding; GABA = gamma-aminobutyric acid.

Receptors and neuropeptide expression in mouse and rat neurons.
Balloon plots of average nuclear mRNA for a selected list of receptors and neuropeptides related to metabolism in mouse and rat neurons (rat n = 4721, mouse n = 49,131). Expression values represent average log-normalized counts for each gene. All cell identities at layer-3 for each species are shown.

The rodent DVC cell hierarchy.
(A) UMAP plot of the integration between the murine hierarchy and our rat dataset using treeArches (cells n = 184,035). (B) Dendrogram of the harmonized hierarchy of our labeled rat dataset and the murine DVC hierarchy. We highlight the incorporated rat cell identities (blue and magenta). Those in magenta are the novel rat identities established by us. The immunity-akin and the ortus-akin classes were not found in the murine datasets. The M3 and M5 rat identities are subgroups of a Ludwig layer-4 cell identity, therefore yielding a fifth layer within the monoamine neuronal class. UMAP = uniform manifold approximation and projection; DVC = dorsal vagal complex; APNTS = area postrema and nucleus of the solitary tract; Lnc = long non-coding; OL = oligodendrocyte; OPC = oligodendrocyte precursor cell; Neur = neuronal; COE = collier/Olf1/EBF transcription factors.

Dendrogram after incorporation of the rat dataset to the murine hierarchy.
Dendrogram (i.e. tree) resulting from incorporating our rat dataset (cells n = 12,167) to the existing murine DVC cell hierarchy using treeArches. The rat cell populations in red that did not show similarity with the previously existing murine identities are appended at the bottom by treeArches (shown in red). DVC = dorsal vagal complex.

Expression of Agrp and Hcrt in DVC neurons.
(A) UMAP plots of the neurons (n = 48,131) from our murine DVC dataset showing scaled expression of Agrp and Hcrt. (B) In situ hybridization images of Agrp and Hcrt in the DVC show no exonic mRNA. We included tissue sections of the arcuate nucleus of the hypothalamus and in the lateral hypothalamus as controls, and Agrp and Hcrt were present in these sites, respectively. For each reaction, tissue from the same mouse was used through two different sections, one including the DVC and one including the hypothalamic region. (C) UMAP plots of the Dowsett (n = 8288) and the Ludwig datasets neurons (n = 49,392) showing scaled expression of Agrp and Hcrt. (D) Genome visualization of reads coverage in the Agrp and Hcrt regions of the Mus musculus refdata-gex-mm10-2020-A reference from two of our mouse sample sequencing replicates. Resulting BAM files from the original murine DVC FASTQ files were obtained after processing. We visualized the reads from these BAM files that mapped to the Agrp and Hcrt genes using IGV v2.16.1 (Thorvaldsdóttir et al., 2013). The coverage of these genes (Agrp region: chr8:105,566,180–105,582,268; Hcrt region: chr11:100,761,693–100,762,931) is shown in gray. The majority of the reads mapped to the intronic regions of Agrp, and to the last part of the 3′-end of exon 2 and the 3′-untranslated region of Hcrt. (E) UMAP plot of rat DVC neurons (n = 4721) showing scaled expression of Agrp and Hcrt. For all UMAP projections, neurons from each dataset were subset based on our labeling for our murine and rat datasets, the treeArches resulting labeling for the Dowsett dataset, and the original labeling from Ludwig and collaborators. Then the neurons were processed using Seurat (Hao et al., 2024) and integrated using Harmony (Korsunsky et al., 2019) to yield the final UMAP embeddings. DVC = dorsal vagal complex; UMAP = uniform manifold approximation and projection; AP = area postrema; NTS = nucleus of the solitary tract; ARH = arcuate nucleus of the hypothalamus; 3v = third ventricle; ME = median eminence; LH = lateral hypothalamus. BAM = binary alignment map.

Meal-related transcriptional changes in the mouse DVC.
(A) Horizontal bar plots of the number of differential genes between refed and ad libitum fed, and refed and fasted mice in glial cells and neurons per layer-2 cell identity (MAST algorithm on log-normalized counts; adj. p < 0.05; refed neurons n = 14,396; refed glial cells n = 18,019; ad libitum fed neurons n = 9885; ad libitum fed glial cells n = 4610; fasted neurons n = 8829; fasted glial cells n = 12,496). (B) Venn diagrams of the differentially expressed genes in magnaclass 1 and 2 neurons between refed and ad libitum fed, and between refed and fasted treatments (MAST algorithm on log-normalized counts). The number of overlapping log2FC ≥1 upregulated genes between the two magnaclasses per treatment is highlighted. The treatment-induced changes in only one cell class are shown as non-overlapping. The percentage is based on the total differential genes surveyed per comparison. (C) Volcano plots of the differential genes in neurons belonging to the magnaclasses 1 and 2 between refed and ad libitum fed, and (D) refed and fasted treatments (MAST algorithm on log-normalized counts; adj. p < 0.05). Each point represents one gene. Only genes upregulated or downregulated with a log2FC ≥2 are labeled. Since we considered a minimum log fold-change of 0.1 between treatments per cell group, those genes with very low variance (i.e. log fold-change <0.1) were excluded from the differential expression analysis; therefore, a variable number of genes are shown per comparison per identity in the volcano plots. DVC = dorsal vagal complex; vr = versus; OPC = oligodendrocyte precursor cell; FC = fold-change; adj. = adjusted; NS = non-significant.

Expression of metabolism-associated receptors in the DVC magnaclasses.
Violin plots of scaled expression of 10 genes (i.e. Glp1r, Gfral, Calcr, Ramp3, Gipr, Lepr, Cckar, Cckbr, Mc4r, and Prlr) coding for metabolism-associated receptors in the magnaclass 1 and 2 neurons of the mouse DVC in fasted (magnaclass 1 n = 2756; magnaclass 2 n = 2789) and refed (magnaclass 1 n = 4586; magnaclass 2 n = 4,423). DVC = dorsal vagal complex.

Proportion of neurons in all cell identities expressing glutamatergic markers alone (dark green), Slc32a1 alone (light green), both glutamatergic markers and Slc32a1 (purple) or expressing neither Slc32a1 or glutamatergic markers (grey).

Balloon plot of Slc32a1, Gad1 and Gad2 across cell types.
The GLP1-expressing neurons express Gad2 but minimal Slc32a1.
Tables
connective | glial | neuron | unspecific | vascular | |
---|---|---|---|---|---|
Elavil3 (HuC) | 0 | 2.662964 | 2.276006 | 0 | 0.504195 |
Elavi4 (HuD) | 0.4995 | 0.375098 | 5.120204 | 0.342021 | 0.723189 |
Additional files
-
Supplementary file 1
Murine brain databases used to map our unlabeled mouse dataset.
We used four databases included in the scRNAseq v2.16.0 (Marques et al., 2016; Campbell et al., 2017; Chen et al., 2017; Romanov et al., 2017; Risso and Cole, 2025) R package and the celldex v1.12.0 built-in mouse RNA-seq reference (which we called ‘SingleR’ database) using SingleR (Aran et al., 2019) function to transfer labels to our dataset in R (R Development Core Team, 2023). The name, brain area, and reference in which each database was made public are specified.
- https://cdn.elifesciences.org/articles/106217/elife-106217-supp1-v1.xlsx
-
Supplementary file 2
Literature-derived markers per central nervous system cell identity.
We thoroughly searched in the literature for established markers at the mRNA and/or protein level. Some markers were found to be described for more than one cell population which are specified in columns ‘notes’ and ‘notes_2’. After initial labeling of our murine dataset with the databases as described, we proceeded to obtain the average expression (on log-normalized counts) of each gene coding for the markers for the labeled cells (i.e. astrocyte, endothelial, microglia, neuron, and oligodendrocyte) and unlabeled cells (i.e. unknown). We also calculated the average expression for each cluster in our dataset, established during our processing pipeline, at resolution 1.0.
- https://cdn.elifesciences.org/articles/106217/elife-106217-supp2-v1.xlsx
-
Supplementary file 3
Description of our DVC layer-3 neuronal cell identities established in rats.
Rat neurons were subset and reprocessed before labeling at layers 2 and 3 of granularity. We used neuronal class gene expression markers to identify clusters at resolution 2.0 that belonged to each neuronal cell class. These cell classes were further subclustered, and markers from Supplementary file 6 and UMAP projections were used to label each population in the dataset.
- https://cdn.elifesciences.org/articles/106217/elife-106217-supp3-v1.xlsx
-
Supplementary file 4
Literature-derived markers per central nervous system cell identity.
Breakdown of the cell identities (except neurons at layer-3 of cellular resolution) we established and corroborated using a combination of markers surveyed in Supplementary file 11 and additional markers established using the FindConservedMarkers function from Seurat v5 (Hao et al., 2024) package, indicating a cell state derived from the literature. For example, ‘s_fibroblasts’ were named after finding high expression of Stk39, a stress-related marker (Kasai et al., 2022).
- https://cdn.elifesciences.org/articles/106217/elife-106217-supp4-v1.xlsx
-
Supplementary file 5
Expression of neurotransmitter-related genes in our murine neurons.
Expression of genes coding for neurotransmitter synthesis enzymes TH, DBH, TPH, GABA, NOS, prohormone SST, and primary neurotransmitter transporters VMAT2, VGLUT1, VGLUT2, GLYT2 per neuron in our dataset. Expression and glutamate and GABA release are specified as binary variables (i.e. 0 = no expression; 1 = expression) in the green and blue columns, respectively. The percentile of the log-normalized expression is shown. ~55% of the neurons display expression of genes related to the synthesis/release of >1 primary neurotransmitter.
- https://cdn.elifesciences.org/articles/106217/elife-106217-supp5-v1.xlsx
-
Supplementary file 6
Description of our DVC layer-3 neuronal cell identities established in mice.
Neurons were subset and reprocessed before labeling at layer-3 granularity. The MAST algorithm (Finak et al., 2015) was used through Seurat v5 (Hao et al., 2024) on scaled data with a log fold-change threshold of 0.25 and detection threshold on ≥40% of cells on each neuronal cluster at resolution 1.0. We searched the resulting gene expression markers in the literature and compiled them. We then based the names for our cell identities on these genes.
- https://cdn.elifesciences.org/articles/106217/elife-106217-supp6-v1.xlsx
-
Supplementary file 7
Expression of the murine DVC identity markers by rat DVC clusters.
Average log-normalized expression of the top 5 gene expression markers of each of our murine layer-3 cell identities (except duplicated genes n = 37, and mouse-specific genes n = 27; total genes included = 186) by rat DVC cluster at resolution 2.0.
- https://cdn.elifesciences.org/articles/106217/elife-106217-supp7-v1.xlsx
-
Supplementary file 8
Labeled cells in our murine dataset.
Number of cells per treatment per identity layer-2.
- https://cdn.elifesciences.org/articles/106217/elife-106217-supp8-v1.xlsx
-
Supplementary file 9
Receptor expression in magnaclasses 1 and 2 between meal-related treatments in mice.
Differential expression analysis (MAST algorithm) (Finak et al., 2015) results for genes coding for metabolism-related receptors between meal-related treatments (i.e. refed versus ad libitum fed, and refed versus overnight fasted) and per neuronal magnaclass. If a receptor is missing for any given comparison for any of the magnaclasses, the log fold-change of that gene between both treatments in that neuronal class was <0.1; therefore, it is not reported (see Methods).
- https://cdn.elifesciences.org/articles/106217/elife-106217-supp9-v1.xlsx
-
Supplementary file 10
Filtering of our and Dowsett datasets during processing.
Compilation of statistics and cells processed for the initial high-quality cell filtering per sample.
- https://cdn.elifesciences.org/articles/106217/elife-106217-supp10-v1.xlsx
-
Supplementary file 11
Cluster filtering of our and Dowsett datasets during processing.
Compilation of low-quality cluster analysis statistics to filter high-quality clusters per sample.
- https://cdn.elifesciences.org/articles/106217/elife-106217-supp11-v1.xlsx
-
Supplementary file 12
Singlets retained per dataset after processing.
Number of doublets and singlets identified by DoubletFinder (McGinnis et al., 2019) after initial filtering of our and Dowset dataset per sample. Doublets were removed from the datasets before sample integration with Harmony (Korsunsky et al., 2019) (this step was not performed in the rat dataset as it was only one sample, n = 2 rat DVC).
- https://cdn.elifesciences.org/articles/106217/elife-106217-supp12-v1.xlsx
-
MDAR checklist
- https://cdn.elifesciences.org/articles/106217/elife-106217-mdarchecklist1-v1.docx