Abstract
The ventrolateral medulla (VLM) is a crucial region in the brain for visceral and somatic control, serving as a significant source of synaptic input to the spinal cord. Experimental studies have shown that gene expression in individual VLM neurons is predictive of their function. However, the molecular and cellular organization of the VLM has remained uncertain. This study aimed to create a comprehensive dataset of VLM cells using single-cell RNA sequencing in male and female mice. The dataset was enriched with targeted sequencing of spinally-projecting and adrenergic/noradrenergic VLM neurons. Based on differentially expressed genes, the resulting dataset of 114,805 VLM cells identifies 23 subtypes of neurons, excluding those in the inferior olive, and five subtypes of astrocytes. Spinally-projecting neurons were found to be abundant in seven subtypes of neurons, which were validated through in situ hybridization. These subtypes included adrenergic/noradrenergic neurons, serotonergic neurons, and neurons expressing gene markers associated with premotor neurons in the ventromedial medulla. Further analysis of adrenergic/noradrenergic neurons and serotonergic neurons identified nine and six subtypes, respectively, within each class of monoaminergic neurons. Marker genes that identify the neural network responsible for breathing were concentrated in two subtypes of neurons, delineated from each other by markers for excitatory and inhibitory neurons. These datasets are available for public download and for analysis with a user-friendly interface. Collectively, this study provides a fine-scale molecular identification of cells in the VLM, forming the foundation for a better understanding of the VLM's role in vital functions and motor control.
Significance statement
The ventrolateral medulla (VLM) is an anatomically complex region of the brain that plays a crucial role in regulating vital functions, including autonomic and respiratory control, sleep–wake behaviors, cranial motor functions, and locomotion. This study comprehensively classifies VLM cell types and neuronal subtypes based on their molecular and anatomical features, by leveraging single-nuclei RNA sequencing, RNA fluorescence in situ hybridization, and axonal tract tracing. We present a dataset comprising 114,805 single-nuclei transcriptomes that identify and validate the precise molecular characteristics of neurons involved in autonomic and motor systems functions. This publicly available dataset offers new opportunities for comprehensive experimental studies to dissect the central organization of vital homeostatic functions and body movement.
Introduction
The ventrolateral medulla (VLM) is a brain region that is critical for the regulation of autonomic function, respiratory control, sleep–wake behavior, cranial motor function, and locomotion (Guyenet, 2014; Scammell et al., 2017; Del Negro et al., 2018; Arber and Costa, 2022; Coverdell et al., 2022; Kleinfeld et al., 2023). In adult mammals, the VLM is characterized by interlacing networks of fiber bundles and lacks easily demarcated cytoarchitectural boundaries (Watson et al., 2019). However, neurophysiological and neuroanatomical studies have identified functionally specialized groups of neurons in the VLM associated with unique molecular profiles (Guyenet et al., 2013; Shi et al., 2017; Coverdell et al., 2022; Veerakumar et al., 2022; Yackle, 2023) and functionally specialized astrocytes (Kasymov et al., 2013). As these studies have shown, establishing the transcriptional landscape of the VLM is important for understanding the VLM's role in vital functions. However, our understanding of the molecular landscape of the VLM is still limited, posing a challenge for experimental studies. Single-nuclei RNA sequencing (snRNA-seq) parses cellular diversity based on the genome-wide expression profiles of individual cells and groups them by similarity into candidate cell types. Applying this approach to discrete brain regions in mice has revealed a remarkable diversity of cell types, raising hypotheses about their functional roles, and providing genetic access points for their further investigation (Habib et al., 2016; Lake et al., 2016; Bakken et al., 2018; Liu et al., 2023). Several studies have applied single-cell or single-nuclei RNA-seq to specific cell populations within the VLM (Hayes et al., 2017; Bachmutsky et al., 2020; Beine et al., 2022; Coverdell et al., 2022; Veerakumar et al., 2022; Winter et al., 2023). However, a complete molecular database of all VLM cell types does not yet exist. We therefore used high-throughput snRNA-seq to comprehensively and systematically characterize the molecular cell types of the mouse VLM. Furthermore, we used anatomically and genetically targeted cell sorting of spinally-projecting and adrenergic/noradrenergic VLM neurons to identify and characterize neurons that provide input to neurons in the spinal cord involved in sympathetic and locomotor function.
Materials and Methods
Animals
All experiments were conducted in accordance with the National Institutes of Health's Guide for the Care and Use of Laboratory Animals and approved by the University of Virginia Animal Care and Use Committee (protocol #4312). Animals were group-housed wherever possible and maintained under a 12 h light/dark cycle at 23–24°C with water and food provided ad libitum.
We obtained H2B-TRAP mice (Roh et al., 2017; RRID:IMSR_JAX:029899) as a gift from Evan Rosen and Linus Tsai (Beth Israel Deaconess Medical Center, Harvard Medical School). H2B-TRAP mice express mCherry-tagged nuclear H2B protein (H2B-mCherry) and an eGFP-tagged ribosomal L10 protein from the Rosa26 genomic locus upon removal of a translational stop cassette by Cre recombination. Dbh-Cre mice [strain STOCK Tg(Dbh-cre)KH212Gsat/Mmcd; RRID:MMRRC_032081-UCD] procured from the Mutant Mouse Regional Resource Centers were maintained as Cre hemizygous (DbhCre/+). DbhCre/+ mice were crossed with H2B-TRAP homozygous mice to generate DBHTRAP mice. Dbh+/+ Cre-negative mice hemizygous for H2B-TRAP were injected with AAVrg-Cre in the spinal cord to generate SpinalTRAP mice. For anatomical tract-tracing studies, we used R26-loxSTOPlox-L10-GFP Cre-reporter mice acquired from The Jackson Laboratory (strain #024750, RRID:IMSR_JAX:024750) and bred as heterozygotes with C57BL/6J mice in house.
Single-nuclei RNA-seq library generation
To avoid stress- and anesthesia-related changes in gene expression, we rapidly decapitated unanesthetized mice upon scruffing and extracted their brains. Each sample batch consisted of tissue from 5–8 mice of each sex (8–20 weeks of age). Each brain was placed in ice-cold Earle's balanced salt solution (catalog #E2888; Sigma-Aldrich) or Hibernate-A (catalog #A1247501; Thermo Fisher Scientific) for 2–5 min, transferred dorsal surface up to a chilled brain matrix (catalog #SA-2175; Roboz), and sectioned in the coronal plane at 1 mm intervals. Keeping the tissue chilled throughout, we microdissected RVLM tissue with a scalpel blade on an ice-filled petri dish. We then processed the brain tissue into cell nuclei suspensions using one of two protocols, as described below. Importantly, tissue for all sample batches was dissected using the same approach, though the postdissection processing of tissue samples differed between batches 1–2 and batches 3–6, as described below (details provided in Table 1).
Information on batch and sample characteristics
For sample batches 1 and 2, we transferred the dissected tissue to an RNase-free microcentrifuge tube, froze it in a slurry of dry ice and methanol, and then stored at −80°C for less than a month. After thawing a batch of tissue samples on ice, we isolated cell nuclei from the thawed tissue using a previously described hypotonic lysis protocol (Matson et al., 2018), modified as follows: trituration with 5 Pasteur pipettes, starting with an unpolished pipette, followed by a series of four pipettes with tips fire polished to progressively smaller diameters; 5–10 strokes for each size pipette, or until the tissue passed smoothly through the pipette; and filtered resuspended nuclei with 20 µm strainer (preseparation filter; catalog #130-101-812; Miltenyi Biotec). We fluorescently labeled nuclei with propidium iodide and then quantified them with a DeNovix CellDrop automated cell counter. We then proceeded with the Chromium Next GEM Single Cell 3′ v3.1 kit (catalog #PN-1000128 or PN-1000121; 10x Genomics) according to the manufacturer's instructions (CG000204 Rev D; 10x Genomics), tagging each sequencing library with a unique 10× Sample Specific Index (catalog #PN-1000213; 10x Genomics) for sequencing libraries in multiplex.
For sample batches 3 through 6, which involved fluorescence-activated cell sorting (FACS), we isolated cell nuclei from preserved RVLM tissue samples based on a published protocol (Habib et al., 2016; Todd et al., 2020). RVLM tissue was microdissected as described above but then transferred to RNase-free microcentrifuge tubes containing ice-cold RNAprotect reagent (catalog #76104; Qiagen). We stored the tissue samples in RNAprotect at 4°C overnight, then removed the liquid, and transferred the preserved tissue samples to −80°C storage for no more than 3 months. When ready to proceed, we thawed a batch of tissue samples on ice, homogenized them in 1 ml of homogenization buffer 25 times with pestle A and 25 times with pestle B (Kimble Kontes Dounce Tissue Grinder, DWK Life Sciences, catalog #885300-0001), and then separated the nuclei from debris by density gradient centrifugation according to the protocol (Habib et al., 2016). We filtered the resuspended nuclei through a 20 µm cell strainer (Miltenyi Biotec preseparation filter), fluorescently labeled DNA with NucBlue Live ReadyProbes Reagent (catalog #R37605, Thermo Fisher Scientific), and sorted the nuclei on a Becton Dickinson Influx cell sorter equipped with an 85 µm nozzle and set to purity mode. All sorting was performed with a Specialist in Cytometry (SCYM; ASCP) certified technical assistance at the University of Virginia Flow Cytometry Core. To separate nuclei from non-nucleated debris, we gated for events with high relative intensity of NucBlue fluorescence. We then gated the nuclei by forward scatter area versus side scatter area to eliminate large aggregates, followed by forward scatter area versus forward scatter height to select for singlets. Finally, from the single nuclei, we selected H2B-mCherry+ and H2B-mCherry-negative nuclei based on their high and low relative fluorescence, respectively, with 561 nm excitation and 610/20 nm collection filter. We separated H2B-mCherry+ and H2B-mCherry-negative nuclei singlets into low-bind microcentrifuge tubes, each containing 18.8 µl of RT Reagent B (from Chromium Next GEM Single Cell 3′ v3.1 kit). We added the remainder of the 10× master mix reagents to each tube of sorted nuclei, before proceeding with the kit according to the manufacturer's instructions (CG000204 Rev D; 10x Genomics), tagging each sequencing library with a unique 10× Sample Specific Index (catalog #PN-1000213; 10x Genomics) for sequencing libraries in multiplex.
Single-nuclei RNA-seq data generation, processing, and analysis
We sequenced the single-nuclei RNA-seq libraries on an Illumina NextSeq 550 with high output, 75-cycle kits (catalog #20024906; Illumina). After sequencing each single-nucleus library to a targeted minimum of 20,000 reads per nucleus, we demultiplexed the sequencing reads, aligned them to the mouse genome mm10 with Cell Ranger v5.0 (10x Genomics) including intronic reads with the command–use introns, and quantified gene-level (feature) expression values for each cell-specific barcode as a feature-barcode matrix. We then processed each feature-barcode matrix with CellBender (v0.2.2; parameters can be accessed through Zenodo, record number: 10048289) to remove ambient RNA reads and random barcode swapping from the dataset, with the exception of batches 1 and 2 where the lack of a clear inflection point in the barcode rank plot undermined CellBender's performance.
We analyzed the data with the Seurat software package (v4.3.0.1) in R (v4.0.2) using default Seurat parameters unless otherwise specified. Notably, the H2B-mCherry+ cell nuclei from A5 noradrenergic neurons were excluded from the analysis because of low yield, preventing a robust analysis. Across all samples, we excluded cells whose transcriptomes contained >5% mitochondrial transcripts. To remove putative doublet samples, the upper quartile (Q3) and the interquartile range (IQR) of the number of genes and the number of UMIs in samples in each batch were calculated. Within each batch, the formula Q3 + 3 * IQR was used to calculate the upper cutoff values for the number of genes and the number of UMIs, and samples exceeding these values were excluded from further analysis. The lower cutoff for the number of genes was different per batch and was chosen based on each batch's distribution of this quality metric. For each batch, we log-normalized the data using Seurat's NormalizeData() function; selected the 2,000 most variable genes (feature selection), scaled expression of each gene, and used principal component analysis (PCA) for linear dimensionality reduction of the transcriptomes in high-variance gene space. Reciprocal principal component analysis-based (RPCA) integration was then used to correct for technical differences between batches. This integration created a new integrated assay, for which we then scaled the expression of each gene; used PCA for linear dimensionality reduction of the transcriptomes in a high variable gene space; clustered the cells using the Louvain algorithm, based on the Euclidean distance in the PCA space comprising the first 30 principal components (PCs) and with a resolution value of 0.5; and performed nonlinear dimensionality reduction by Uniform Manifold Approximation and Projection (UMAP) for visualization of the dataset in two dimensions.
To match each sample to a known cell type, we checked each cluster's expression for cell type marker genes (Table 2) and annotated them accordingly. We identified marker genes for each cluster using the Seurat FindAllMarkers function, returning only positive markers but with otherwise default parameters, which used a Wilcoxon rank-sum test and Bonferroni-corrected significance threshold to compare gene expression of cells within a given cluster against all other cells.
Cell type identity genes
We identified clusters representing cell doublets (i.e., transcriptome from more than one cell) based on their (1) coexpression of marker genes for multiple clusters and (2) otherwise lacking specific marker genes and excluded these doublet clusters from further analysis. We identified ambiguous clusters, likely representing cells with degraded or contaminated transcriptomes, based on our lack of detection of specific marker genes in these clusters. It is possible that these clusters expressed specific marker genes which we failed to detect due to gene dropout or other technical limitations of single-cell transcriptomics. After removing doublet and ambiguous clusters, we repeated the postintegration data processing with the same parameters. For each subsequent cell-type level analysis (astrocytes, all neurons, 5-HT neurons, and Cartpt/Dbh neurons), we repeated this pipeline of separately processing each batch, RPCA integration, followed by reprocessing the new integrated assay and identification of marker genes with the Seurat FindAllMarkers function. In these analyses, if doublet or ambiguous clusters were removed after the first postintegration analysis, a second round of this pipeline was performed until all doublet and ambiguous clusters were removed. This iterative workflow and parameters can be accessed through Zenodo, record number: 10048289.
To further analyze the neurons in our dataset, we subsetted the clusters labeled as neurons based on their expression of classical neuron marker genes. To ensure only high-quality neuron samples were included in this subsequent analysis, cells within these clusters that did not individually express Syn1, Rbfox3, or Syp were excluded. Neurons in which we detected expression of the oligodendrocyte marker gene Prr5l were suspected as contaminated or doublets and were also excluded.
Integration with Winter et al. (2023)
To compare our spinally-projecting neuron subtypes to those described in Winter et al., we analyzed the datasets with the Seurat software package (v4.3.0.1) in R (v4.0.2) using default Seurat parameters unless otherwise specified. We selected the clusters within our dataset that contained substantial amounts of SpinalTRAP samples (n06–n10, n16, and n20) and the clusters in the Winter et al. dataset that were annotated with a “MED” prefix. For each dataset, we log-normalized the data using Seurat's NormalizeData() function, selected the 2,000 most variable genes (feature selection), scaled expression of each gene, and used PCA for linear dimensionality reduction of the transcriptomes in high-variance gene space. RPCA-based integration was then used to correct for technical differences between the datasets. This integration created a new integrated assay, for which we then scaled the expression of each gene; used PCA for linear dimensionality reduction of the transcriptomes in a high-variance gene space; clustered the cells using the Louvain algorithm, based on the Euclidean distance in the PCA space comprising the first 30 principal components (PCs) and with a resolution value of 0.5; and performed nonlinear dimensionality reduction by UMAP for visualization of the dataset in two dimensions.
Intracardial perfusion for immunofluorescence and in situ hybridization
Mice (mixed sex, 8–12 weeks of age) were deeply anesthetized with a mixture of ketamine (100 mg/kg) and dexmedetomidine (0.2 mg/kg) given intraperitoneally and perfused transcardially with 4% paraformaldehyde (19208, Electron Microscopy Sciences), pH 7.4, in 100 mM phosphate buffer or pH-buffered 10% formalin, pH 7.0 (HT501128, Sigma-Aldrich). Brains were removed and postfixed in the same fixative for 12–24 h at 4°C. Brains were sectioned (30–50 µm) on a vibratome (VT-1000S, Leica Biosystems), and sections were stored in a cryoprotectant [30% ethylene glycol (v/v), 20% glycerol (v/v), 50% 100 mM phosphate buffer, pH 7.4] at −20°C. Immunohistochemistry was performed on free-floating sections at room temperature unless noted otherwise. Serial 1-in-3 sections were rinsed; then blocked in a solution containing 100 mM Tris, 150 mM saline, 10% horse serum (v/v) and 0.1% Triton X (v/v); and incubated with primary antibodies for 60 min at room temperature then 4°C overnight. Sections were subsequently rinsed and then incubated with secondary antibodies for 60 min at room temperature and rinsed again before mounting on slides. Slides were mounted with ProLong Gold antifade reagent with DAPI (P36931, Thermo Fisher Scientific).
Multiplex fluorescent in situ hybridization (FISH) was performed using RNAScope (V1 kit, Advanced Cell Diagnostics). Serial 1-in-3 sections were washed in RNase(ribonuclease)-free phosphate-buffered saline, mounted on charged slides, dried overnight, and processed according to the manufacturer's instructions. When necessary, immunohistochemistry was performed on the slide after the RNAScope procedure using antibodies against mCherry or GFP along with their respective secondary antibodies. Immediately following the RNAScope procedure, sections were rinsed and then incubated in blocking solution for 10 min followed by incubation in primary antibodies at room temperature for 60 min, rinsed and incubated in secondary antibodies for 30 min, and rinsed and then dried overnight before coverslipping. Slides were coverslipped with Prolong Gold antifade mountant with DAPI (P36931, Thermo Fisher Scientific). Primary and secondary antibodies and FISH probes used are detailed in Tables 3 and 4, respectively.
Antibodies used for immunohistochemistry
Probes used for fluorescent in situ hybridization
Microscopy, imaging, and cell mapping
Neuronal mapping was conducted using Neurolucida software (MBF Bioscience) with an AxioImager M2 microscope (Carl Zeiss Microscopy). Digital images were acquired in grayscale using a Hamamatsu C11440 Orca-Flash 4.0LT digital camera (Hamamatsu). Max projections of z-stack images were exported (8 bit), and further image modifications were performed in Fiji software. Representative images were pseudocolored and optimized for presentation, and brightness and contrast were adjusted equally in all pixels of the image. Cell counting was performed manually and unblinded using a 1-in-3 series that was subsequently aligned to Paxinos and Franklin's 4th edition (Paxinos and Franklin, 2012). In Figures 9 and 10, cell counts are reported for the rostral, caudal VLM, A2, and C2/C3 region. Rostral VLM (i.e., C1 region) was considered to be a region within 500 µm caudal of the facial nucleus, posterior to a bregma level −6.47 mm. Caudal VLM (i.e., A1 region) was considered to be a region within 500 µm rostral of the spinal decussation (bregma level: −7.55 to −8.05 mm). A2 neurons were counted in the ventrolateral nucleus of the solitary tract (bregma level: −7.55 to −7.67 mm); counts were based on neurons expressing intense Dbh or Slc6a2. The C2/C3 region was considered the region of the dorsal medial medulla (bregma level: −6.75 to −6.95 mm). To generate these counts of cells in the aforementioned regions, all cells expressing the gene/genes of interest were counted in 2 or 3 adjacent sections in a 1 in 3 series. Cells were considered to express a gene of interest if the cell contained at least five particles of fluorescent signal that were not in multiple fluorescent channels.
Results
Transcriptional analysis confirms major cell classes and uncovers astrocyte heterogeneity
To classify cell types of the VLM by their gene expression profiles, we dissected VLM tissue from adult male and female mice for single-nuclei RNA-sequencing (snRNA-seq). We used naïve C57Bl/6J mice and H2B-TRAP(mCherry) Cre-reporter mice which were either crossed with DBH-Cre mice, resulting in H2B-mCherry expression in VLM adrenergic/noradrenergic neurons (DBHTRAP mice; Fig. 1A–C), or injected in the thoracic spinal cord with a retrograde AAV-Cre, resulting in H2B-mCherry expression in spinally-projecting neurons (SpinalTRAP mice) (Fig. 1D–G). In DBHTRAP mice, nuclear H2B-mCherry expression was observed in 76 ± 8% of adrenergic/noradrenergic VLM neurons identified by expression of Dbh by RNA FISH (Fig. 1B). Furthermore, nuclear H2B-mCherry expression was observed in 84 ± 8% tyrosine hydroxylase (TH)-expressing VLM neurons, determined by immunofluorescence (Fig. 1C). Importantly, Dbh expression was detected in all cells that expressed H2B-mCherry (based on 313 H2B-mCherry+ cells in three cases) and 89 ± 3% of H2B-mCherry+ cells expressed TH. Hence, DBHTRAP mice provide a selective and efficient targeting approach to capture adrenergic/noradrenergic VLM neurons. In SpinalTRAP mice, we observed H2B-mCherry expression in neurons in the VLM and ventromedial medulla (VMM; Fig. 1E–G), including TH+ neurons in the rostral VLM, representing 36% of all TH+ cells between −6.59 and −7.07 mm of bregma (55 of 150 neurons in two cases; Fig. 1F,G). To control for the possibility that cells were transduced due to the leakage of AAV into the CSF, we pipetted 1 µl of AAV onto the exposed surface of the thoracic spinal cord followed by wound irrigation and closure. This control procedure did not produce any detectable H2B-mCherry expression in the brainstem.
Approach for collection of ventrolateral medulla (VLM) noradrenergic/adrenergic and spinally-projecting neurons for single-nuclei RNA-seq. A, Genetic approach to label Dbh-expressing neurons for collection, i.e., generation of DBHTRAP mice. Figure made with BioRender. B, Expression of H2B-mCherry in the VLM of DBHTRAP mice colocalizes with Dbh+ neurons by RNA fluorescence in situ hybridization (FISH). Yellow arrows identify double-labeled neurons; white arrows identify Dbh+ neurons that did not contain a H2B-mCherry-positive nuclei. Scale bars: 100 µm for the main panel and 50 µm for the bottom panels. C, Expression of H2B-mCherry in the VLM of DBHTRAP mice colocalizes with TH immunofluorescence. Yellow arrows identify double-labeled neurons; white arrows identify TH+ neurons that did not contain a H2B-mCherry-positive nuclei. Scale bar, 50 µm. D, Approach to label spinally-projecting neurons for collection, i.e., generation of SpinalTRAP mice. Figure made with BioRender. E, Spinal cord section from a SpinalTRAP mouse showing the location and approximate spread of AAV injections. Scale bar, 100 µm. F, SpinalTRAP cells in the VLM express TH. Bregma level: −6.64 mm. Yellow arrows identify double-labeled neurons; white arrows identify TH+ neurons not containing H2B-mCherry positive nuclei. Scale bar, 50 µm. G, Distribution of neurons expressing H2B-mCherry in the medulla of SpinalTRAP mice and the region of interest (ROI) for tissue collection for single-nuclei RNA-seq. H, Brain slice showing remaining tissue after bilateral ROI was collected for tissue processing.
To prepare samples for snRNA-seq of the VLM, we collected 1-mm-thick transverse sections of the medulla between −6 and −8 mm caudal from bregma. We then dissected bilateral “wedges” of the VLM bounded medially by the pyramids and laterally by the spinal trigeminal tract with the tip of the wedge in the intermediate reticular nucleus (Fig. 1H). Care was taken to avoid collecting midline raphe nuclei (raphe pallidus and obscurus) and the trigeminal region. Based on reference anatomical atlases, the cells collected for sequencing likely originated from regions including paragigantocellular reticular nucleus, parapyramidal nucleus, nucleus ambiguus, facial nucleus, lateral reticular nucleus, magnocellular reticular nucleus, gigantocellular reticular nucleus, intermediate reticular nucleus, and the inferior olive (Fig. 1H).
We isolated cell nuclei from the dissected VLM tissue and profiled their genome-wide mRNA content by snRNA-seq. Nuclei from DBHTRAP and SpinalTRAP tissue samples were sorted by flow cytometry to separately purify single mCherry-positive and mCherry-negative cell nuclei, whose poly-adenylated RNA were then processed in parallel into sequencing libraries. After sequencing the libraries and quantifying genome-wide expression values, we filtered out low-quality transcriptomes and putative cell doublets. This procedure left 114,805 single-nuclei transcriptomes (∼cells) for further analysis. Annotated Seurat code for the dataset can be accessed through Zenodo (record number: 10048289).
We integrated the cells across six sample batches, clustered them based on similarity of expression of high-variance genes, and annotated each cluster according to its expression of canonical cell class marker genes (Fig. 2A–F, Table 2). Clusters of neurons, identified by their expression of genes encoding NeuN (Rbfox3), Synaptophysin (Syp), and Synapsin 1 (Syn1), comprised the largest proportion of cells at 41.0%. A smaller cluster of neurons (Syn1, Rbfox3) distinguished by high expression Foxp2 likely originated from the inferior olive (Foxp2+ neurons; 1.8% of cells; Fig. 2G). Oligodendrocytes, marked by their expression of the oligodendrocyte transcription factor gene Olig1, made up the second largest cluster at 36.6% of all cells. The relative abundance of oligodendrocytes was expected given the high white matter content of the VLM. Astrocytes, identified by their expression of marker genes (Agt, Gfap, Slc13a3, Slc1a2), made up the third largest cell cluster (12.0% of cells). Interestingly, reclustering the astrocytes apart from the other cell types revealed five molecularly distinct subtypes (Fig. 2H). Polydendrocytes, also known as oligodendrocyte precursor cells (OPCs), were distinguished by their expression of the oligodendrocyte transcription factor genes Olig1 and Olig2 and the chondroitin sulfate proteoglycan 4 (NG2) gene Cspg4 and represented 3.3% of all cells. Approximately 2.3% of the cells formed a macrophage cluster, containing mostly microglia (P2ry12+) and a relatively smaller number of perivascular macrophages (Mrc1+). The remainder of the dataset comprised two distinct clusters of vascular and leptomeningeal cells (Ranbp3l+; 1.0% of cells) and one cluster each of pericytes (Kcnj8+, Abcc9+; 0.5% of cells), endothelial cells (Slco1c1+; 1.3% of cells), and choroid plexus cells (CP; Kl+, Folr1+; 0.1% of cells). Overall, our molecular taxonomy of VLM cells includes all expected cell classes and reveals an unexpected degree of molecular heterogeneity among VLM astrocytes.
Single-nuclei RNA-seq of all cells in the VLM. A, UMAP of all cells, labeled according to cell type. B, Genes differentially expressed to a significant extent across all cell types. C, Cluster correlation matrix for all cells. D, Dendrogram, number of cells, unique molecular identifiers, genes detected, and the expression of cell type marker genes for cell types shown in B. E, UMAP of all cells colored by sample batch. F, Batch composition per cluster for all cells. G, Expression of Foxp2 from Allen Mouse Brain Atlas (https://mouse.brain-map.org/experiment/show/72079884). Note the intense expression in the inferior olive. H, UMAP of astrocyte subtypes. I, Dot plot of differentially expressed genes in astrocyte subtypes.
Subclustering neurons reveals 23 candidate molecular subtypes
To better resolve the diversity of VLM neurons, we subsetted all neurons except those putatively from the inferior olive and reclustered them apart from the non-neuronal cells. After removing potential doublets and molecularly ambiguous cells, we grouped 15,342 neurons into 23 clusters (Fig. 3A). On average (± standard deviation), we detected 2,277 ± 1,213 genes in each neuron, based on 4,806 ± 3,760 transcripts per neuron (Fig. 3B).
Diversity of neurons in the VLM. A, UMAP of all neurons, labeled by cluster ID. B, Dendrogram showing neuron cluster relatedness; numbers of cells per cluster, UMI per cell, and genes detected per cell for each of the neuron clusters shown in A. C, Dot plot of selection of differentially expressed genes for neuron clusters shown in A. D, Heat map of neurotransmitter phenotype related genes for neuron clusters shown in A.
Comparing gene expression across neuronal clusters revealed candidate marker genes for each cluster (Fig. 3C). To identify the anatomical origins of each neuronal cluster, we cross-referenced its top marker genes with publicly available ISH data from the Allen Brain Atlas (Lein et al., 2007). Table 5 provides a summary of the putative anatomical origin of each population based on this analysis. To predict the neurochemical phenotype of each neuron cluster, we assessed the expression of genes involved in neurotransmitter synthesis and cell identity (Fig. 3D): excitatory neurons (Slc17a7/VGluT1, Slc17a6/vGluT2), inhibitory neurons (Slc32a1/VGAT, Slc6a5/GlyT2, Gad1/GAD67, Gad2/GAD65), cholinergic neurons (Chat/ChAT, Slc18a3/VAChT), serotonergic neurons (Slc6a4/SERT1, Tph2/Tryptophan Hydroxylase 2, Fev/Pet-1), and adrenergic/noradrenergic neurons (Dbh/Dopamine Beta-Hydroxylase, Slc17a6/vGluT2 and Slc6a2/NET). This analysis identifies nine clusters of inhibitory neurons, seven clusters of excitatory neurons, two clusters of serotonergic neurons (n07, n08), one cluster of adrenergic/noradrenergic neurons (n09), and one cluster of cholinergic neurons (n21). Two clusters (n10, n23) expressed markers for both excitatory and inhibitory neurons but rarely in the same cells (35/1,001 cells for n10 and 6/232 cells for n23), suggesting these clusters are heterogeneous.
Presumptive location of neuron cluster based on differentially expressed genes and neurotransmitter phenotype
SpinalTRAP classifies seven molecular subtypes of bulbospinal neurons
To identify neurons in our snRNA-seq dataset that innervate the spinal cord, we tracked SpinalTRAP cells throughout the analysis and found that these cells distributed across several neuron clusters but were especially enriched in n06–n10, n16, and n20 (Fig. 4A–C). Of note, SpinalTRAP-positive cells clustered with SpinalTRAP-negative cells in clusters n06–n10 and n16, while n20 was predominantly composed of SpinalTRAP-positive cells. Among these were serotonergic neurons (n07 and n08) and adrenergic/noradrenergic neurons (n09), both of which are a source of innervation for sympathetic preganglionic neurons in the spinal cord (Bowker et al., 1981; Tucker et al., 1987). Clusters n06, n16, and n20 expressed novel marker transcripts that were detectable by ISH in the VLM and VMM (Fig. 4D). The excitatory cluster n20 was marked by C1ql2, a gene expressed at high levels in VMM neurons, more specifically the ventral part of the gigantocellular nucleus and medial portion of the lateral paragigantocellular nucleus. The location and morphology of C1ql2+ neurons (Fig. 4D), neurotransmitter phenotype (Fig. 3D), and transcription factors (Fig. 4F) of cells in this cluster suggest that n20 represents excitatory reticulospinal neurons involved in locomotion and limb control (Cepeda-Nieto et al., 2005; Bretzner and Brownstone, 2013; Bouvier et al., 2015; Capelli et al., 2017; Leiras et al., 2022; Winter et al., 2023). The other two populations of bulbospinal neurons (n06 and n16) in our dataset were inhibitory. Cluster n16 may represent a population of bulbospinal inhibitory neurons that are broadly distributed in the VLM based on the expression of marker genes (Lypd6b; Figs. 3C, 4E). This cluster also expressed GlyT2 (Slc6a5) and VGAT (Slc32a1; Fig. 3D), similar to neurons identified in rats with a putative sympathetic function (Schreihofer et al., 1999; Stornetta et al., 2004). The other inhibitory bulbospinal cluster was defined by high expression of Crhbp (Fig. 3C), which encodes a protein involved in the negative regulation of corticotropin-releasing hormone, CRH (Potter et al., 1991, 1992). We mapped the distribution of bulbospinal Crhbp neurons by combining axonal tracing from the spinal cord with FISH. Neurons expressing Crhbp were located in both the VLM and VMM, with Crhbp+ cells in the VMM tending to have larger cell bodies than those in the VLM (Fig. 4G,H). Bulbospinal Crhbp+ neurons were located primarily in the ventral part of the gigantocellular nucleus and medial lateral paragigantocellular nucleus (Fig. 4G), suggesting that that expression of Crhbp in the VMM may identify inhibitory neurons involved in behavioral arrest (Capelli et al., 2017).
SpinalTRAP identifies clusters of spinally-projecting neurons. A, Sample ID composition of neuron clusters. Note that clusters containing SpinalTRAP neurons were integrated in clusters with cells from other samples. B, All-neuron UMAP (from Fig. 2A) colored to show the distribution of SpinalTRAP neurons. C, Distribution of SpinalTRAP samples across neuron clusters, shown as a proportion of the total number of SpinalTRAP cells in the dataset. D, Expression of C1ql2 (https://mouse.brain-map.org/gene/show/86401) in the VLM from Allen Mouse Brain Atlas. E, Expression of Lypd6b (http://mouse.brain-map.org/experiment/show/77332684) in the VLM from Allen Mouse Brain Atlas. F, Dot plot of differentially expressed genes in n20 (C1ql2 and Shox2) and genes expressed by neurons in the VMM that are involved in motor control (Lhx4, Lhx3, and Vsx2/Chx10; see text for references). G, Expression of Crhbp by FISH in the VLM and overlap with spinally-projecting neurons (Spinal-GFP). Yellow arrows in the right panel indicate expression of Crhbp in spinally-projecting neurons. Note that the images in the left and right panel are from different tissues. Scale bars, 100 µm. H, Distribution of Crhbp neurons in the medulla identifying the location where Crhbp was found in spinally-projecting neurons. I, Heat map of top differentially expressed genes in SpinalTRAP -positive and -negative samples. J, Gene ontology analysis of genes differentially expressed between SpinalTRAP-positive and SpinalTRAP-negative samples.
To compare spinally-projecting neurons and non-spinally-projecting neurons in our dataset, we first identified genes which were differentially expressed between SpinalTRAP-positive and SpinalTRAP-negative neurons and then categorized these genes by gene ontology. Our differential expression analysis identifies 423 and 401 genes significantly enriched in SpinalTRAP-positive and SpinalTRAP-negative neurons, respectively (Fig. 4I). Gene ontology analysis (Fig. 4J) of the differentially expressed genes revealed enrichment of transcripts controlling cytoplasmic translation and nucleoside triphosphate metabolic process in SpinalTRAP-positive neurons, consistent with differences in translational capacity and metabolism. On the other hand, SpinalTRAP-negative neurons were enriched with transcripts controlling gamma-aminobutyric acid signaling pathway and adult behavior. Interestingly, both SpinalTRAP-positive and SpinalTRAP-negative neurons were enriched with expression of genes in the synapse organization GO category, though the gene identities differed between groups. This suggests that the profile of synapse organization genes differs significantly between centrally versus spinally-projecting neurons.
Molecular identification of VLM respiratory neurons
The VLM contains neuronal circuits that underlie respiratory and orofacial motor activity patterns (Smith, 2022; Kleinfeld et al., 2023). We reasoned that these networks might be represented in our dataset as one or more discrete neuron clusters. We examined the distribution of genes known to be expressed by respiratory-related neurons (Tacr1, Sst, Oprm1, Penk; Gray et al., 1999; Stornetta et al., 2003b; Stornetta, 2008; Cui et al., 2016; Bachmutsky et al., 2020), which were enriched in several clusters in our snRNA-seq dataset but overlapped consistently in clusters n05 and n12 (Fig. 5A). Interestingly, Sst, which is expressed in both inhibitory and excitatory respiratory-related neurons (Stornetta et al., 2003a; Tupal et al., 2014; Sherman et al., 2015; Cui et al., 2016), was enriched in n05, a cluster of inhibitory neurons, and expressed at relatively high levels in a small subset (∼5%) of neurons in n12, a cluster of excitatory neurons. Furthermore, Cdh9, a putative marker of excitatory respiratory-related neurons in mice (Yackle et al., 2017; Vagnozzi et al., 2022), was enriched in excitatory clusters n12, n16, and n18 (Fig. 5A). Reelin, which is expressed in respiratory-related neurons (Tan et al., 2012), was also enriched in n05, n12, and along with seven other clusters. Pax2, which is reportedly expressed by excitatory respiratory-related neurons (Wu et al., 2017; Xia et al., 2022), was enriched in excitatory cluster n12 supporting the possibility that this cluster contains excitatory respiratory-related neurons. Finally, neurons regulating diaphragm activity are reported to be enriched with Cdh6 and Cdh10, in addition to Cdh9, in embryonic mouse medulla (Vagnozzi et al., 2022). In our snRNA-seq dataset based on adult mice, Cdh6 and Cdh10 expression was enriched in multiple clusters, suggesting that these genes do not effectively discriminate between respiratory neurons and other neurons in the VLM in adult mice. Altogether, the distribution of genes associated with respiratory-related neurons in our snRNA-seq dataset indicates clusters n05 and n12 are most likely to contain neurons involved in generating respiratory motor activity. Consistent with this, genes enriched in clusters n05 and n12 such as Gpc3 and Glra3 were localized to the VLM and VMM (Fig. 5B,C). Crhbp, which was enriched in n05, is expressed by primarily non-bulbospinal neurons in respiratory-related regions of the VLM (Fig. 4F). Also, both n05 and n12 had relatively high expression of Trpc4, which was also evident in cells localized to the respiratory column (Fig. 5D). Of note, a relatively small number of SpinalTRAP neurons was present in both n05 and n12, whereas n16 had a relatively large number of SpinalTRAP neurons (Fig. 4B,C). However, n16 expresses markers for inhibitory neurotransmission (Fig. 5A), which is unlike phrenic premotor neurons (Wu et al., 2017). Finally, we examined our dataset for genes expressed by central chemoreceptor neurons in the retrotrapezoid nucleus (RTN), which regulate breathing in relation to brain pH (Guyenet et al., 2019). In mice, there are 600–700 RTN neurons, which can be distinguished from neighboring neurons by the expression of Nmb and Gal (Shi et al., 2017; Souza et al., 2023). In our dataset, Nmb and Gal were detected in only a few cells and only codetected in four cells across the entire dataset. This shows that RTN neurons do not form a unique cluster in our dataset, which may be attributed to these cells being a relatively small population of cells in the VLM.
Identification of novel marker genes for VLM respiratory neurons. A, Violin plot of marker genes enriched for respiratory neurons (see text for references) and for neurotransmitter phenotypes, and differentially expressed genes for neuron cluster n05 and n12. Refer to Figure 3 for a more comprehensive list of differentially expressed genes for neuron cluster n05 and n12. B, Expression of Gpc3 (http://mouse.brain-map.org/experiment/show/71020431) in the VLM from Allen Mouse Brain Atlas. C, Glra3 (https://mouse.brain-map.org/experiment/show/73788474) in the VLM from Allen Mouse Brain Atlas. D, Trpc4 (http://mouse.brain-map.org/experiment/show/1306) in the VLM from Allen Mouse Brain Atlas.
Molecular subtyping of serotonergic neurons
Previous single-cell transcriptomic analysis has revealed molecular diversity among serotonin neurons in the distinct raphe serotonergic cell groups (Spaethling et al., 2014; Okaty et al., 2015; Huang et al., 2019; Ren et al., 2019). To learn whether similar heterogeneity exists among the parapayramidal 5-HT neurons, we subclustered the 5-HT neurons in our snRNA-seq dataset. First, we selected neuron clusters expressing the phenotypic markers of 5-HT cells, Slc6a4 and Tph2, which were restricted to n07 and n08 (Fig. 6A). Fev (i.e., Pet1) was expressed at high levels in <10% of cells in n07 and n08. Importantly, Slc6a4, Tph2, and Fev expression was restricted to neurons in clusters n07 and n08. Moreover, these clusters differed from each other in their expression of genes including Trh, Cpne7, Cpa6, and Gm43154 (Fig. 6A). Transcripts for marker genes of n07 and n08 were detected in both the parapyramidal region and midline raphe (pallidus and obscurus; Fig. 6B,C) and both n07 and n08 had a similar proportion of SpinalTRAP neurons. This suggests that n07 and n08 represent molecularly distinct but anatomically intermixed subpopulations of 5-HT neurons. Subclustering of 1,226 cells from n07 and n08 resulted in six distinct subgroups of 5-HT neurons (Fig. 6F,G) which differed in their expression of many novel marker genes (Fig. 6D), as well as genes that differentiate 5-HT neurons elsewhere in the brain (Fig. 6G; Okaty et al., 2015; Huang et al., 2019; Xiao et al., 2021; Schneeberger et al., 2022). Cell identity genes (Tph2, Slc6a4) were present in all 5-HT subgroups, although their expression appeared to be highest in cells in n07/08.s01, n07/08.s05, and n07/08.s06. SpinalTRAP neurons were common in 5-HT subclusters n07/08.s02, n07/08.s03, and n07/08.s04 (Fig. 6E, Table 6). A comparison of feature plots show SpinalTRAP cells had variable expression of Tph2 (Fig. 6E,F). Collectively, this analysis indicates that the parapyramidal 5-HT cell group comprises distinct populations of spinally-projecting neurons defined by a combination of newly identified marker genes as well as genes that discriminate 5-HT neurons in other brain regions.
Identification of subtypes of serotonergic neurons in the parapyramidal region. A, Violin plot of phenotypic genes for serotonergic neurons (Slc6a4, Tph2) and expression of differentially expressed genes for clusters n07 and n08 in “all-neuron” clusters. B, Expression of Cpne7 in VLM from Allen Mouse Brain Atlas (http://mouse.brain-map.org/experiment/show/73817432). C, Expression of Cpa6 in VLM from Allen Mouse Brain Atlas (http://mouse.brain-map.org/experiment/show/75197606). D, UMAP generated by subclustering n07 and n08 based on differentially expressed genes. E, UMAP from panel D indicating the distribution of SpinalTRAP samples and pie chart showing the proportion of 5-HT SpinalTRAP samples in each subcluster. Values for pie chart provided in Table 6. F, Feature plot showing the expression of Tph2 across 5-HT subclusters. G, Expression of phenotypic genes and genes that are expressed in 5-HT subpopulations in other brain regions (left panel) and differentially expressed genes for 5-HT neuron subclusters (right panel).
Values for pie charts in Figure 8
Molecular subtyping of noradrenergic and adrenergic neurons
Adrenergic/noradrenergic neurons in the ventrolateral medulla regulate many important aspects of physiology, including the sympathetic control of blood pressure, glucose homeostasis, immune function, and the neuroendocrine responses to inflammation and stress (Guyenet et al., 2013). Our snRNA-seq analysis identified cluster n09 as adrenergic/noradrenergic neurons based on expression of Cartpt, Dbh, Th, and Maoa (Fig. 3C). Both DBHTRAP and SpinalTRAP cells were prevalent in n09 and a majority of cells in cluster n09 expressed Dbh (Fig. 7A). Conversely, Pnmt, a marker for adrenergic neurons, and Slc6a2, the norepinephrine transporter, were expressed in different regions of the n09 feature plot (Fig. 7A). This observation aligns with the established expression patterns of Pnmt and Slc6a2 expression in adrenergic/noradrenergic VLM neurons in rats (Guyenet et al., 2013). Therefore, we subclustered n09 to further characterize the heterogeneity of adrenergic/noradrenergic neurons. Our analysis of 894 cells resulted in nine molecularly distinct subclusters (Fig. 7B). By tracking the cells derived from DBHTRAP and SpinalTRAP samples, we identified six of nine subclusters containing cells that were genetically tagged for Dbh expression (Fig. 7C, Table 7), whereas six of nine subclusters contained SpinalTRAP cells (Fig. 7D). Both Cartpt and Dbh were present in most subclusters (Fig. 7E–G). Consistent with this result, Dbh was detected by FISH in 52 ± 4% of all Cartpt+ VLM neurons, whereas Cartpt was detected in 81 ± 2% of Dbh+ VLM neurons (Fig. 8A,B). Of note, Slc17a6 (VGluT2) and Slc6a2 (NET) transcripts were present in nonoverlapping Cartpt/Dbh subclusters (Fig. 7G), and subclusters expressing Slc6a2 (n09.s06 and n09.s04) expressed Dbh and Th at relatively high levels (Fig. 7G). Of note, SpinalTRAP neurons were present in n09.s04, but not in n09.s06 (Fig. 7D). Collectively, this indicates that there are potentially two subtypes of noradrenergic neurons in the VLM.
Identification of subtypes of VLM Cartpt/Dbh neurons. A, All-neuron UMAP colored to show the distribution of SpinalTRAP neurons (far left). Insets show the distribution of SpinalTRAP and DBHTRAP samples the expression of phenotypic genes expressed by noradrenergic/adrenergic neurons in cluster n09. B, UMAP generated by subclustering n09 neurons based on differentially expressed genes. C, UMAP from panel B colored to show the distribution of DBHTRAP neurons and a pie chart showing the proportion of DBHTRAP neurons in Cartpt/Dbh subclusters. Values for pie chart provided in Table 7. D, UMAP from panel B colored to show the distribution of SpinalTRAP samples and a pie chart showing the proportion of SpinalTRAP samples in Cartpt/Dbh subclusters. Values for pie chart provided in Table 7. E, Feature plot showing the distribution of cells expressing Cartpt in Cartpt/Dbh subclusters. F, Feature plot showing the distribution of cells expressing Dbh in Cartpt/Dbh subclusters. G, Expression of phenotypic genes and genes that are reported to be expressed in VLM Dbh neurons (left panel) and differentially expressed genes for Cartpt/Dbh subclusters.
Distribution and colocalization of Cartpt and Dbh in the VLM of mouse. A, Expression of Dbh and Cartpt in the VLM based on FISH. Scale bar, 100 µm. B, Co-expression analysis of Dbh and Cartpt in the VLM based on FISH.
Values for pie charts in Figure 7
Presympathetic C1 neurons are known to express the gene Slc17a6 (Stornetta et al., 2002), which encodes the vesicular glutamate transporter 2. We detected Slc17a6 expression in the subclusters that expressed Cartpt at high levels, as well as one subcluster (n09.s03) that did not express Dbh or Cartpt. Pnmt was present in cells in several subclusters, and we observed an enrichment of Dbh and Th along with Slc17a6 in 5/9 subclusters (Fig. 7G) consistent with descriptions of VLM adrenergic neurons (i.e., the C1 neurons). Furthermore, the coexpression of these genes (Dbh, Th, and Slc17a6) across multiple Cartpt/Dbh subclusters suggests the presence of several molecularly distinct subpopulations of C1 neurons. Several neuropeptides that are known to be expressed in C1/A1 neurons in rats were enriched in Dbh+ subpopulations in our dataset; Npy was expressed in three clusters, one of which expressed high Slc6a2. Npy is enriched in C1/A1 neurons in the rat and rabbit (Blessing et al., 1987; Tseng et al., 1993; Stornetta et al., 1999). Npy expression was greatest in n09.s09, which contains a substantial number of SpinalTRAP neurons, and expressed at lower levels in n09.s08 and n09.s06. We also found expression of the gene encoding pituitary adenylate cyclase-activating polypeptide, PACAP (Adcyap1), in subcluster n09.s08, which is expressed by spinally-projecting catecholaminergic neurons in the rostral VLM in rats (Farnham et al., 2008). Consistent with this, n09.s08 was an excitatory cluster of Dbh+ neurons that included SpinalTRAP neurons. As such, the molecular profile of the mouse Cartpt/Dbh neuron subtypes identified by our dataset conforms with expectations based on anatomical studies in other species.
Importantly, our analysis also identified a number of novel genes that differentiate the subclusters of Cartpt/Dbh neurons (Fig. 7G, right). To provide validation for this analysis, we examined the expression of EYA transcriptional coactivator and phosphatase 1 (Eya1), which was a marker gene for subcluster n09.06. We first assessed Eya1 expression across all neurons in our dataset (Fig. 9A), finding only sparse expression including within the Cartpt/Dbh cluster. This indicates that Eya1 is expressed by a limited number of neurons in the VLM of the adult mouse. Despite this, Eya1 expression was highly selective for the Cartpt/Dbh subcluster n09.s06 (Fig. 9B), which expressed high levels of Dbh and Slc6a2, and no Slc17a6 (Fig. 7G). Interestingly, the expression of Eya1 was strongly correlated with Slc6a2 at a single cell level, indicative of coexpression (Fig. 9C). Based on this, we predicted that we would find a high degree of colocalization between the expression of Eya1, Dbh, and Slc6a2, which we evaluated by FISH (n = 2 mice). Indeed, Eya1 colocalized with Dbh and Slc6a2 consistently in the caudal VLM, aligning with the location of the A1 neurons, whereas Eya1 signal was diffuse in the rostral VLM (Fig. 9D–G). In the caudal VLM (i.e., the A1 region), in which virtually all Dbh neurons expressed detectable Slc6a2, we observed Eya1 expression in 76% of Dbh+/Slc6a2+ neurons (Fig. 9D,E,G,H). Conversely, very few Dbh+ neurons in the rostral VLM (i.e., the C1 region) expressed Slc6a2 (12/130 neurons) or Eya1 (2/130 neurons; Fig. 9F–H). Interestingly, Eya1 was also coexpressed with Dbh and Slc6a2 in the A2 and A6 region, but not in the C2/C3 region, in which Slc6a2 expression was also rare (Fig. 9G,H,I). These observations show that Eya1 is enriched in brainstem noradrenergic neurons and provides a putative marker for the A1 neurons.
Eya1 is enriched in VLM A1 neurons. A, Violin plot of Eya1 expression in “all-neuron” clusters. B, UMAP from Figure 7B (Dbh/Cartpt subclusters) colored to show expression of Eya1. C, UMAP from Figure 7B (Dbh/Cartpt subclusters) colored to show codetection of transcripts for Eya1 and Slc6a2. Codetection of these transcripts was primarily concentrated in subcluster n09.s06. D, Expression of Dbh, Slc6a2, and Eya1 in the caudal VLM (A1 region) by FISH. Yellow arrows in the lower panel indicate cells coexpressing Dbh, Slc6a2, and Eya1. Scale bar, 100 µm. E, High-magnification image demonstrating colocalization of Dbh, Slc6a2, and Eya1 in the caudal VLM (A1 region). Scale bar, 100 µm. F, Expression of Dbh, Slc6a2, and Eya1 in the rostral VLM (C1 region) by FISH. Scale bar, 100 µm. G, Distribution and overlap of Dbh, Slc6a2, and Eya1 expression in the VLM. H, Cell counts for the expression of Dbh, Slc6a2, and Eya1 in medullary catecholaminergic cell groups based on two separate cases. I, Expression of Dbh, Slc6a2, and Eya1 in the locus ceruleus (A6) by FISH. Scale bar, 100 µm.
To further validate our results, we examined the distribution of Hexokinase 2 (Hk2), a gene that encodes a kinase that phosphorylates glucose to regulate glucose metabolism, that is present in Dbh+ neurons based on whole mouse brain single-cell RNA-seq datasets (Langlieb et al., 2023; Yao et al., 2023). In our snRNA-seq dataset, high expression of Hk2 was largely restricted to the Cartpt/Dbh cluster (Fig. 10A) and was enriched in 4 Cartpt/Dbh neuron subclusters (n09.s01, n09.s02, n09.s07, n09.s08; Fig. 10B), of which three contained SpinalTRAP neurons (Figs. 7D and 6). Furthermore, Cartpt/Dbh subclusters in which Hk2 expression was detected exhibited low Slc6a2 and relatively high Slc17a6 (Fig. 7G), indicative of a glutamatergic rather than noradrenergic neurotransmitter phenotype. This information suggests that Hk2 is expressed by subpopulations of Cartpt/Dbh VLM neurons, including neurons that innervate the spinal cord. We tested this prediction by mapping Hk2 in the VLM combined with staining for Dbh and retrograde tracing from the spinal cord. Based on FISH (n = 2 mice), Hk2 was expressed by 60% of Dbh+ neurons in the rostral VLM (132/221 neurons) and only 1/138 Dbh+ neurons in the caudal VLM (Fig. 10D–F). Moreover, 59% of spinally-projecting Dbh+ neurons in the rostral VLM expressed detectable Hk2 (45/76 neurons; Fig. 10D–F). Notably, we also observed Hk2 expression in both spinally-projecting and non-spinally projecting Dbh+ neurons in the C2/C3 cell group regions (Fig. 10D,F), whereas Hk2 was not detectable in A1 or A2 (Fig. 10D,F) and was not observed in A5 or A6 neurons (data not shown) regardless of whether these neurons projected to the spinal cord or not. This led us to hypothesize that Hk2 was enriched in adrenergic neurons, which are defined by the expression of Pnmt (Guyenet et al., 2013). Based on FISH, Pnmt was detected in 56% of Hk+ neurons in the rostral VLM, whereas 78% of Pnmt+ neurons also expressed Hk2 (Fig. 10G,H). Hk2 also colocalized with Pnmt in the C2/C3 region, with 78% of Hk2+ neurons expressing detectable Pnmt and 98% of Pnmt neurons expressing Hk2 (Fig. 10G,H). Hk2 and Pnmt were rarely expressed by noradrenergic cells in the A1 and A2 region based on the expression of Slc6a2 (Fig. 10G,H). This collectively shows Hk2 expression is enriched in adrenergic neurons as judged by conventional criteria irrespective of the cell's location in the medulla.
Hk2 is enriched in VLM adrenergic neurons. A, Violin plot of Hk2 expression in “all-neuron” clusters. B, UMAP from Figure 7B (Dbh/Cartpt subclusters) colored to show Hk2 expression. C, Locations of images in D and G. D, Expression of Dbh and Hk2 by FISH, and detection of spinally-projecting neurons (Spinal-GFP) by immunofluorescence in the rostral VLM/C1 region (larger image and left panels), C2/C3 region (middle panels), and the caudal VLM/A1 region (right panels). Yellow arrows indicate cells coexpressing Dbh, Hk2, and Spinal-GFP; white arrows indicate cells expressing Hk2 and Dbh but not spinal-GFP. Scale bars, 50 µm for all. E, Distribution and overlap of Dbh, Hk2, and spinal-GFP in the VLM. F, Cell counts for the expression of Dbh, Hk2, and spinal-GFP in medullary catecholaminergic cell groups based on two separate cases. G, Expression of Hk2, Pnmt, and Slc6a2 by FISH in the rostral VLM/C1 region (larger image and left panels), C2/C3 region (middle panels), and the caudal VLM/A1 region (right panels). Yellow arrows indicate cells coexpressing Pnmt, Hk2, and Slc6a2; white arrows indicate cells expressing Slc6a2 only. Scale bars, 50 µm for all. H, Cell counts for the expression of Dbh, Hk2, and Slc6a2 in medullary catecholaminergic cell groups based on two separate cases.
Comparison of VLM dataset with Winter et al. (2023)
Single-cell RNA sequencing of spinally-projecting neurons in the medulla has characterized 37 distinct subpopulations (Winter et al., 2023). In addition to genes that are well established to label spinally-projecting neurons, we noted overlap in the genes identified in spinally-projecting neurons in the dataset of Winter et al. those in the VLM, such as Cpa6, C1ql2. Given this observation, we sought to characterize the convergence of cell identity between spinally-projecting clusters identified in the present with those identified by the dataset of Winter et al. To achieve this, we partitioned and integrated VLM clusters containing spinally-projecting neurons (total cells, 3,901) with medullary spinally-projecting clusters (i.e., cluster with the MED prefix) from Winter et al. (total cells: 14,894), clustering cells in the integrated dataset based on their expression of high-variance genes (Fig. 11A–F). This analysis resulted in 32 clusters (denoted with a prefix of i00 through i31; Fig. 11B,C), of which eight clusters had >25% of total cells originating from the VLM dataset after discounting cells from cluster n10, which was spread across 25 distinct clusters (Fig. 11D,E). VLM noradrenergic/adrenergic cells (n09.) separated into two clusters (i7 and i27) solely with cells originating from the MED-Lmx1b-Pnmt subpopulation identified by Winter et al. Both i7 and i27 expressed many of the same genes (Fig. 11F,G), although only i27 was enriched for Ttr (Fig. 11D) and Eya1 (19.7% of cells in i27 had detectable expression vs 2.3% of cells in other clusters; p = 7.30 × 10−44). VLM 5-HT clusters (n07, n08) grouped into four clusters (i9, i10, i12, i30) with cells originating from four subpopulations of neurons identified by Winter et al. (i.e., MED-Lmx1b-Tph2-Trh, MED-Lmx1b-Gad2-Cpa6, MED-Lmx1b-Tph2-Gad2-Cpa6, MED-Lmx1b-Gad2). The 5-HT clusters in the integrated dataset were defined by the expression of Tph2 and Slc6a4 (Fig. 11F) and distinguished from each other by differential expression of genes including Cpa6 (i9), Trh (i10), and Calcr (i11 and i30; Fig. 11F,G). Cells from VLM cluster n20, which was characterized by enrichment of C1ql2, grouped in cluster i14 predominantly with cells from MED-Lhx3/4-Kit-C1ql2, with a small number of cells from VLM cluster n20 grouping with other Lhx3/4 expressing subpopulations in cluster i4 and i13. Finally, cells from VLM cluster n16 grouped in cluster i11 with a majority of cells originating from the MED-Lhx1/5-Glyc-Ntsr1 subpopulation identified by Winter et al., while VLM cluster n06 grouped in cluster i23 with a small number of cells from the MED-Lhx1/5-Glyc subpopulation. Collectively, this data shows that the bulbospinal cell types defined by our VLM dataset correspond to those subpopulations identified by Winter et al. but that our dataset provides an improved resolution for VLM noradrenergic/adrenergic populations owing to the enrichment strategy used for cell collection. The presence of clusters that contain few cells from the dataset of Winter et al. could reflect subpopulations of bulbospinal neurons that are rare and/or restricted to the VLM or alternatively the presence of non-bulbospinal cells in the VLM dataset. The latter interpretation is supported by the enrichment of Eya1 in cluster i27.
VLM spinally-projecting clusters map to established identities for spinally-projecting neurons in the medulla. A, Procedure for partitioning and integration of VLM spinally-projecting neurons with a previously published dataset for spinally-projecting neurons in the medulla (Winter et al., 2023). B, UMAP for integrated dataset clustered according to high-variance genes. C, UMAP colorized based dataset of origin to identify the overlap of datasets. D, Number of neurons from VLM spinally-projecting dataset within the cell clusters in the integrated dataset. E, River plot showing the correspondence between the original spinally-projecting VLM neuron clusters and cluster in the integrated dataset. F, Patterns of gene expression across integrated dataset clusters for high-variance genes identified by clustering all neurons in the VLM dataset (Fig. 3). G, Expression of top marker genes for each cell cluster in the integrated dataset.
Discussion
The study provides a detailed transcriptional atlas of VLM cell types and a dissection of the major neuronal subtypes that innervate the spinal cord. Collectively, we characterized the nuclear transcriptomes of 114,805 cells and determined specific molecular markers that can be used to identify distinct subtypes of neuronal and non-neuronal cells, infer their function, and gain genetic access to them for further investigation. Among our findings, we identify several distinct populations of spinally-projecting neurons in the VLM that are potentially involved in the control of autonomic function and motor coordination. Together our results provide new insight into molecular and cellular organization of the VLM.
SpinalTRAP identifies the molecular features of spinally-projecting neurons in the VLM
A goal of this study was to characterize neuronal subpopulations in the VLM that innervate the spinal cord. To achieve this we profiled the transcriptome of spinally-projecting neurons labeled with injections in the thoracic spinal cord of AAV that is efficiently internalized by axons (AAVrg-Cre; Tervo et al., 2016). An important consideration in interpreting our results is that it was not possible to determine the effective spread of AAV injections in the spinal cord. Moreover, damage to axons en route to the lumbar and sacral cord may lead to transduction of neurons projecting outside of the upper thoracic segments of the spinal cord. Given these factors, it was not possible to determine the segmental targeting of SpinalTRAP neurons, which will require precise projection field mapping of each cell type. Furthermore, viral tropism could affect which cells are transduced by AAVrg, which could skew the proportions of each spinally-projecting population in on our dataset. Nevertheless, our analysis identifies novel spinally-projecting clusters in the VLM as well as provides the first detailed transcriptome for bulbospinal Cartpt/Dbh and 5-HT neuron types, which were expected be present in the SpinalTRAP populations (Bowker et al., 1981; Tucker et al., 1987). This was despite reports that AAVrg is inefficient at transducing serotonergic and noradrenergic cell types (Tervo et al., 2016; Ganley et al., 2021). This discrepancy may be explained by the higher efficiency of the Cre-reporter expression systems compared with promoter-based expression systems when expressing H2B-mCherry or a similar reporter. The novel clusters of spinally-projecting neurons in our dataset included one excitatory cluster (n20) and two inhibitory clusters (n06 and n16). Based on the integration of our dataset with previous work (Winter et al., 2023), cluster n20 reflects medullary Lhx3/4-expressing excitatory reticulospinal cells that overlap with Chx10 (Vsx2) neurons, which are therefore likely to be involved in voluntary movement, and involuntary postural and gait control (Arber and Costa, 2022; Leiras et al., 2022). Both cluster n06 and n16 appear to reflect two related subpopulations of Lhx1/5-expressing inhibitory reticulospinal neurons. Currently, the function of these cells is unknown, and we expect the genes identified as possible markers for these subpopulations could be useful to study these cells in the future.
Molecular characterization of adrenergic/noradrenergic neurons in the VLM
Monoaminergic neurons in brainstem were characterized in the 1960s and 1970 s by Dalhstrom Fuxe and others (Dahlstrom and Fuxe, 1964; Fuxe et al., 1970; Hökfelt et al., 1974). Since then, the expression of monoaminergic cell markers for noradrenergic, adrenergic, and serotonergic neurons in the VLM has provided an essential neuroanatomical handle for exploring the integrative function of this brain region. Catecholaminergic neurons in the VLM are a highly conserved group of cells classified as C1 or A1 based primarily on the expression of Pnmt and their relative location in the medulla (Blessing et al., 1978; Howe et al., 1980; Reiner and Vincent, 1986; Minson et al., 1990; Halliday and McLachlan, 1991; Ruggiero et al., 1992; Dormer et al., 1993). In addition, experimental studies indicate that subpopulations of C1 and A1 cells exist based on efferent connectivity (Tucker et al., 1987; Stornetta et al., 1999) and gene expression (Phillips et al., 2001; Pilowsky and Goodchild, 2002; Guyenet et al., 2013). Our data indicate there are seven molecularly distinct subpopulations of Dbh neurons in the VLM, including five adrenergic subtypes and two noradrenergic subtypes. We histologically verified two novel markers of the adrenergic and noradrenergic neurons in the VLM, one generated by our own data (Eya1) and the other (Hk2) identified in other studies (Langlieb et al., 2023; Yao et al., 2023). This study indicates that Eya1 is a unique marker for a subpopulation of noradrenergic neurons in the caudal VLM and also other noradrenergic subpopulations in the dorsal medulla and pons, but not adrenergic neurons. In contrast, Hk2 was expressed in four excitatory Cartpt/Dbh subclusters and neither of the noradrenergic subclusters, suggesting this gene was expressed primarily by glutamatergic/adrenergic cells. Hk2 expression was found primarily in Dbh+ neurons in the rostral VLM, including both bulbospinal and non-bulbospinal neurons, as well as majority of Pnmt+ neurons in the rostral VLM and C2/C3 regions. Hence, Hk2 expression appears to be a novel and fairly selective positive selection criterion for adrenergic neurons in the medulla. In addition to validating our snRNA-seq dataset, our identification of Eya1 and Hk2 reaffirms the notion that molecularly distinct adrenergic and noradrenergic neurons are spatially intermixed in the VLM, emphasizing the need for precise molecular and genetic approaches to dissect the function of each subtype.
SpinalTRAP identifies putative marker genes for noncatecholaminergic presympathetic neurons in the rostral VLM
In rats, adrenergic neurons represent a majority of neurons in the rostral VLM that project to sympathetic preganglionic neurons in the spinal cord, accounting for up to 70% of the total presympathetic neurons in this region of the brain (Lipski et al., 1995; Schreihofer and Guyenet, 1997). These cells have a well-established role in the sympathetic control of blood pressure (Souza et al., 2022). The remainder of the presympathetic population in the VLM in rats are non-catecholaminergic excitatory neurons, although a distinguishing marker for the “non-C1 neurons” remains elusive. In our dataset, neurons expressing adrenergic/noradrenergic genes (Dbh, Th) also expressed the neuropeptide gene Cartpt, similar to what has been found in rats (Koylu et al., 1998; Dun et al., 2002; Burman et al., 2004). Moreover, Burman et al. (2004) suggested that Cartpt identifies presympathetic VLM neurons involved in cardiovascular function. Interestingly, by subclustering the Cartpt/Dbh cells, we identified a subtype of glutamatergic neurons that express Cartpt and are largely devoid of catecholamine cell markers (neuron subcluster n09.s01). Many cells in this subcluster were SpinalTRAP cells, and therefore, this cluster could represent non-catecholaminergic presympathetic neurons in the VLM. Another candidate for non-catecholaminergic presympathetic neurons is subcluster n09.s03 as it contained a substantial number of SpinalTRAP cells and expressed Penk (Stornetta et al., 2001), but no detectable Cartpt or Dbh. Future genetically targeted functional experiments using the markers for VLM Cartpt/Dbh neurons identified in this study can be used to elucidate which neurons are critical for the sympathetic control of cardiovascular function.
Molecular characterization of serotonergic neurons in the parapyramidal region
Serotonergic neurons in the parapyramidal region are implicated in thermoregulation through sympathetic control of cutaneous blood flow, brown adipose tissue thermogenesis, cardiac and sudomotor function (McAllen, 1986; Dampney and McAllen, 1988; Bamshad et al., 1999; Pelaez et al., 2002; Cano et al., 2003; Nakamura et al., 2004), as well as respiratory functions (Ribas-Salgueiro et al., 2005; da Silva et al., 2010; Brust et al., 2014). This study provides a detailed transcriptomic analysis of 5-HT neurons in the parapyramidal region, a region that has been overlooked by previous studies. We discovered genes that distinguished two broad classes of 5-HT neurons, as well as genes that distinguished six putative subtypes within these classes. Furthermore, genes that are expressed in distinct types of 5-HT cells in other brain regions were also evident in different subtypes of 5-HT neurons in our dataset (e.g., Tacr1, Gad1/2, Calcr, Adra1a, Trh). This suggests that parapyramidal 5-HT neurons are composed of molecularly distinct subpopulations, similar to 5-HT groups in other brain regions.
Identification of putative markers for excitatory and inhibitory respiratory-related neurons
The VLM contains the neural network that generates breathing (Feldman et al., 2003; Del Negro et al., 2018). The classification of VLM neurons that generate the rhythm and pattern of breathing is principally based on their firing properties in relation to the respiratory cycle (Richter and Smith, 2014; Smith, 2022). However, studies have identified genes that are expressed by respiratory-related neurons with varying degrees of selectivity. We leveraged this information to identify putative respiratory-related clusters in our dataset (n05 and n12). These clusters appear to represent two classes of neurons distinguished from each other by their neurotransmitter phenotype (i.e., excitatory and inhibitory) and from other cells in the dataset by the differential expression of specific genes. The expression of the highly variable genes of n05 and n12 in the brainstem indicate these clusters are unlikely to directly relate to a specific functional compartment within the ventral respiratory column, such as the pre-Bötzinger complex. Future studies are needed to confirm that the genes identified in cluster n05 and n12 label respiratory-related neurons. Interestingly, this study shows that Trpc4, a nonselective calcium-permeable cation channel regulated by Gq and Gi/o G-protein-coupled receptor signaling (Schaefer et al., 2000; Thakur et al., 2016; Tian et al., 2022), is expressed at relatively high levels in both n05 and n12 compared with other neurons in our dataset. Trpc4 expression has been reported in functionally identified auto-rhythmic respiratory neurons in the pre-Bötzinger complex of neonatal mice (Kallurkar et al., 2022), providing additional support for our conclusion that n05 and n12 contains respiratory-related neurons. Given the abundance of Trpc4 expression in n05 and n12, we propose that the activity of TRPC4 in VLM respiratory neurons integrates the coactivation of Gq and Giio signaling and as such, contributes to the state-dependent control of respiratory motor output (Ramirez et al., 2012).
Conclusions and implications for future studies
Our study highlights the diversity of cell types in the VLM, many of which overlap in anatomical distribution, and provides a valuable resource for future studies involving the VLM. Leveraging this data, which can be explored using the Broad Institute Single Cell Portal, will facilitate a rigorous delineation of cells involved in the regulation of autonomic, respiratory, and motor function. This dataset also provides a snapshot of the transcriptional state of the VLM cells in healthy adult mice that can be referenced to identify and characterize signaling pathways that may be affected in diseases that involve changes in VLM function, such as neurodegenerative disease (Benarroch, 2019) and neurogenic hypertension (Guyenet et al., 2020), or through comparisons with genome-wide association studies. Overall, these data provide a basis for generating novel hypotheses concerning VLM function and offer the means to test them.
Footnotes
We thank Carla Winter and Zhigang He for sharing the clustered and annotated dataset from their Winter et al. (2023) published study. This work was supported by National Institutes of Health (NIH) R01 HL148004 to S.B.G.A. and R01 HL153916, Pathway to Stop Diabetes Initiator Award of the American Diabetes Association 1-18-INI-14, and UVA 3Cavaliers Grant to J.N.C.
↵*D.C.S., D.S.S., and R.-J.A.-F. contributed equally as first authors.
↵**J.N.C. and S.B.G.A. contributed equally as senior authors.
The authors declare no competing financial interests.
- Correspondence should be addressed to Stephen B. G. Abbott at sba6t{at}virginia.edu or John N. Campbell at jnc4e{at}virginia.edu.
