Transcriptional profiling of dental sensory and proprioceptive trigeminal neurons using single-cell RNA sequencing

Abstract

Dental primary afferent (DPA) neurons and proprioceptive mesencephalic trigeminal nucleus (MTN) neurons, located in the trigeminal ganglion and the brainstem, respectively, are essential for controlling masticatory functions. Despite extensive transcriptomic studies on various somatosensory neurons, there is still a lack of knowledge about the molecular identities of these populations due to technical challenges in their circuit-validated isolation. Here, we employed high-depth single-cell RNA sequencing (scRNA-seq) in combination with retrograde tracing in mice to identify intrinsic transcriptional features of DPA and MTN neurons. Our transcriptome analysis revealed five major types of DPA neurons with cell type-specific gene enrichment, some of which exhibit unique mechano-nociceptive properties capable of transmitting nociception in response to innocuous mechanical stimuli in the teeth. Furthermore, we discovered cellular heterogeneity within MTN neurons that potentially contribute to their responsiveness to mechanical stretch in the masseter muscle spindles. Additionally, DPA and MTN neurons represented sensory compartments with distinct molecular profiles characterized by various ion channels, receptors, neuropeptides, and mechanoreceptors. Together, our study provides new biological insights regarding the highly specialized mechanosensory functions of DPA and MTN neurons in pain and proprioception.

Fig. 1

Single-cell RNA sequencing (scRNA-seq) identified dental primary afferent (DPA) neuron subtypes. a Schematic workflow of scRNA-seq for DPA neurons. b A t-distributed stochastic neighbor embedding (t-SNE) plot showing six distinct clusters, C4, C6, C7, C8, C10, and C13, identified in a total of 83 DPA neurons from 10 mice. The cluster labels are referenced from a previously published scRNA-seq dataset of the TG. Separated clusters are shown in different colors. c Violin plots showing expression of subtype-specific enriched genes [e.g., Trpm8, Cd34, S100b, Piezo2, Calca, Trpv1, Nppb, and Mrgprd, presented in TG] in DPA neurons. d Two new sub-clusters exhibiting in C6, separated by C6A or C6B. e, f Most enriched genes in C6A [e.g., Chodl, Pvalb, Gfra3, and Tmem233] and C6B [e.g., Gpr83, Ntrk3, Rgs9, and Scn11a] (for more details, refer to Table S1)

Fig. 2

Transcriptome analysis of DPA neurons. a Distribution of five major clusters of DPA neurons (excluding C13) through multiscale bootstrap hierarchical clustering analysis. The top node of the hierarchical tree is denoted by as a blue-colored dot. Numbers in red indicate the approximately unbiased (AU) p-value for the corresponding cluster. b Heatmap displaying the top 45 differentially expressed (DE) genes distinguishing C4/C6 from C8/C10. Gene of interest is highlighted in red (refer to Table S2). c–h Heatmaps depicting genes in categories relevant to intrinsic properties of somatosensory neurons, including mechanoreceptors/myelination, neurotropic factor receptors, pain, voltage-gated sodium channels, voltage-gated calcium channels, and taste receptors. Data represent the log transform of the mean DESeq2-normalized counts of transcripts in each cluster. Gene names are shown with both gene symbols (italicized, first letter uppercase) and protein symbols (not italicized, all letters uppercase). i Transcript levels of Ntrk1 across DPA clusters. j Summary table presenting reference genes and functional annotations of each DPA cluster. LTMRs low-threshold mechanoreceptors, PEP peptidergic, MNs mechanosensitive nociceptors

Fig. 3

P2X₃ expression is frequently observed in IB₄-negative DPA neurons. a–d Representative fluorescence images (20X magnification) and quantitative results showing DPA neurons labeled with Fluoro-Gold (FG, blue-dotted outline) and each marker gene, Mrgprd (C13) and P2rx3, shown in green after RNAscope assay. Arrowheads indicate DPA neurons expressing each marker gene. Scale bars: 25 μm. Bar graph represents the cell body size distribution (µm²) of DPA neurons expressing Mrgprd (n = 1 of 62 FG⁺ DPA neurons from n = 2 mice) or P2rx3 (n = 12 of 36 FG⁺ DPA neurons from n = 2 mice) as the average of two mice. e A representative image displaying an IB₄-negative DiI-labeled DPA neuron (asterisk). A scale bar: 20 μm. f The number of IB₄-positive (n = 3 of 96 neurons, green) or IB₄-negative (n = 93 of 96 neurons, gray) DPA neurons during the collection process. g Transcript levels of P2Y (G-protein coupled) [e.g., P2ry1, P2ry2, and P2ry14] and P2X (ligand-gated) [e.g., P2rx3, P2rx4, P2rx5, P2rx6, and P2rx7] purinergic receptor family genes in DPA neuron transcriptomes. h Representative immunofluorescence images of FG-labeled DPA neurons (blue) showing P2X₃ (red) and IB₄-binding (green) in TG sections. Nissl stain (white) was used for cell body identification. Arrowheads indicate P2X₃-expressing DPA neurons but do not bind to IB₄. Scale bars: 25 μm. i Proportion of 4 categories [e.g., IB₄⁻/P2X₃⁺, IB₄⁻/P2X₃⁻, IB₄⁺/P2X₃⁺, IB₄⁺/P2X₃⁻] among all FG-labeled DPA neurons. Numbers represent mean ± SD. n = 100 neurons from n = 4 mice

Fig. 4

scRNA-seq analysis reveals heterogeneity of mesencephalic trigeminal nucleus (MTN) neurons. a Schematic diagram of MTN neuron preparation for scRNA-seq. b Representative fluorescence images for a couple of MTN neurons, which were retrogradely labeled with DiI (red) from the masseter muscle, in the brainstem. A scale bar: 100 µm. Insets represent the high-magnification images showing DiI (red) and Nissl staining (white). Scale bars: 20 µm. n = 1 mouse (P21). c Representative immunofluorescence images showing Advillin positivity (green) in cultured DiI-labeled MTN neurons. A scale bar: 20 µm. n = 1 mouse (P21). d Heatmap representing the distribution of two major subgroups of MTN neurons: those with high levels of both Ntrk3 and Ntrk2 (designated as Ntrk2^high) and those with high levels of Ntrk3 but low levels of Ntrk2 (designated as Ntrk2^low). e–h Heatmaps showing genes categorized according to intrinsic properties, including neurotropic factors and receptors, proprioceptor/mechanoreceptors, glutamate release and receptors, and acetylcholine receptors. Data represent the log transform of mean DESeq2-normalized counts of transcripts in each group (refer to Table S3). i–k Volcano plots showing the DE genes with a minimum log2 (fold change; FC) of 0.5 and an adjusted P-value < 0.05 that are enriched in the DPA neurons (left side) or MTN neurons (right side). For analysis, gene lists from Gene Ontology (GO) categories including sensory perception (GO:0007600), response to stimuli (GO:0050896), and a custom gene list of mechanosensitive ion channels were used. Genes of interest within the top 20 DE genes indicated with text labels), ranked by log2 FC in each GO term, are highlighted in blue (refer to Table S4)

Fig. 5

Postnatal changes in MTN neurons. a The number of DE genes enriched in P21 (left side) or P28 (right side) of MTN neurons (refer to Table S5). The gene list used for this analysis is within GO categories, including nervous system development (GO: 0007399), programmed cell death (GO: 0012501), muscle structure development (GO: 0061061), synapse organization (GO: 0050808), and translation at synapse (GO: 0140241). b Cell body size (μm²) of MTN neurons during the collection procedure. P21: n = 66 neurons and P28: n = 48 neurons. Bar graph presents as the mean ± SD. Unpaired nonparametric test (Kolmogorov-Smirnov test) was performed, *P= 0.015 6. c Representative fluorescence images showing DiI-labeled MTN neurons (red) and Nissl stain (green) located in the rostral or caudal part of the brainstem according to mouse age. Scale bars: 100 µm. d Violin plots quantitatively depicting neuronal cell body size (μm²) by age. n = 3 mice (P21), n = 2 mice (P28), and n = 2 mice (5-week). One-way analysis of variance with Bonferroni post hoc test was conducted, *P= 0.031 4