1. Introduction

Regulating appetite in response to changing environmental conditions to maintain energy balance requires the intricate interplay of multiple organs. Within the central nervous system, multiple discrete brain regions and cell types control appetite[1]. More specifically, the dorsal vagal complex (DVC) of the brainstem, comprising the area postrema (AP), nucleus of the solitary tract (NTS) and dorsal motor nucleus of the vagus (DMV), acts as a key nexus point for metabolic signals, evidenced by the dense expression of many relevant hormone receptors[2]. Separately, the DVC also serves as the primary site of gut-based signals via vagal afferents[3–5]. Functionally, the DVC was long considered a site promoting meal termination[6] but growing evidence suggests the DVC also regulates long-term energy balance[7, 8].

Much of the interest in DVC neurons’ role in appetite and energy balance stems from their role as therapeutic targets for obesity and anorexia in cancer cachexia[9, 10]. For many years, most pharmacotherapies failed to produce efficacious and long-lasting weight loss in obesity and increase appetite in cancer cachexia[11]. More recently however, the glucagon-like peptide 1 receptor (GLP1R) agonists are eclipsing the weight loss achieved by previous generations of therapeutics[12]. While expressed across the central nervous system, GLP1R expression is particularly enriched in DVC neurons[13] and it is the DVC that mediates many of the effects of GLP1R agonists[14]. Separately, an emerging strategy for treating cancer cachexia is inhibiting growth differentiation factor 15 (GDF15) signaling in the DVC by blocking the action of GDNF Family Receptor Alpha Like (GFRAL)[15], the receptor for GDF15[16]. Given the clinical significance of the DVC for this spectrum of conditions, understanding the molecular heterogeneity and function of this region is essential to therapeutically altering appetite while avoiding side effects.

Traditionally, identification and ensuing study of DVC cell types was limited to known neurotransmitter and receptor-expressing cell types. For instance, the cholecystokinin and monoamine-expressing cells, demarcated by tyrosine hydroxylase (Th) expression, were previously shown to be non-overlapping cell populations with unique roles in appetite suppression[17]. More recently, the advent of single-cell mRNA sequencing approaches made unbiased molecular censuses of heterogenous populations possible. Such studies focusing on the mouse have revealed the complexity of the mouse DVC, comprised of dozens of transcriptionally distinct cell types[2, 18, 19]. Moving forward, unifying information from multiple single-nuclei RNA sequencing (snRNA-seq) datasets of the DVC into a reference atlas is needed to more fully understand the cells within the DVC and generate a common and scalable classification for cell identities. Such a reference atlas is also needed to contrast the cell types found within the mouse to other research models, including the rat, which are prominent models for DVC biology[2, 20].

To address this, we generated a comprehensive snRNA-seq based atlas of the mouse and rat DVC. This atlas features mouse and rat labelled cells with increasing granularity and an accompanying computational toolbox including the mouse and rodent atlas label transfer environment using treeArches. The mouse DVC atlas reveals a higher degree of cellular complexity than previously appreciated. The cellular census of the rat DVC exhibits many similarities to the mouse but also rat-specific cells including a novel population of leptin receptor and Pdgfra-expressing neurons localized to the AP. Using this unified atlas, we also uncover the cell types that transcriptionally respond to meal consumption which may hold the potential of suppressing appetite without gastrointestinal side effects.

2. Materials and methods

2.1. Experimental design

To generate a snRNA-Seq based atlas of the mouse dorsal vagal complex (DVC), we isolated the DVC from 30 adult C57BL/6J mice and subjected them to 10x Genomics-based single-nuclei RNA barcoding and sequencing. Mice were either fasted overnight (n=5), fasted overnight and then refed 2 hours prior to euthanasia (n=10), fed ad libitum (n=5), injected with LiCl (n=5) or injected with vehicle (n=5). Conserved cell clusters across these groups were then identified.

2.2. Animals

All animals were bred in the Unit for Laboratory Animal Medicine at the University of Michigan and the Research Institute of the McGill University Health Centre. Procedures performed were approved by the University of Michigan Committee on the Use and Care of Animals, in accordance with Association for the Assessment and Approval of Laboratory Animal Care and National Institutes of Health guidelines. Alternatively, procedures were approved by the animal care committee of McGill University.

2.2.1. Mouse tissue for sn-RNA sequencing assay

Wild-type C57BL/6J mice (Jackson laboratories) were given ad libitum access to food (Purina Lab Diet 5001) and water in temperature-controlled (22°C) rooms on a 12-hour light-dark cycle with daily health checks. Mice were euthanized by isoflurane anesthesia and decapitated, then the brain was removed from the skull and aligned in a brain matrix. The cerebellar cortex was removed and the DVC was micro-dissected using AP-centric targeting and flash frozen in liquid nitrogen (Fig.S1). To limit any circadian-driven gene expression changes, all animals were euthanized 4 hours following the onset of the light phase. Frozen tissue was stored at -80°C until use.

For assessment of feeding-related gene expression, mice were either ad libitum fed or food was removed at the onset of the dark cycle and the animals were fasted overnight. Animals were then euthanized the next morning under ad libitum fed or fasting conditions. A subset of fasted animals (n=5) were provided food and allowed to refeed for 2 hours prior to euthanasia.

2.2.2. Mouse and rat tissue collection for in-situ hybridization and immunofluorescent assays

Mice were provided with ad libitum access to food (Purina Lab Diet 5001) and water in temperature-controlled (22°C) rooms on a 12 h light-dark cycle with daily health status checks. C57BL/6J mice (Jackson laboratories) or Sprague Dawley (Jackson laboratoties) were euthanized by isoflurane anesthesia followed by CO2 asphyxiation. Animals were then trans-cardially perfused with phosphate buffered saline for 3 minutes followed by 5 minutes perfusion with 10% formalin. Brain tissue was then collected and postfixed for 24 hours in 10% formalin before transfer into 30% sucrose for a minimum of 24 hours. Brains were then sectioned as 30 μm thick, free-floating sections.

2.2.3. Rat tissue collection for snRNA sequencing assays

Wild-type Sprague Dawley rats (Charles River) were provided with given access to food (Purina Lab Diet 5001) and water in temperature-controlled (22°C) rooms on a 12-hour light-dark cycle with daily health checks. Rats were euthanized by intraperitoneal injections of pentobarbital and decapitated, then the brain was removed from the skull and aligned in a brain matrix. The cerebellar cortex was removed and the DVC was dissected out in an AP-centric manner then flash frozen in liquid nitrogen. Frozen tissue was stored at -80°C until use.

2.3 Isolation of nuclei from DVC tissue

Frozen tissue was pooled (5 pooled DVC from mice and 2 DVCs from rats) and homogenized in Lysis Buffer (EZ Prep Nuclei Kit, Sigma-Aldrich) with Protector RNAase Inhibitor (Sigma-Aldrich) and filtered through a 30mm MACS strainer (Miltenyi). Filtered samples were centrifuged at 500 rcf for 5 minutes at 4°C, and pelleted nuclei were resuspended in fresh-made wash buffer (10 mM Tris Buffer, pH 8.0, 5 mM KCl, 12.5 mM MgCl2, 1% BSA with RNAse inhibitor) before undergoing a second filtration and centrifugation. The pelleted nuclei were resuspended in wash buffer with propidium iodide (Sigma-Aldrich) and underwent fluorescence-activated cell sorting (FACS) on a MoFlo Astrios Cell Sorter to remove debris. PI+ nuclei were collected and centrifuged at 100 rcf for 5 minutes at 4°C and resuspended in wash buffer to obtain a concentration of 750-1,200 nuclei/μL in preparation for sequencing.

2.4. Single-nuclei RNA-sequencing

Library preparation was handled by the Advanced Genomics Core at the University of Michigan. RT mix was added to target approximately 10,000 nuclei recovered per sample and loaded onto the 10× Chromium Controller chip. The Chromium Single Cell 3′ Library and Gel Bead Kit v3, Chromium Chip B Single Cell kit, and Chromium i7 Multiplex Kit were used for subsequent RT, cDNA amplification, and library preparation, as instructed by the manufacturer. Libraries were sequenced on an Illumina NovaSeq 6000 (pair ended with read lengths of 150 nt).

2.5. Processing of the snRNA-sequencing files

From the snRNA-seq, we obtained four FASTQ files (two per sequencing lane, containing forward and reverse sequences respectively) per mouse treatment and two files with the forward and reverse sequences for the rat samples. We pre-processed the murine FASTQ files using either CellRanger 5.0.1 with inclusion of introns or CellRanger 7.0.1 with default parameters (which include intronic sequences) to align sequencing reads to the murine genome[21]. We aligned to the Mus musculus refdata-gex-mm10-2020-A reference genome. For the rat files, we used CellRanger 3.0.1 with no introns inclusion and aligned using the Rattus norvegicus genome assembly Rnor_6.0[21]. We obtained 110,167 and 13,360 cells from the murine and rat experiments, respectively. In addition to these two datasets, we also used the following steps to process the murine Dowsett dataset[18] publicly available as a raw gene expression matrix set of files. We processed the gene expression matrices in R[22] with Seurat v5 and SeuratObject v5 to remove low quality cells[23]. For processing and analysis, R v4.3.1[22] was used. We filtered cells with ≥500 number of genes mapped and >1.22 RNA unique molecular identifiers (UMIs) counts/genes mapped ratio (based on the first 2 percentiles per sample) (Table S1). Subsequently, we merged the data per treatment, calculated the median RNA content of clusters at resolution 1.0 and those below the first quartile of the median distribution were removed (Table S2). Additionally, we used DoubletFinder v2.0.3 based on the expected doublets rates by 10X Genomics after adjusting for homotypic doublets modeled as the sum of squared annotation frequencies[24] (Table S3). On each iteration we processed the data with libraries scaled to 10,000 UMIs per cell and log-normalized. We identified the most variable genes computing a bin Z-score for dispersion based on 20 bins average expression. We regressed UMI counts and used principal component (PC) analysis for dimensionality reduction on to the top 2,000 most variable genes. We used the first 30 PCs for k-nearest neighbors clustering and for uniform manifold approximation and projection (UMAP) projections using Seurat v5[23] default parameters. Murine samples integration in a single matrix was done with Harmony[25], after which clustering and UMAP projections were performed using the harmony embeddings instead of PCs.

2.6. Mapping of our murine data to existing brain databases

We mapped the resulting murine integrated high quality singlets to the celldex v1.12.0 built-in mouse RNA sequencing reference and other four mouse databases[26–29] (Table S4) using SingleR v2.4.1[30] and scRNA-seq v2.16.0[31] packages in R[22]. Since each database has its own cell type nomenclature, in our murine dataset we established a ‘likely-cell type’ to be unanimous among all databases, for example if all databases labeled a cell as ‘neuron’. We initially labeled the cells with conflicting or unassigned label by one or more databases as ‘unknown’. In order to confirm the cell type of all the cells labeled and to identify the unknown cell identities, we mapped 473 markers from the literature. For this, the average expression per marker gene was obtained on scaled data for each cluster at resolution 1.0 (i.e. 48 clusters) (Table S5). We also manually visualized UMAP projections from all these markers in our dataset. Some clusters (e.g. cluster 27) showed in the previous step to contain cells belonging to different cell types, so we subset and reprocessed these clusters using the Harmony embeddings to assign their cell identity based on marker expression at the lowest resolution level that allowed to separate them.

2.7. Assigning high resolution clustering-based cell identity in the murine dataset

In the murine dataset, for all non-neuronal clusters at resolution 1.0 gene expression markers were obtained using the FindConservedMarkers function from Seurat v5[23] package in R[22]. The function was run using the MAST algorithm[32] on scaled data with a log fold-change (FC) threshold of 0.25. Only upregulated marker genes detected on 40% of all cells in that cluster, were retained. Since the FindConservedMarkers gives a log2FC output per sample, we calculated the average log2FC across samples and the proportion of samples for which a gene had an adjusted p-value ≥0.05. A marker gene was considered as such if it had an average log2FC>1 across samples and was significant (i.e. had an adjusted p-value ≥0.05) in more than 80% of the samples. Furthermore, neurons were subset, reprocessed and subjected to new clustering through Seurat v5[23]. We obtained the gene markers per cluster at resolution 1.0 as described. We named the clusters based on the known functions of the upregulated/downregulated genes (e.g. myelinating oligodendrocytes), peculiarities of the cell groups (e.g. Ca⁺-permeable astrocytes) or their upregulated gene expression markers if the previous methods were unsuitable (e.g. Ano2 neurons) (Table S6). All visualizations, unless otherwise stated, were done using ggplot2 v3.5.0 and ComplexHeatmap v2.18.0 packages in R[22].

2.8. Oligodendrocytes trajectory inference

We subset the oligodendrocytes (i.e. pre-myelinating, myelinating intermediate and myelinating) in our murine dataset and placed each cell on an inferred cellular trajectory using their transcriptomic data in R[22]. We used the SCORPIUS v1.0.9 package[33] with default parameters to perform dimensionality reduction by PCA on log-normalized counts and to obtain the trajectory image, with a random seed of 1,000. We additionally performed the analysis on the Harmony embeddings from the original data integration.

2.9. Gene expression files format interconversions

To convert Seurat (‘RDS’) DVC datasets to ‘h5ad’ format and vice versa, we used reticulate v1.35.0[34] and anndata 0.7.5.6[35] packages in R[22] with import of scipy, scanpy, numpy and anndata modules from Python v3.9[36]. To convert .txt, .csv and .tsv datasets (i.e. the Ludwig dataset[2], the spinal cord dataset[37] and the cortex/hippocampus database[38]), we converted those to ‘h5ad’ format in Python v3.9[36].

2.10. Murine DVC hierarchy construction

We converted our murine labeled count matrix to h5ad. We used treeArches[39] in Python v3.9[36] to create a manual tree using the three layers of cellular granularity in our database. We reprocessed the samples using the treeArches pipeline[39] to normalize the count data (counts normalized per cell and then log-normalized), identify the top 2,500 most variable genes and integrate the samples using scVI[40]. This reference latent space obtained after integration was used to generate the UMAP embeddings. We trained our manual tree based on the cell identity layer-3 labels using default parameters.

We subset the Ludwig dataset[2] to match the initial 2,500 most variable genes from our datasets and performed surgery to incorporate labels from the Ludwig dataset[2]. The two layers of granularity established by Ludwig[2] were combined in the variable “identity_layer3” to contain the highest granularity of cell identity for neurons and non-neuronal cell types. Next, we normalized the count data as specified for our murine dataset, and mapped the Ludwig dataset on the reference latent space using scArches. We then used the learn_tree function with default parameters for hierarchical progressive learning from scHPL v1.0.5[41] which yielded a hierarchy of harmonized cell labels from both datasets. We printed all trees using matplotlib v3.8.2[42] and scHPL v1.0.5[41]. All steps performed and scripts used in this phase are available in the GitHub repository: https://github.com/LabSabatini/DVC_cell_atlas.

2.10.1. Corroboration of the validity of our murine cell hierarchy using treeArches prediction modality

Using treeArches[39] and its dependencies in Python v3.9[36] we compared the original cell labels in our and Ludwig datasets with the resulting prediction of the label for each cell using our learned hierarchy. We used the predict_labels function from scHPL v1.0.5[41] with default parameters, and the latent representation obtained when we constructed our murine hierarchy. We compared the original and the predicted labels through a heatmap visualization. Next, we predicted the labels of new data, a murine dataset by Dowsett which comprises the NTS, a part of the DVC[18]. We processed the Dowsett data similarly to the Ludwig dataset and obtained the query latent representation. We then used the predict_labels function from scHPL v1.0.5[41] with default parameters to transfer the labels in our hierarchy to the Dowsett dataset based on gene expression similarity among the 2,500 initial most variable genes. We then compared the resulting labels from this prediction to the manual mapping of our cell identities based on maker expression (Table S7) using UMAP embeddings in R[22].

2.11. Labeling of our rat DVC dataset

We processed our rat dataset as previously described to obtain high-quality singlets. Based on the three layers of cell identity that we created to label the murine data from our experiments, we calculated the average expression of the top 5 marker genes for each of the layer-3 identities (except duplicated genes and mouse-specific genes) on every rat cluster to identify those that were equivalent and could be labeled as the murine data (Tables S6 and S8). Since this dataset is smaller than the mouse data, we decided to use clusters at resolution of 2.0 to capture more heterogeneity and delineate better the non-neuronal cell populations. Clusters 27 and 29 were further subset and reclustered since we found them to contain a mix of endothelial cells, mural cells, tanycytes and fibroblasts based on marker expression. Microglia and mixed neurons were sub-clustered as well. The lowest resolution level that allowed to separate populations was used to label the cells. Neurons were also subset and we mapped the murine neuronal classes markers to resolution 2.0 markers which were manually curated to delineate better the different populations within them. For all cell identities, and because we only applied one treatment to the rats (i.e. fed ad libitum), we used the FindMarkers function from Seurat v5[23] package in R[22] to obtain gene expression markers of each rat cell population. The function was run using the MAST algorithm[32] with same parameters used for the murine data (Table S9).

2.12. Construction of the rodent DVC hierarchy

We then integrated our rat labeled dataset to the murine samples in Python v3.9[36] using the initial top 2,500 most variable genes initially defined in our mouse data. Incorporation of rat labels was done as previously described for the Ludwig dataset using the learn_tree function for hierarchical progressive learning from scHPL v1.0.5[41]. The cell populations that were not rejected but placed at the bottom of the tree were moved manually to their belonging parent branches. For example, the ‘smooth_muscle_rat’ label was moved to be a children node of the mural cells node. If the parent branch did not exist, we created it. For example, the ‘immunity_related_rat’ label belonged to a new rat neuronal class, therefore, we incorporated this parent node within the neurons node and then we moved the ‘immunity_related_rat’ label to be a children node. After manual curation, we re-trained the tree. We printed all trees using matplotlib v3.8.2[42] and scHPL v1.0.5[41]. All steps performed and scripts used in this phase are available at the GitHub repository https://github.com/LabSabatini/DVC_cell_atlas for reproducibility.

2.13. Mapping gene expression in other brain-site snRNA-seq datasets

We downloaded labeled datasets from six brain regions different from the DVC: the HypoMap, spinal cord, cortex, hippocampus, pons and forebrain databases[37, 38, 43, 44]. For the HypoMap data, we subset the database by randomly selecting cells belonging to 32 included hypothalamic samples. We processed the datasets in Seurat v5[23] as previously described and obtained Harmony[25] embeddings to obtain UMAP projections. The labels used in the visualization are the original curated labels provided on each dataset. Neurons were further subset from the processed datasets and reprocessed as previously described. The other datasets were processed in the same manner after subsetting astrocytes. The pons and forebrain data was downloaded using in scRNAseq v2.16.0 packages in R through SingleR v2.4.1[30].

2.14. Differential gene expression analysis in murine treatments

We subset the murine dataset neurons and glial cells at layer 2 granularity of identity and obtained differential gene expression for pair-wise comparison of treatments (i.e. refed versus ad libitum fed, refed versus fasted and LiCl versus vehicle injected mice). We used the FindMarkers function from Seurat v5[23] package in R[22] using the MAST algorithm[32] on log-normalized counts with a minimal log fold-change threshold of 0.1 and default parameters. The number of cells per treatment can be found in Table S10.

2.15. Immunofluorescent staining

Free-floating brain sections from mouse and rat were washed with PBS, three times for five minutes. The sections were then blocked for one hour in PBS containing 0.1% Triton X-100 and 3% normal donkey serum (Fisher Scientific). The sections were incubated overnight at room temperature with Goat anti-PDGFRA (Invitrogen, CAT#, 1:200), HuC/HuD (Invitrogen, A-21271, 1:30). The following day, sections were washed and incubated for two hours with Alexflour488 and Cy3 conjugated secondary antibody (Jackson immunoresearch, 1:500). Tissues were then washed three times in PBS before being mounted onto glass slides, covered in mounting media (Fluoromount-G, Southern Biotechnology) and coverslipped. Imaging was performed using an Olympus BX61 microscope and a Zeiss LSM780-NLO laser scanning confocal microscope equipped with IR-OPO lasers at the Molecular Imaging Platform, Research Institute of the McGill University Health Centre (RI-MUHC), Montreal, CA.

2.16. In-situ hybridization and imaging

Brain sections containing the DVC from four wild-type C57BL/6J mice obtained as described above, were fixed on glass slides and stored at -20°C for no more than 36 hours. We followed the manufacturer’s ACD 323100 user manual for the RNAscope® Multiplex Fluorescent Reagent Kit v2 Assay for fixed frozen tissue samples. Briefly, we rinsed the slides with 1X PBS and incubated them at room temperature for 10 minutes after adding ∼1-2 drops of H2O2 per section. We removed the H2O2 from slides and rinsed them twice with distilled water (diH2O). Next, we submerged the slides in 200 mL of hot diH2O for 10 seconds followed by 200mL of RNAscope® 1X Target Retrieval Reagent for 5 minutes, using a steamer with lid. After briefly transferring the slides to room temperature diH2O, we submerged them in 100% ethanol for 3 minutes and allowed them to air dry. We applied approximately 1-2 drops of RNAscope® Protease III to each section and incubated them at 40°C for 30 minutes in a HybEZ^TM oven with diH2O wet paper in the tray. After rinsing the slides with diH2O, we hybridized the probes at 40°C for 2 hours.

From this point on, we performed all steps by incubating the slides at 40°C in the HybEZ^TM oven with humid tray followed by rinse using RNAscope® 1X Wash Buffer. We applied approximately 1-2 drops of RNAscope® reagents AMP1, AMP2, AMP3 and then intercalated HRP channel (15 minutes), fluorophore (30 minutes) and HRP blocker (15 minutes) for each channel present in the probe mix. The mix of probes used for the Th/Cck co-expression assay was RNAscope® Mm-Th-C2 (317621-C2) and Mm-Cck-C3 (402271-C3) diluted in probe diluent as specified in the user manual. For Hcrt we used RNAscope® Mm-Hcrt-C2 (490461-C2) and for Agrp, Mm-Agrp (400711), in a dilution as specified in the user manual. The used fluorophores were Cy3 for Cck and Hcrt (1:1,500), and Cy5 for Th and Agrp (1:1,000) diluted in RNAscope® TSA buffer. DAPI dye (1:1,000) was applied at the end of the assay, and the samples were stored for ∼48 hours at -4°C, protected from light, before imaging. We imaged the sections using a Zeiss LSM780-NLO laser scanning confocal microscope with IR-OPO lasers at the Molecular Imaging Platform at the RI-MUHC, Montreal, CA.

3. Results

3.1. The murine dorsal vagal complex cell classes

To generate a de novo snRNA-seq atlas of the mouse DVC, we isolated the DVC from 30 adult C57BL/6J mice and subjected them to 10X Genomics-based snRNA barcoding and sequencing (Fig.1A). From the CellRanger pre-processing pipeline[21] we obtained 110,167 cells of which 99,740 were retained after processing. Cell identities were initially mapped to 5 databases from different brain regions and then curated manually using 473 cell identity markers (Tables S4-S5). We distinguished clusters of neurons and glial, vascular and connective tissue cells, a first layer of cell identity (Fig.1B; Table S9). We further describe a second layer of cellular resolution, which are subgroups of the first layer that contain the cell classes (Fig.1C; Table S7).

Figure 1:

Among the ten neuronal cell classes identified at layer 2 resolution, we found two major groups which we called ‘magnaclass1’ and ‘magnaclass2’ (Fig.S2). Magnaclass1 neurons predominately express genes required for inhibitory neurotransmission (e.g. Gad1, Gad2) and share transcriptomic markers with the Sall3-class of neurons, whereas magnaclass2 neurons express mostly excitatory transcripts (e.g. the glutamate transporter VGLUT2) also found in the Lhx2-class and monoamine-class (Fig.1D-E; Fig.S2). However, as described in other brain areas[45, 46], evidence of co-transmission/co-release of multiple primary neurotransmitters was found in at least 55% of neurons regardless of their cell class (Fig.1D-E; Table S11). We also identified 8 transcriptionally distinct neural classes in addition to the magnaclasses. Among these 8 classes we found two neural classes expressing traditionally non-neuronal genes; the Golli and pH-related classes. The Golli-class neurons are a relatively small cluster, comprising ∼1% of all neurons. Interestingly, Golli-class neurons express multiple myelin-associated genes such as Mag and Mog[47–49] normally expressed by oligodendrocytes (Fig.1D; Table S6). Additionally, we identified a class of neurons with genes normally upregulated in astrocytes (e.g. Slc6a11) and related to clearance and uptake of glutamate, gamma-aminobutyric acid (GABA), Ca⁺ and bicarbonate [50, 51], which we named ‘pH-related class’ (Table S6). Similar to Golli-class neurons, pH-related neurons were a relatively rare cluster (<1% of all neurons). While the corelautus-class neurons express very few genes at higher levels, among them are Cdk8, Lars2, Gabra6 and Fat2. Expression of Cdk8 is largely absent from the DVC but is enriched within the DVC-adjacent hypoglossal nucleus (Fig.S3) suggesting these cells are hypoglossal neurons. We also identified the Lhx2-Class which was enriched for the Parvalbumin (Pvalb) gene whose expression is largely excluded from the DVC. As Pvalb is highly expressed in the neighbouring cuneate nucleus (Fig.S3), we posit this class may belong to this anatomically adjacent area which was captured in our dissection. When interrogating snRNA-seq data from the hypothalamus[43] (Fig.S4), another brain site important for energy balance, we could not find clusters sharing all markers with our neuronal classes, indicating DVC-specific neuronal programs (Fig.S4C).

3.2. The murine DVC cells at their highest resolution

To better define the heterogeneity of the mouse DVC, we further assigned a clustering-based third layer of cellular resolution resulting in fifty cell identities, of which 15 are non-neuronal (Fig.2A). At this resolution we differentiated between two major groups of astrocytes. Differential expression in synapse-related factors revealed a group of highly excitable astrocytes with high expression of the inwardly rectifying GIRK1 (Kcnj3) channel[52] and lack of expression of the AMPAr subunit GluA2 (Gria2), additionally rendering them K⁺/Na⁺/Ca⁺-permeable[53, 54] (Fig.2B). We speculate these astrocytes are specific to the DVC as we could not find Kcnj3+; Gria2-astrocytes in the hypothalamus, forebrain, cortex or hippocampus[38, 43, 44] (Fig.S4D; Fig.S5). Furthermore, although some regions close to the DVC such as pons and the spinal cord show some co-expression of Kcnj3 and Gfrap, the patterns are not to the same extent as observed in the DVC [44] (Fig.S5F-G). In addition, we distinguished a spectrum of DVC mature oligodendrocytes, including a group resembling the intermediate oligodendrocytes described in the other brain sites that decrease in number with age[37, 55] (Fig.2C). This gradient of expression of markers including Fyn, Opalin and Anln in oligodendrocytes (Fig.2C), do not respond to a developmental trajectory (Fig.S6). Although ‘myelinating-intermediate’ and ‘myelinating’ oligodendrocytes include cells with overlapping gene expression programs, there is no evidence that each oligodendrocyte will follow a single trajectory (Fig.S6). Furthermore, to properly separate some cell identities, we performed sub-clustering in three DVC clusters (i.e. clusters 23, 26 and 27) (Fig.S7A). This permitted us to label pre-myelinating oligodendrocytes, endothelial and mural cells, and distinguish two fibroblast sub-types (Fig.S7). One of these subtypes expresses high Stk39 previously described in fibroblast stress-response[56] and so was termed ‘s-fibroblasts’ for stress-associated fibroblasts (Fig.S7B), which we failed to detect in the hypothalamus (Fig.S4E). Additionally, we excluded the possibility of monocytes in our microglial data (Fig.S8A), which had no evidence of microglial activation as we mapped minimal levels of the activation genes Cxcl10, Cd5, Cxcl9 and Zbp1[57], therefore we posit that our DVC microglial cells are in a basal state (Fig.S8B). To further evaluate possible activation states in microglia, we performed sub-clustering and evaluated their marker genes (Fig.S8C). We have identified a group of microglia expressing oligodendrocyte markers and resembling a subpopulation of hippocampal microglia in an Alzheimer’s disease model which is disease-protective and synaptic-function-supporting[58], and a cortical/subcortical phagocytic subpopulation in a multiple sclerosis model[59]. A second group expresses Ndrg2 and receptors for GABA and glutamate which have been related to microglia that may induce neuronal damage[60, 61] (Fig.S8D).

Figure 2:

We also obtained the DVC gene expression markers of the main cell types in the central nervous system and neurons (Fig.2D). To more thoroughly classify DVC neural types, we subset and re-clustered the 10 neural classes from layer 2 which resulted in 35 neuronal identities as a third layer of granularity (Fig.3A; Table S6). These identities include ‘preganglionic’ or Chat-expressing neurons and a large number of neurons with unspecific markers which we called ‘mixed neurons’ (Table S6). Of the mixed neuron classes, two are largely excitatory (Mixed neurons 3 and 4) and two largely inhibitory, however, cells in the four mixed neuron classes share similar transcriptional profiles with many other neural classes. When we performed sub-clustering on mixed neurons, we could distinguish a total of 10 subtypes from the original 4 mixed neuronal identities (Fig.S9); however, transcriptional differences between these were subtle. Interestingly, a mixed neurons-like group of cells was previously described in the murine hypothalamus[26], suggesting multiple brain sites may contain neurons which lack highly differentiable transcriptional programs from snRNA-sequencing data.

Figure 3:

Additionally, we split the monoamine-class into four layer-3 cell identities (i.e. M0, M1, M2 and M3) with specific cell markers (Fig.3B; Table S6). M0 and M1 neurons express genes for enzymes synthesizing norepinephrine and serotonin, but M0 expresses minimal levels of the transporter VMAT2 (Slc18a2) (Fig.3C; Table S6). Due to the low proportion of cells expressing Th and Tph2 in M2 and M3 neurons, we consider them monoamine-modulators as they seem to synthesize trace amines (Fig.3C; Table S6). Of note, the receptor for GDF15, GFRAL, is expressed heavily by a subset of M0 neurons, some of which co-express the calcium sensing receptor (Casr) gene (Fig.3C).

We named ‘GLP1’ a neuronal cluster with high expression of the pre-proglucagon gene (Gcg) (Fig.3D; Table S6), which also expresses high levels of both the prolactin and leptin receptors[62] (Fig.3E). Additionally, these cells have the highest proportion of DVC neuronal co-expression of glutamate and GABA related genes (Fig.3F).

In our analysis we found a subset of monoaminergic neurons that express the cholecystokinin (Cck) gene (Fig.4A). While the Cck and Th-expressing DVC neurons are considered to be non-overlapping cell types with distinct physiologies[17, 63], we confirmed co-expression of Th and Cck results by in-situ hybridization (Fig.4B). The majority of this overlap occurs in cells belonging to the COE-Atp10a and monoamine class (i.e. M0, M3) cell identities (Fig.4A).

Figure 4:

3.3. The murine DVC cell atlas

To generate a comprehensive murine DVC database from multiple datasets, we utilized treeArches[39] to harmonize our snRNA-seq DVC labels with that from a publication by Ludwig and collaborators (i.e. the ‘Ludwig dataset’)[2] (Fig.5A-B). After building our initial cell hierarchy (Fig.S10), we integrated our dataset with the Ludwig dataset, giving a total of 171,868 cells (Fig.5B). Next, the labeled cell identities from the Ludwig dataset[2] were incorporated into our tree based on the similarity between the transcriptomic programs of our labeled identities and the Ludwig-labeled identities using progressive learning through scHPL[41] (Fig.5A). Some of the Ludwig dataset identities failed to be incorporated into our tree, like the oligodendrocyte precursor cells (OPCs), which contain a mix of OPCs and pre-myelinating oligodendrocytes (Fig.S11A). The mix of cell programs makes it impossible for treeArches to correctly add one branch of the Ludwig OPCs to the tree as ideally it would be added as a subpopulation of both the OPCs and pre-myelinating oligodendrocytes, which is not biologically possible. The resultant murine hierarchy has a fourth layer of cellular resolution as some of Ludwig identities, mainly neurons, are subgroups of our cell identities (Fig.S11B).

Figure 5:

We confirmed the validity of our cell hierarchy using three methods. First, by using the learned model to predict the labels of the cells we used to develop them (i.e. this study and the Ludwig dataset), for which we obtained highly concordant cell labeling (Fig.5C). Secondly, we used the learned model to predict the labels of another publicly available dataset from Dowsett and collaborators (i.e. the ‘Dowsett dataset’)[18] (Fig.5D). Using this approach, we were able to successfully label 95% of these cells at high resolution (i.e. layer 3 or 4 of cell identity granularity); further validating our combined DVC labels. Finally, we manually mapped the marker genes from Ludwig[2] and Dowsett[18] neuronal clusters in our dataset and confirmed the sub-groups of neurons among our identities pointed out by treeArches (Fig.S11B).

3.4. An integrated rodent DVC cell hierarchy

To identify rat DVC cell clusters and subclusters and contrast these with the mouse DVC atlas, we generated a de novo ad libitum-fed rat dataset (cells n=12,167) (Fig.6A) and labeled it by mapping the top 5 gene expression markers of each of our murine layer-3 cell identities (Table S8). After manual curation at three granularities (Fig.6B-C; Table S9), we obtained 52 DVC rat identities at high cellular resolution (Fig.6C) including a small group of ‘unspecific’ cells (n=35) that we could not confidently label because its top gene expression markers have unknown functions and none of our known identities’ genes clearly marked this population (Fig.6B-C; Table S8).

Figure 6:

We established 10, layer-3 neuronal identities not previously found in the murine data (Fig.6C; Fig.S12A; Table S9). Notably, two of those did not belong to any of the existing murine neuronal classes (i.e. layer-2 of cellular granularity), therefore we created two new layer-2 classes (Fig.6B). One of these, the immunity-related neurons, have high cystatin C, Inpp5d and Csf1r expression, genes commonly associated with immunity-related processes in microglia and with immunomodulation in cancer[64–66] (Fig.7; Table S10). The second neural class identified in rat but not mouse DVC were the ‘Pdgfra neurons’ (Fig.7). These neurons express Pdgfra, Arhgap31 and Itga9, known markers for progenitor cells for which we named this cell class ‘ortus-akin’ (Latin: ‘origin’) because they are similar to cells giving origin to other cells[67] (Fig.7A; Table S10). Interestingly, this neuronal class not only expressed Pdgfra, but also other OPC signature genes (Fig.7A). Mapping the expression levels of metabolic-associated receptors revealed that both immunity-akin and ortus-akin express GFRAL, prolactin receptor and leptin receptor (Fig.7B), although they seem to have a varied neurotransmitter profile (Fig.S12B). To confirm the presence of Pdgfra-expressing neurons in the rat DVC, we performed immunostaining in the mouse and rat brainstem. Within the mouse, we found PDGFRA-immunoreactivity(IR) in the two morphologically distinct cells; matching OPC and blood vessel-associated fibroblast morphologies(Fig.7C). However, within the rat, we identified PDGFRA-IR in cells with three morphologies, one of which resembled neurons with an enlarged cell body and protruding cellular processes (Fig.7C). To determine if these PDGFRA-IR cells were neurons, we co-stained with the neural marker, HuC/D and noted several PDGFRA and HuC/D copositive cells (Fig.7D); confirming our sequencing data that a PDGFRA-expressing neural class exists in the rat DVC and localizing these cells to the area postrema.

Figure 7:

To construct an integrated rodent DVC cell atlas, we incorporated the labeled rat dataset into our hierarchy using treeArches[39] (Fig.8). The rodent data was projected onto the mouse datasets (Fig.8A) and the 42 high resolution rat cell identities were appended to the tree (Fig.8B; Fig.S13). The final hierarchy has 123 labels of which 99 are high resolution according to our cellular granularity (i.e. layers 3 to 5) (Fig.8B). In addition it has 20 layer-2 identities and 7 cell identities that are novel from the rat dataset (Fig.8B).

Figure 8:

3.5. Hcrt and Agrp expression in the DVC

Since Agrp-expressing cells have been recently described within the DVC[68], we mapped this and other neuropeptide genes in DVC neurons. Surprisingly, we found Hcrt and Agrp mRNA in neurons in our dataset (Fig.S14A). We then performed in-situ hybridizations for these transcripts in both the mouse hypothalamus and DVC. While both Hcrt and Agrp were readily detectable in the hypothalamus, we failed to observe signal within the DVC (Fig.S14B). Since we also found expression of Hcrt and Agrp in neurons of the Dowsett and Ludwig datasets[2, 18] (Fig.S14C), we then mapped the sequenced reads from these transcripts and confirmed that most of them mapped to intronic and non-coding regions of these genes (Fig.S14D). Furthermore, our rat DVC dataset was pre-processed excluding introns from the reads mapping, and the Hcrt and Agrp expression levels observed is minimal (Fig.S14E), which is in agreement with the non-coding location of the reads in mouse DVC. These data explain absence of signal from our in-situ hybridizations and suggests DVC cells do not produce HCRT or AGRP neuropeptides.

3.6. Meal-related transcriptional programs in the murine DVC

As the DVC acts as a primary node for meal-related signals[69], we wanted to define which of our cell identities responded to meal consumption and thus performed differential gene expression analysis between DVC cell types from mice that were euthanized under fasting conditions, with ad libitum access to food or two hours post-refeeding following an overnight fast (Fig.1A; Fig.9). This analysis revealed widespread transcriptional changes across nearly every neural cell identity, the exceptions being corelautus, Golli class and pH classes (Fig.9A). Additionally, the size of the transcriptional changes were relatively small, with approximately 90% of the differentially expressed genes having increases or decreases between 0.5 and 1.0 log2 fold-change, with only ∼10% of the gene changes being ≥1 log2 fold-change between conditions (Fig.9A). Of the responding neurons, magnaclasses 1 and 2, which represent large populations of predominately inhibitory and excitatory neurons (Fig.S2), showed the greatest number of differentially expressed genes following refeeding for both comparisons refed versus fasting, and refed versus ad libitum food access (Fig.9A). To understand how refeeding transcriptionally alters inhibitory compared with excitatory DVC neurons, we compared differentially expressed genes between these two classes across conditions. Intriguingly, despite being in transcriptionally distinct classes, we found meal consumption altered many of the same genes (Fig.9B-D), suggesting these cells may be receiving similar postprandial signals which are propagated via these cells in unique manners. This is supported by the similar receptor expression profile of magnaclass 1 and 2 (Fig.S15). One notable exception is Cckbr which is upregulated in magnaclass 2 and downregulated in magnaclass 1 in refeeding when compared with ad libitum fed mice(Table S12). We also note that oligodendrocytes contained the highest number of differentially expressed genes in refeeding (Fig.9A). While the roles of oligodendrocytes in refeeding responses within the DVC have yet to be determined and the biologic significance of these changes are presently unclear, DVC oligodendrocytes also have robust transcriptional responses in response to fasting[18]; suggesting these cells are sensitive to both positive and negative energy balance associated-cues and may have important roles in regulating energy homeostasis..

Figure 9:

4. Discussion

The recent success of DVC-based therapies[12, 70] has moved DVC biology to the forefront of weight control therapies and there is now clear clinical impetus to fully define DVC cell types and their functions. Previous studies detailing DVC cell types describe a high degree of heterogeneity[2, 18, 19], however, these studies lacked an in-depth analysis of neurotransmitter profiles of neurons, any cross-species evaluation of DVC cell types, or an analysis of how meal consumption transcriptionally regulates these cells. We therefore developed de novo mouse and rat DVC cell atlases to address these points.

To more fully describe the cellular complexity within the DVC, we combined a pipeline of cell label transfer using brain existing published datasets and thorough curation of cell identities, which resulted in a 4-layer mouse cell hierarchy and a 5-layer rodent cell hierarchy of the DVC. Within this hierarchy we found several unique features of DVC glial cells. First, we found highly excitable Ca⁺-permeable astrocytes in the DVC in which Gfap is almost exclusively expressed and they coexpress Kcnj3 (GIRK1 potassium channel). These cells largely lack Gria2, consistent with findings in Bergmann glial cells in which colocalization of Gria2 and GFAP was found to alter neuron-glial interactions[71]. Comparing Ca⁺-permeable DVC astrocytes to hypothalamic astrocytes[43], we found the same expression pattern with minimal overlap between Gria2 and Gfap. However, DVC astrocytes seem to have some unique properties when compared with hypothalamic, cortical and hippocampal astrocytes. Principally DVC astrocytes express the GIRK1 channel, likely rendering them more excitable. Additionally, while we identified hypothalamic fibroblasts which express high Aldh1a1 similar to the DVC p-fibroblasts, we failed to find a group equivalent to s-fibroblasts. Other cell populations seem to be consistent to what has been found in other brain sites, like the intermediate oligodendrocytes that are found enriched in the hypothalamus and the corpus callosum in mice[55]. In summary, these findings show certain DVC glial cells share identities with glial cells from other CNS sites, there are DVC-specific glial cells.

Within the mouse neural classes, we found overlap between Cck and Th-expressing cells. As these cells have previously been described as non-overlapping, this was unexpected[17]. Our in-situ hybridization confirmed the presence of Cck/Th positive cells within the NTS and AP. We believe the discrepancy may be due to detection methods as these cells were initially defined as non-overlapping using immunofluorescent detection of TH protein and a Cck genetic reporter[17]. Functional studies of Cck- and Th-expressing neurons within the NTS demonstrate these cells have somewhat distinct roles; while both reduce appetite, Cck neurons are strongly aversive whereas Th-labelled neurons are non-aversive[72]. While no functional characterization of Cck and Th neurons of the AP has been performed, we posit that as some Th-expressing AP neurons also express Gfral and Casr; both having potently aversive ligands (e.g. GDF15, deoxynivalenol)[16, 73]; they are likely to mediate aversive anorexia. The role of dual Cck/Th positive neurons remains to be determined.

While AGRP expression has long been considered to be hypothalamic-exclusive, a recent report describes Agrp-expressing cells within the mouse DVC[68], i n part from the detection of Agrp mRNA in scRNA-Seq data. Similarly, we also detected transcripts for both Agrp and Hcrt in our snRNA-Seq data but failed to detect any transcripts within the DVC by in-situ hybridizations. One possible explanation for the discrepancy between the snRNA-Seq data and our in-situ hybridizations data is that the snRNA-Seq mRNA reads forAgrp and Hcrt gene overwhelmingly mapped to non-coding regions, of the respective genes. Combined, these data suggests to us that neither AGRP or HCRT neuropeptides are produced by DVC cells in rodents.

We show the presence of neurons with seemingly mixed programs in the DVC, which we called ‘mixed neurons’. Similar groups have been described in the hypothalamus[26] suggesting that some neurons may have hybrid-functionalities. These mixed neurons share many transcriptomic markers expression with many other cell identities, and may require the use of other technologies to fully identify them anatomically in the brain.

To perform a cross-species analysis, we generated a rat-based DVC atlas and compared it with the mouse DVC. Our rat dataset was smaller than the mouse data, it is possible some of the of the new neuronal cell identities (e.g. M_sema, M4) may contain cells belonging to multiple cell types, as treeArches did not appended these identities to the tree. Despite this, we corroborated two new cell types in the rat DVC that are not present in mice the Pdgfa-expressing (ortus-akin class) and the immunity-related neurons. Pdgfra is a well described marker gene for OPCs and fibroblasts within the central nervous system[67] but, to our knowledge, has not been associated with neurons. Notably, although the PDGFRA immuno-signal in HuC/D+ neurons is clearly discernible, it is less intense than in OPCs and fibroblasts, in agreement with expression levels from our snRNA-Seq data. Since it was previously demonstrated that OPCs can give rise to neurons in the adult brain[74] and there is evidence of neurogenesis in the AP following amylin administration in rats[75], it is possible that these PDGFRA+/HuCD+ cells in the AP represent a transient stage of PDGFRA+ OPCs differentiating into neurons.

Regardless of their maturity, we speculate that the Pdgfra neurons (ortus-akin class) have roles in appetite regulation as they express leptin receptor and Gfral. Interestingly, while multiple studies have failed to detect Lepr via Cre-based reporters[62, 76] or leptin-induced phosphorylated STAT3 in the AP of mice[77], Lepr and leptin-induced phosphorylated STAT3 are detectable within the rat AP[78–80]. We believe one explanation for this species discrepancy in Lepr expression is the absence of the ortus class neurons in mice. Future studies are needed to characterize these cells and determine whether their presence is found in other mammals including humans.

DVC-targeting obesity therapeutics, including amylin and GLP1 mimetics, produce clinically relevant weight loss but are associated with multiple gastrointestinal adverse events and there is a need to develop new obesity therapies that can reduce appetite without unwanted gastrointestinal side effects[81, 82]. Indeed, the DVC cells that regulate appetite suppression and nausea are separable[72], making targeting such cells an attractive target for pharmacotherapy. As meal consumption is generally not associated with gastrointestinal distress, we sought to identify refeeding responsive cells within the mouse DVC. From this analysis, we found nearly all neural cell identities are transcriptionally altered by meal consumption. We also found the magnitude of transcriptional changes was relatively small, regardless of the cell class examined. While relatively small, the scale of the transcriptional changes in response to refeeding are similar to many of the gene expression changes observed in Agrp-neurons following 16 hours food deprivation[43]. Regardless of whether this is a technical limitation of snRNA-Seq methods or accurately reflects the scale of transcriptional changes in the DVC, utilizing Therefore, identifying neurons that promote non-aversive anorexia via the transcriptional response to refeeding to identify neural classes as targets for as the basis of pharmacotherapy is unlikely to yield specific and targetable pathways.

In general, one limitation in our study is that it is based on nuclear RNA which may reflect a fraction of biological changes. Using nuclear isolation facilitates generating to have single-neuron data[83] since obtaining viable whole neurons (including all their processes) is very challenging. Previous studies on the relationship between nuclear and cellular RNA content in brain show that nuclear sequencing reads may map to intronic regions of genes in higher proportion[84]. Although this does not appear to impact gene expression analysis for neurons in the current pipelines like the ones we used for our mouse data[84, 85], it is possible that some biological features such as microglial activation state, are not fully recovered using nuclear RNA[86]. Although we performed sub-clustering in the mouse microgia and assessment of the resulting marker genes for each group, we did not find evidence of activation possibly due to this matter. Separately, for meal consumption, we were unable to analyze responses in individual mice as pooling DVCs from individual mice was the most amenable option to maximize nuclei for in our pipeline. As we pooled DVCs from both sexes, this limited our ability to detect sexual dimorphic cell class proportions and gene expression levels.

5. Conclusion

As the amount of snRNA-Seq data performed on the DVC continues to grow, there is a need for a common, scalable label set that can be applied to new studies and integrate new sequencing data. To develop such a label set, we started with a de novo DVC atlas, we generated three layers of cellular granularity for both neural and non-neuronal cell types. Then, using the treeArches pipeline[39] we harmonized our snRNA-Seq labels with the Ludwig dataset[2] and increased our cellular granularity from three to four layers of resolution. Manually curating and harmonizing our dataset enhanced the ability of our cell atlas to capture the complexity of the DVC neurons since not all groups seem to contain easily delineable subgroups of cells. For example, eleven of the labels incorporated into our cell hierarchy in this pipeline belong to the monoaminergic neuronal class, which seem to encompass multiple subgroups of specialized cells. We further constructed the rodent DVC cell atlas which resulted in a fifth layer of cell identity to the monoaminergic M3 neurons. Although the number of identities increases as layers are more specific, this gives this unified atlas greater adaptability to different research needs. Additionally, our labels can be accurately harmonized with new databases regardless of size and new DVC snRNA-Seq data can be added to the tree to generate new branches. Thus, our unified rodent DVC atlas is highly amenable to future research elucidating new biologic and therapeutic insights into DVC cell types.

Data statement

All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. FASTQ files have been deposited in the NCBI Sequence Read Archive (Accession: PRJNA1161292). In addition, our murine and rat datasets have been made available for visualization and analysis through the Broad Institute Single Cell Portal (https://singlecell.broadinstitute.org/single_cell/study/SCP2773). The labels from our cell atlas can be transferred to new datasets through treeArches using the available hosted Jupyter Notebook in Google Colab for the mouse data (https://colab.research.google.com/drive/1v983gcQxIwBN-4lGx18JXPweZoKqMbwX#scrollTo=SJTZppZWa-63), and for the rodent data (https://colab.research.google.com/drive/10tAmfRosMscuBak6wHX80TgdEBsHeWQB#scrollTo=941ab5c1-8910-4e67-a632-4d18601d45b4). Original code used for constructing the murine and rodent cell hierarchies and predicting labels using treeArches are available in the GitHub repository https://github.com/LabSabatini/DVC_cell_atlas.

Acknowledgements

We would like to thank authors of the Ludwig and the Dowsett publications for making publicly available their datasets. We thank M. Myers Jr. for the valuable help provided; S. Hebert and S. Bailey for their insights and suggestions in regards to the bioinformatic analysis; members of the Kokoeva lab for their insights on the in-situ hybridization and immunohistochemistry pipelines; S. Feng and M. Fu for technical assistance for images acquisition at the Molecular Imaging Platform at the RI-MUHC, and A. Duensing for assistance with relevant file transfer.

Additional information

Author contributions

Conceptualization: CH and PS

Methodology: CH, AB, LM, MK and PS

Sample collection: AB and PS

Software: CH, LM and HJM

Data analysis: CH, LM, MK and PS

Data curation: CH and HJM

Validation of findings: CH, FS, HJM and PS

Imaging: CH and PS

Visualization: CH, MK and PS

Supervision: CH and PS

Writing—original draft: CH and PS —review and editing: CH, AB, LM, FS, HJM, MK and PS

Funding acquisition: MK and PS

Funding

This research was also supported by grants from Canadian Institutes for Health Research (PJT180590 for PVS and 202010PJT for MK) the Natural Sciences and Engineering Research Council of Canada (RGPIN-2022-03390). HJM was funded in part by the Canada First Research Excellence Fund and Fonds de recherche du Quebec, awarded to the Healthy Brains, Healthy Lives initiative.

Abbreviations

• 3v: third ventricle

• adj.: adjusted

• AMPAr: alpha-amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid receptor

• AP: area postrema

• APNTS: Area postrema and nucleus of the solitary tract

• ARH: arcuate nucleus of the hypothalamus

• BAM: binary alignment map

• CC: central canal

• COE: collier/Olf1/EBF transcription factor

• diH2O: distilled water

• DMV: dorsal motor nucleus of the vagus

• DVC: dorsal vagal complex

• FACS: fluorescence-activated cell sorting

• FC: fold-change

• GABA: gamma-aminobutyric acid

• GDF15: growth differentiation factor 15

• GIRK1: G protein-coupled inwardly rectifying potassium channel 1

• GLP1R: glucagon-like peptide 1 receptor

• GFRAL: GDNF family receptor alpha-like

• IR: immunoreactivity

• LH: lateral hypothalamus

• Lnc: long non-coding

• ME: median eminence

• Neur: neuronal

• NG: Ng2 glial cell

• NS: non-significant

• NTS: nucleus of the solitary tract

• OL: oligodendrocyte

• OPC: oligodendrocyte precursor cell

• PC: principal component

• PDGFRA: Platelet-derived growth factor receptor A

• Res.: resolution

• RNA-seq: RNA sequencing

• smc: smooth muscle cells

• snRNA-seq: single-nuclei RNA sequencing

• UMAP: uniform manifold approximation and projection

• UMI: unique molecular identifiers

• v: version

• VGLUT2: vesicular glutamate transporter-2

• vr: versus

Additional files

Supplemental tables

Supplemental figures