Performing scRNA-seq in DPA neurons innervating the tooth pulp

We performed scRNA-seq on DiI-labeled DPA neurons manually collected from the culture of the TG using SMART-Seq protocol, as illustrated in Fig. 1a. To ensure the DPA neuron identity, we selected cell bodies of DPA neurons with DiI labeling in the TG, specifically innervating the maxillary first molars (Fig. 1a and Fig. S1a). Additionally, we confirmed that all DiI-labeled cells expressed Advillin, a marker for somatic sensory neurons20 (Fig. S1b).

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.8 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 TG8] 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)

We first assigned the uniquely mapped reads to the RefSeq annotated genes to generate gene counts (Fig. S1c), then normalized the dataset to mitigate differences in sequencing library depth.21,22 The transcriptomes of a total of 83 DPA neurons collected from 10 adult mice showed an average of 11 936 detected genes per neuron. We observed no significant cross-sample contamination (data not shown) and noted no noticeable correlation among the number of detected genes (Fig. S1d). The dataset also passed our quality control criteria (Fig. S1e, f) and underwent validation with reference genes, including housekeeping genes [e.g., Gapdh and Actb (β-actin)], sensory neuronal markers [e.g., Tubb3 (β-tubulin III), Uchl1 (PGP9.5), Avil (Advillin), and Rbfox3 (NeuN)],12,13 glial markers [Gfap (Glial fibrillary acidic protein), Kcnj10 (Kir4.1), and Gja1 (Connexin-43)],23 and motor neuronal markers [e.g., Chat (ChAT) and Neurog2 (Neurogenin-2)]24 (Fig. S1g).

Identification of DPA neuron subtypes

Since it is generally accepted that somatosensory neurons are a heterogeneous cell population that encompasses touch (or low-threshold mechanosensitive)- and pain-related neurons, we leveraged a previously published scRNA-seq dataset of the TG,8 of which DPA neurons are a subset. This TG dataset used droplet-based sequencing on 3 500 neurons to identify and characterize 13 distinct classes of cells within the TG. Using Seurat’s TransferAnchor function,25 we assigned cell subtype cluster labels to our DPA neuron datasets based on previously established cluster identities, predicted to have the following functional classifications: C1-C2 cool-responsive neurons, C3 C-low threshold mechanoreceptors (C-LTMRs), C4-C6 Aβ-LTMRs/Aδ-LTMRs/Aδ-nociceptors respectively, C7-C10 different subtypes of polymodal peptidergic (PEP) nociceptors, C12 itch-transmitting neurons, and C13 non-peptidergic (NP) nociceptors.8 The 2-dimensional visualizations of well-separated clusters with different colors are displayed in Fig. 1b.

Our DPA neurons were annotated to six of the 13 TG clusters: namely C4, C6, C7, C8, C10, and C13, as per the cluster numbers of the TG dataset (Fig. 1b), according to the highest Seurat prediction score (Fig. S2a). Of note, since only a single DPA neuron out of 83 was annotated to C13 (Fig. 1b), we could not show its density plot in Fig. S2a. This single neuron assigned to C13 exhibited especially high expression of Mrgprd, which is a mas-related G protein-coupled receptor channel implicated in noxious mechanical pain transduction (Fig. 1c and Fig. S1h). This gene Mrgprd was also suggested as a marker gene for C13, predicted to be NP mechanosensitive nociceptors.8 The other 82 DPA neurons were annotated to five of the 13 previously defined TG clusters, with C6 (putative Aδ-nociceptors) being particularly abundant (42.2%, n = 35 of 83 neurons, as shown in Fig. 1b). In addition, these five DPA neuron clusters consisted of neurons with distinct cell body sizes (C7/C8/C10, predictive of polymodal PEP nociceptors appearing smaller than C6, predictive of Aδ-nociceptors on average) that supported further classification (Fig. S1i).

To compare our clusters with the clusters reported for the TG, we investigated the transcript levels of representative marker genes of the 13 TG clusters8,9 (Fig. 1c and Fig. S1h). These markers included a menthol receptor, Trpm8 (C1-C2, cool-responsive); a cell surface antigen, Cd34 (C3, C-LTMRs); a calcium-binding protein, S100b and a mechanosensitive ion channel, Piezo2 (C4–C6, respectively Aβ-LTMRs/Aδ-LTMRs/Aδ-nociceptors); calcitonin gene-related peptide, Calca, a heat and capsaicin receptor, Trpv1, and a wasabi receptor, Trpa1 (C7-C10, polymodal PEP nociceptors); natriuretic polypeptide B, Nppb (C11, itch-related neuropeptide); and Mrgprd (C13, NP nociceptors). While Trpm8 and Cd34 were highly expressed in C7 and C6 neurons respectively, the expression patterns of other representative genes, such as S100b, Piezo2, Calca, Trpv1, Nppb, and Mrgprd, were comparable with those of TG neurons in general (Fig. 1c and Fig. S1h). Mrgpra3 (C12, itch related), another mas-related G protein-coupled receptor channel mediating chloroquine-evoked itch, was filtered out as it was expressed in fewer than three cells. Therefore, our findings indicate that retrograde tracer-labeled DPA neurons fall within a subset of all TG neurons but seem to exclude some cluster markers that are possibly more associated with oral sensory afferents from the gingiva, palate, or tongue, facial skin afferents, or other non-dental afferents.

We also observed that the C6 neurons, presumed to be Aδ-nociceptors,8,9 exhibited notable expression of both Piezo2 and Calca (Fig. 1c), which could contribute to rapid nociception elicited by low-threshold mechanical stimulation. However, given the relatively lower expression of Trpv1 in C6 compared to that in C7/C8/C10 (Fig. 1c), we inferred that C6 is more associated with Aβ-fibers than with Aδ-fibers. More specifically, the DPA neurons belonging to C6 were sub-divided into two new sub-clusters, namely C6A and C6B, based on distinct expression profiles (Fig. 1d). The C6A sub-cluster displayed significantly elevated expression of Chodl (Chondrolectin), Pvalb (Parvalbumin), Gfra3 (GDNF family receptor alpha 3), and Tmem233, which encodes a transmembrane protein involved in the regulation of a voltage-gated sodium channel26 (Fig. 1e and Table S1), which suggested the C6A sub-cluster are Aβ-fiber slowly adapting (SA) LTMRs.27 The C6B sub-cluster showed enrichment of Gpr83 (G protein-coupled receptor 83), Ntrk3 (TRKC), Rgs9 (Regulator of G protein signaling 9), and Scn11a (Nav1.9) (Fig. 1f and Table S1), suggesting its classification as Aβ-nociceptors.27 Together, these findings demonstrate that DPA neurons consist of five major types that correspond to C4, C6, C7, C8, and C10.

Subtype-specific gene signatures of DPA neurons

To assess the robustness of the assigned clusters, we conducted a multiscale bootstrap hierarchical clustering using Euclidean distance on a subset of the 3 000 most variable genes. Bootstrap analysis revealed that the top node of the hierarchical tree distinctly separated C6 from the other clusters, with C4 also forming a separate subtree, distinct from two subtrees containing mixtures of C7/C8 and C7/C10 at a bootstrap replication p-value level of 0.05 (Fig. 2a). Thus, while the bootstrap tree did not show significant separation of some of the finer cluster relationships, the separation of C4/C6 from C8/C10 was robustly reproduced. We identified many differentially expressed (DE) genes enriched in C4/C6, such as Pvalb, Kcna1 (Kv1.1), Scn1a (Nav1.1), Scn8a (Nav1.6), Ldhb (Lactate dehydrogenase B), Nefh (Neurofilament heavy chain), S100b, and Piezo2, which are all known for their high expression in myelinated A-fibers13 (Fig. 2b and Table S2). Conversely, pain-related DE genes, such as Trpa1, Gal (Galanin), Calcb (Calcitonin gene-related peptide 2), Trpv1, Tac1 (Substance P), Npy1r (Neuropeptide Y receptor type 1), Tlr4 (Toll-like receptor 4), Calca, and Scn11a, were enriched in C8/C10 (ref. 28) (Fig. 2b and Table S2).

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

To further investigate the subtype classification of DPA neurons, more detailed gene families in somatosensory neurons were examined using the log transform of the mean (DESeq2 normalized) transcript counts in each cluster (Fig. 2c-h). By profiling genes specifically expressed or enriched in each cluster, we assigned clusters of DPA neurons to putative functional subtypes. C4 was assigned to Aβ-fiber LTMRs, marked by high expression of Calb1 (Calbindin 1, 28 kDa),11,29Sfrp1 (Secreted frizzled-related protein 1),27 and Ntrk3 (Ref. 12,27,29) (Fig. 2c, d), along with low expression of Scn10a and Scn11a, specifically expressed in nociceptive C-fibers among the voltage-gated sodium channel subtypes27,28,29 (Fig. 2f). This indicates that C4 is not involved in nociceptive functions. Interestingly, with the exception of C4, all other DPA neuron clusters showed high expression of Ntrk1 (TRKA), which is associated with nerve growth factor (NGF) (Fig. 2d, i). NGF signaling plays an important role in determining the fate of sensory neurons that become nociceptors and increases mechanically activated currents specifically in PEP nociceptive Aδ- and C-fibers.30,31 Hence, this result suggests that most DPA neurons possess the potential to function as mechanosensitive nociceptors (MNs) under elevated NGF conditions. Since C6 neurons have high levels of canonical markers for myelinated neurons [e.g., Nefh, S100b, Ldhb, and Thy1 (Fig. 2c)] and Calca (Fig. 2e), but low expression of Trpv1 (Fig. 2e), we classified them as heat-insensitive Aβ-fiber MNs. This category could further be divided into SA LTMRs and nociceptors, based on their sub-clusters (Fig. 1d-f). C7 neurons similarly expressed myelination markers for myelinated neurons (Fig. 2c), but featured higher expression of nociceptive markers [e.g., Runx1, Trpa1, Trpv1, and Tac1] than C6, consistent with the previously established C7 cell type in the TG8 (Fig. 2d, e). Consequently, we classified them as intermediate/polymodal Aδ-fiber nociceptors. Finally, C10 and C8 exhibited features associated with polymodal C-fiber nociceptors. Given that C8 and C10 featured high levels of Gal and Tac1 (Fig. 2e), they might represent subclasses of C-fibers transmitting different pain modalities (e.g., nerve injury pain or nociceptive sensitization).12 Interestingly, only C8 showed pronounced expression of Chrna3 (Cholinergic receptor nicotinic alpha 3 subunit), a known marker gene for silent mechanosensitive nociceptors27,31 (Fig. 2e), suggesting that this cluster could play a critical role in mechanical nociception under high NGF conditions. With respect to mechanoreceptors, the DPA neurons had low expression of Piezo1 and Tmc1 (Transmembrane channel like 1) & 2, which are associated with itch32 and hearing33 respectively, but high expression of Asic3 (Acid-sensing ion channel subunit 3)34 and Piezo2 (ref. 35), consistent with previous studies9,36,37 (Fig. 2c). Several voltage-gated calcium channels were moderately expressed across all clusters (Fig. 2g). Taste receptors genes were lowly expressed and similarly not specific to any cluster (Fig. 2h). In summary, most DPA subtypes, excluding C4, can be characterized as potent PEP MNs that are affected by NGF signaling (Fig. 2i, j).

P2X3 expression in the PEP subtype of DPA neurons

Our FISH assay confirmed that several reference genes from the general TG population were also observed in DPA neurons retrogradely labeled with Fluoro-Gold (FG) (Fig. S3a–h and Fig. 3a–d). In line with nociceptive neurons, Tac1 and Trpv1 were predominantly expressed in small to medium-sized DPA neurons (Fig. S3a–d), whereas Calca and Piezo2 displayed high expression across DPA neurons of varying sizes (Fig. S3e–h). Moreover, we consistently observed a low proportion of Mrgprd-positive DPA neurons (Fig. 3a, b), in agreement with our sequencing results (Fig. 1b, c and Fig. S1h). Mrgprd-positive neurons, which bind to IB4, are required for response to noxious mechanical stimuli.38 Consistent with previous reports,9,36 we observed that the majority of collected DPA neurons did not bind FITC-conjugated IB4 (n = 3 positive of 96 neurons, Fig. 3e, f), supporting the notion that DPA neurons are primarily composed of PEP nociceptors (Fig. 2j). We also found that a relatively high proportion of DPA neurons exhibited expression of the extracellular ATP ionotropic receptor P2rx3 (P2X3), showing a bimodal distribution across medium to large sizes (Fig. 3c, d). Notably, since P2rx3 is generally present in NP nociceptors that are distinct from Trpv1-positive neurons,39 its relatively high expression in DPA neurons, which are Mrgprd- and IB4-nagative and Calca-positive PEP subtypes, is an atypical observation in DRG neurons. Our DPA neuron dataset also showed that P2rx3 expression predominated over other P2X (ligand-gated) and P2Y (G-protein coupled) purinergic receptors40 (Fig. 3g). Our immunohistochemical assays further demonstrated that the majority of P2X3-positive DPA neurons did not bind IB4-negative (Fig. 3h, i). These findings thus support the notion that ATP, acting through the P2X3 receptor, could be a crucial modulator of nociceptive tone in DPA neurons that could contribute to mechanical sensitivity in the context of injury or infection.41

Fig. 3

P2X3 expression is frequently observed in IB4-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 (µm2) 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 IB4-negative DiI-labeled DPA neuron (asterisk). A scale bar: 20 μm. f The number of IB4-positive (n = 3 of 96 neurons, green) or IB4-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 P2X3 (red) and IB4-binding (green) in TG sections. Nissl stain (white) was used for cell body identification. Arrowheads indicate P2X3-expressing DPA neurons but do not bind to IB4. Scale bars: 25 μm. i Proportion of 4 categories [e.g., IB4−/P2X3+, IB4−/P2X3−, IB4+/P2X3+, IB4+/P2X3−] among all FG-labeled DPA neurons. Numbers represent mean ± SD. n = 100 neurons from n = 4 mice

Molecular diversity of MTN neurons innervating the masseter muscle spindles

As with DPA neurons, we performed scRNA-seq with the same protocol, but we specifically used juvenile mice to enhance the viability of MTN neurons cultured from brainstem slices (Fig. 4a). For sample preparation, we utilized the retrograde tracer DiI to label cell bodies originating from the masseter muscle spindles in a total of 8 juvenile mice (P21: 5 mice and P28: 3 mice) (Fig. 4a, b) and manually collected DiI-labeled MTN neurons expressing Advillin20 (Fig. 4c). Similar to DPA neurons, after quality control, we detected an average of 11 364 genes per MTN neurons from 108 MTN neurons (Fig. S4a–d). The MTN dataset was validated with the expression of specific marker genes associated with proprioceptive neurons [e.g., Pou4f1 (BRN3A), Ntrk3, Etv1(ER81), Pvalb, and Whrn (Whirlin)]42,43,44 in addition to the same set of reference genes that were previously employed for the validation of the DPA dataset12,13,23,24 (Fig. S4e). We conducted dimensionality reduction (t-distributed stochastic neighbor embedding; t-SNE) using expression data for the top 3 000 most variable genes, generated by Seurat, in pooled MTN samples of all ages. There was little batch effect between sequencing runs (Fig. S4f).

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 Ntrk2high) and those with high levels of Ntrk3 but low levels of Ntrk2 (designated as Ntrk2low). 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)

Proprioceptive neurons innervating muscle spindles are typically classified into two major types: group Ia and group II.45 A third type of proprioceptive neuron, group Ib, innervates Golgi tendon organs.45 A recent scRNA-seq study focusing on proprioceptive neurons in limb skeletal muscle spindles identified three major marker genes, Lmcd1, Fxyd7, and Chad, representing group Ia, II, and Ib afferents, respectively.42 MTN neurons that innervate masticatory muscles are also known to comprise two main afferents, group Ia and II.7 Given these similarities, we investigated whether MTN neurons could be similarly classified based on these three marker genes. Our analysis confirmed the presence of comparable subtypes in MTN neurons, namely Lmcd1-enriched (similar to group Ia afferents, designated as Lmcd1high:Fxyd7low) and Fxyd7-enriched (similar to group II afferents, designated as Lmcd1low:Fxyd7high) subtypes (Fig.