Somatosensory neuron types identified by high-coverage single-cell RNA-sequencing and functional heterogeneity

Abstract

Sensory neurons are distinguished by distinct signaling networks and receptive characteristics. Thus, sensory neuron types can be defined by linking transcriptome-based neuron typing with the sensory phenotypes. Here we classify somatosensory neurons of the mouse dorsal root ganglion (DRG) by high-coverage single-cell RNA-sequencing (10 950 ± 1 218 genes per neuron) and neuron size-based hierarchical clustering. Moreover, single DRG neurons responding to cutaneous stimuli are recorded using an in vivo whole-cell patch clamp technique and classified by neuron-type genetic markers. Small diameter DRG neurons are classified into one type of low-threshold mechanoreceptor and five types of mechanoheat nociceptors (MHNs). Each of the MHN types is further categorized into two subtypes. Large DRG neurons are categorized into four types, including neurexophilin 1-expressing MHNs and mechanical nociceptors (MNs) expressing BAI1-associated protein 2-like 1 (Baiap2l1). Mechanoreceptors expressing trafficking protein particle complex 3-like and Baiap2l1-marked MNs are subdivided into two subtypes each. These results provide a new system for cataloging somatosensory neurons and their transcriptome databases.

Figure 1

Neuron sampling, RNA-sequencing and gene clustering. (A) Triple-immunofluorescent staining showing three major neuron subsets labeled by IB4, CGRP and NF200 in mouse lumbar DRG. Scale bar, 50 μm (left) and 10 μm (right). (B) A schematic of the workflow depicting rapid dissociation and isolation of individual DRG neurons for single-cell RNA-seq profiling. Images show a single neuron positive or negative for IB4 in the tip of a glass pipette. Scale bar, 20 μm. (C) Distribution of mapped reads in 197 neurons. (D) Number of detected genes in 197 neurons. (E) In samples obtained from individual DRG neurons, the number of genes detected at a given sequencing depth (in this case, 52.7 million mapped reads) was correlates highly with that obtained at a triple depth (181 million mapped reads). (F) Transcript expression levels (x and y axes: log2-scale) in two samples obtained by equally dividing cDNA from DRG neuron No. 72 were correlate highly to each other. (G) The cluster dendrogram of 2 043 differential genes represents the co-expression modules identified by WGCNA. Modules of highly interconnected groups of genes correspond to branches and are labeled in different colors. (H) ISH shows the distribution of the representative module genes including Tmem176b, Nxph1, Il31ra, Baiap2l1, Cpne6, S100b, Tmem45b and Fam19a4 of the red, purple, yellow, green, blue, turquoise and black modules, respectively, in lumbar DRG. Scale bar, 50 μm.

Figure 2

Clustering of DRG neurons. (A) A heatmap of size-based hierarchical clustering shows the color-coded correlation matrix of 197 neurons based on the mRNA profiles of 1 745 differentially expressed genes (> 5 fold change). Ten clusters of strongly correlated neurons are marked with black frames and shown with the given names and their averaged area of neurons to the left. (B) The WGCNA heatmap shows the expression profiles of module genes. The eigengenes of neuronal clusters are indicated by yellow frames.

Figure 3

DRG neuron subclusters and the “hybrid states” of subclusters. (A) An enlargement of C4 on the heatmap of the correlation matrix and WGCNA eigengenes (Figure 2A and 2B) reveal two distinct subclusters, C4-1 and C4-2, outlined by black frames. (B) Single-cell real-time PCR confirming the differential expression of Mrgpra3, Mrgprb4 and Mrgprd in C4-1 and C4-2 neurons (n = 3). (C) Double fluorescent ISH showing the co-expression of Mrgprb4 in a subpopulation of Mrgpra3-positive small DRG neurons (arrows). Scale bars, 100 μm (left) and 20 μm (right). (D) Differential distribution of Mrgpra3, Mrgprb4 and Mrgprd in C4 and C5. (E) Gene co-expression network identified by WGCNA shows that the Mrgpra3-containing network is separate from the Mrgprb4 network in the pink module. (F) An enlargement of C2-1 and C2-2 in the heatmap of the correlation matrix and WGCNA (Figure 2A and 2B). (G) Single-cell PCR confirming the differential expression of representative genes in C2-1 and C2-2 neurons (n = 3). (H) Single-cell PCR showing expression of Pvalb in an Nppb-positive, IB4-negative neuron. (I) Double fluorescent ISH showing co-expression of Il31ra with S100b or Cpne6 in a small DRG neuron (arrow). Scale bar, 20 μm. (J) The correlation among Nppb, Il31ra, S100b and Cpne6 expression in clusters of small and large neurons and their predicted relationships with different types of afferent fibers.

Figure 4

Evaluation of neuron clusters markers. (A) Single-cell real-time PCR confirms the expression of marker genes in C1 (n = 6), C2 (n = 6), C3 (n = 6), C4 (n = 3) and C5 (n = 6) neurons. (B) Single-cell PCR reveals the expression of marker genes in C1 (n = 3), C2 (n = 6), C3 (n = 3), C4 (n = 3) and C5 (n = 6) neurons. (C) Single-cell PCR reveals the expression of marker genes in C5 (n = 3), C8 (n = 8) and C9 (n = 4) neurons. (D) A heatmap shows the expression patterns of cluster markers. New markers are indicated in blue. (E) A heatmap shows the expression patterns of Calca, Tac1 and Nefh, and the neurotrophin receptors Ntrk and Gfra in the neuron clusters. IB4 labeling is marked with a green square. (F, G) Double fluorescent ISH shows the cluster-specific expression of markers in DRG neurons. Il31ra is co-expressed with Nppb but not Gal in small neurons. Scale bar, 50 μm. (H, I) Th is co-expressed with both Fam19a4 and Fxyd6 but not Gal and Nppb in small neurons. (J) Mrgpra3 is not co-expressed with Nppb and Gfra2 in small neurons. (K) Mrgprd is not co-expressed with Il31ra in small neurons. (L) S100b and Th are not co-expressed in DRG neurons. (M) Pvalb-positive large neurons and Th-positive small neurons do not contain Ntrk2.

Figure 5

Functional characteristics of neuron clusters. (A) A diagram illustrating the potential correlation between neuron clusters and the differential expression of functional molecules involved in thermoception, mechanoception, nociception, chemicoception and pruriception. Representative functions of molecules are mainly based on the studies conducted in gene knockout mice. T, thermal nociception; M, mechanical nociception; F, formalin nociception. (B) A diagram illustrating the procedures of in vivo electrophysiological recording and simultaneous stimulus application to the plantar skin of the mouse hindpaw as well as neuron classification by single-cell PCR. (C) Examples of neurons responding to various cutaneous stimuli and their cluster identities. The application of a peripheral stimulus results in a large inward Na⁺ current and a small outward K⁺ current, and generates an action potential in the responding neuron. (D) Table showing the number of responsive neurons and their classification.

Figure 6

Neuron types and correlated somatosensory receptors. (A) The morphological characteristics of DRG neurons. (B) Classification of the types and subtypes of DRG neurons, their markers and the type hierarchy of DRG neurons. New types, subtypes and markers are indicated in red. (C) Traditionally classified subsets and their markers (upper block) are not suitable for the classification of neuron types in DRG (lower block). (D) The diagram shows the functional probability of the neuron types as somatosensory receptors. (E) A schematic showing the proposed framework of DRG neuron types and their proportions based on both transcriptomic analysis and ISH, as well as functional annotations suggested by electrophysiological analysis results and published data. C8 and C10 are predicted to be mechanoreceptors or/and proprioceptors.