1. Neuroscience
Download icon

Single-nucleus transcriptomic analysis of human dorsal root ganglion neurons

  1. Minh Q Nguyen
  2. Lars J von Buchholtz
  3. Ashlie N Reker
  4. Nicholas JP Ryba  Is a corresponding author
  5. Steve Davidson  Is a corresponding author
  1. National Institute of Dental and Craniofacial Research, National Institutes of Health, United States
  2. Department of Anesthesiology, College of Medicine, University of Cincinnati, United States
Tools and Resources
  • Cited 0
  • Views 1,075
  • Annotations
Cite this article as: eLife 2021;10:e71752 doi: 10.7554/eLife.71752

Abstract

Somatosensory neurons with cell bodies in the dorsal root ganglia (DRG) project to the skin, muscles, bones, and viscera to detect touch and temperature as well as to mediate proprioception and many types of interoception. In addition, the somatosensory system conveys the clinically relevant noxious sensations of pain and itch. Here, we used single nuclear transcriptomics to characterize transcriptomic classes of human DRG neurons that detect these diverse types of stimuli. Notably, multiple types of human DRG neurons have transcriptomic features that resemble their mouse counterparts although expression of genes considered important for sensory function often differed between species. More unexpectedly, we identified several transcriptomic classes with no clear equivalent in the other species. This dataset should serve as a valuable resource for the community, for example as means of focusing translational efforts on molecules with conserved expression across species.

Editor's evaluation

The manuscript by Nguyen et al. describes the assignment of neuronal cell types to human lumbar dorsal root ganglion (DRG) neurons based on the sequencing of individual nuclei. Bioinformatic comparison of these data to single nucleus sequencing results previously reported for mouse lumbar DRG is also described. The findings begin to close a gap in our understanding of how neuronal cell types and the expression of key genes differs between human and mouse. This kind of information is critical for targeted therapeutic efforts.

https://doi.org/10.7554/eLife.71752.sa0

Introduction

The somatosensory system responds to a wide range of mechanical, thermal, and chemical stimuli to provide animals with critical information about their environment and internal state. For example, our sense of touch is mediated by mechanosensory neurons with soma located in the dorsal root and trigeminal ganglia that innervate the skin (Abraira and Ginty, 2013). In addition to the skin, somatosensory neurons target specialized sensory environments like the cornea and conjunctiva or meninges (von Buchholtz et al., 2020; Huang et al., 2018), the internal organs (Hockley et al., 2019) as well as bones and muscles to provide rich perceptual experiences and trigger appropriate behavioral, reflex, and autonomic responses (Gatto et al., 2019). Among their many roles, somatosensory neurons provide input for the conscious perception of pain and itch (Basbaum et al., 2009; Mishra and Hoon, 2013) and the subconscious coordination of muscles and limbs known as proprioception (Chesler et al., 2016). Peripheral neurites of somatosensory receptor cells must adapt to growth, reinnervate targets after injury and are also affected by inflammation (Pongratz and Straub, 2013).

Studies in model organisms have characterized a range of sensory and growth factor receptors and ion channels that contribute to the properties and selectivity of somatosensory neurons (Le Pichon and Chesler, 2014; Coste et al., 2010; Ranade et al., 2014). Some of these, like the cooling and menthol sensing receptor (Trpm8) appear to define functional classes of cells (Bautista et al., 2007). By contrast, the sense of touch appears to use a complex distributed code involving several different types of cells (von Buchholtz et al., 2021b) to achieve its remarkable discriminatory power. For the most part, the human somatosensory system expresses the same range of functional genes as rodents (Ray et al., 2018) and exhibits similar responses to many types of stimulus (Adriaensen et al., 1983; Basbaum et al., 2009; Chesler et al., 2016; Le Pichon and Chesler, 2014; Ranade et al., 2014; Bautista et al., 2014; Davidson et al., 2016). Moreover, rare individuals with loss of function variants of several of these genes have deficits that recapitulate key effects of knocking out that gene in mice (Chesler et al., 2016; Murthy et al., 2018; Szczot et al., 2018; Drenth and Waxman, 2007; Chen et al., 2015). However, despite the identified similarities between mice and humans, the success of translating new therapeutic strategies that are effective for treating pain in mice has often been disappointing when tested in human subjects (Mogil, 2019; Yezierski and Hansson, 2018).

Recently, various directed genetic strategies have been used in mice to characterize the response properties and anatomical features of a variety of interesting classes of large diameter, fast conducting Aβ and Aδ subtypes (Abraira and Ginty, 2013). Interestingly, these neurons generally have complex peripheral endings that often target hair follicles. Human skin hairs are quite different from those in mice, suggesting that there may be significant differences between the large diameter neurons in mice and humans. By contrast, most types of small diameter, slow conducting c-fibers terminate as free nerve endings both in mice and humans (Basbaum et al., 2009). Single-cell sequencing approaches have produced a transcriptomic classification for mouse somatosensory neurons that corresponds well with their anatomy and function (Gatto et al., 2019; von Buchholtz et al., 2021b; Sharma et al., 2020; Nguyen et al., 2019). In mice, at least two classes of small diameter neurons are best defined by different members of the Mrgpr family of GPCRs (Sharma et al., 2020; Nguyen et al., 2019). However, Mrgprs have undergone massive genetic expansion in rodents, not seen in other animals, often making it difficult to identify true orthologs in humans (Dong et al., 2001; Liu et al., 2009) and raising questions as to whether similar cells would have distinct molecular markers in the two species. A map of human somatosensory neuron transcriptomic classes would help uncover selective differences between the sensory neurons in mice and humans and provide clues as to how similar somatosensory input is in the two species. Finally, such analysis may provide important new targets to consider for translational approaches to treat both pain and itch. Here, we used nuclei-based single-cell transcriptomics to generate an initial description of human cell types, highlight similarities and surprising differences between somatosensory neuron classes in humans and mice that are reflected not only in terms of individual genes but can be discerned in co-clustering. We used multigene in situ hybridization (ISH) to help confirm these conclusions and present evidence for anatomic organization of functionally distinct neuronal classes in the human dorsal root ganglion.

Results

Generating a representative transcriptomic map of human somatosensory cell types

Single lumbar L4 and L5 human dorsal root ganglia (DRG) were rapidly recovered from transplant donors within 90 min of cross-clamp and were immediately stored in RNAlater. Nuclei from individual ganglia were isolated and samples were enriched for neuronal nuclei by selection using an antibody to NeuN (see Figure 1—figure supplement 1). Five ganglia from one male and four female donors with ages ranging from 34 to 55 were subjected to droplet based single nucleus (sn) capture, barcoding, and reverse transcription (×10 Genomics). Combinatorial clustering methods (Stuart et al., 2019) allowed co-clustering of neuronal nuclei into a well-defined set of distinct transcriptomic groups that are well separated from their non-neuronal counterparts (Figure 1A, Figure 1—figure supplement 2). After removal of non-neuronal nuclei from the dataset, reclustering the DRG-neuron data from 1837 cells identified a range of about a dozen diverse transcriptomic classes of human somatosensory neurons (Figure 1B, Figure 1—figure supplements 2 and 3).

Figure 1 with 5 supplements see all
Diverse classes of human dorsal root ganglia (DRG) neurons revealed by single nuclear transcriptomics.

(A) Universal manifold (UMAP) representation of graph-based co-clustered snRNA sequences from human DRG nuclei reveal two well separated groups corresponding to sensory neurons (colored) and non-neuronal cells (gray). To the right, a dotplot highlights the expression of markers that help distinguish these groups of cells (see also Figure 1—figure supplement 2A for more information about preliminary analysis). (B) Reanalyzing 1837 neuronal nuclei clusters human DRG neurons into transcriptomically distinct groups that have been differentially colored. (C) Similarity in expression of differentially expressed genes between human and mouse neuronal types may help functional classification of neuronal types: UMAP representation of human DRG neurons showing relative expression level (blue) of diagnostic markers. For comparison UMAP representation of mouse neurons (Renthal et al., 2020) showing the relative expression patterns of the same markers. In combination, the expression patterns of these and other genes (Figure 1—figure supplements 2 and 4, Figure 1—figure supplement 5) were used to tentatively match several human and mouse transcriptomic classes (Figure 1—figure supplement 4).

One of the best studied groups of somatosensory receptors in mice are nociceptive peptidergic neurons that coexpress a variety of neuropeptides including substance P, calcitonin gene-related peptide (CGRP), and pituitary adenylate-cyclase-activating polypeptide (PACAP). These neurons are typically small soma diameter, nonmyelinated, slow conducting c-fibers, but also include faster conducting lightly myelinated Aδ neurons (Sharma et al., 2020; Nguyen et al., 2019). In the human DRG dataset, TAC1 (substance P), CALCA and CALCB (CGRP), and ADCYAP1 (PACAP), are expressed in several transcriptomic classes (H1, H2, H3, H5, and H6, Figure 1C, Figure 1—figure supplements 4 and 5). For comparison the expression of the same genes in mouse DRG neurons is shown (Figure 1C, Figure 1—figure supplements 4 and 5) using data from single nuclei sequencing (Renthal et al., 2020). Just as in mice, the putative human peptidergic nociceptors express the high affinity nerve growth factor receptor NTRK1, the capsaicin and mustard oil-gated ion channels TRPV1 and TRPA1 but generally only low levels of the stretch-gated ion channel PIEZO2 (Figure 1C, Figure 1—figure supplements 4 and 5).

Although previous localization studies have suggested that in humans the neurofilament protein NEFH is expressed in all sensory neurons (Rostock et al., 2018), this gene showed graded expression in our data (Figure 1C) and marks several classes of cells just as in mice (Figure 1—figure supplement 5). Some of these (including H3 and H6) also express peptidergic markers and the pain-related voltage-gated sodium channel SCN10A (Figure 1—figure supplement 5) and thus have molecular hallmarks of Aδ nociceptors (von Buchholtz et al., 2020). However, the neuronal classes H14 and H15 expressing the highest levels of NEFH are distinct from the peptidergic neurons (Figure 1C, Figure 1—figure supplements 4 and 5), likely representing different types of large diameter, fast conducting myelinated Aβ neurons. These cell types are neurotrophin three receptor NTRK3 positive, some also contain the brain derived neurotrophic factor receptor NTRK2 but exhibit little expression of NTRK1 (Figure 1—figure supplements 4 and 5). In mice, proprioceptors are a subtype of Aβ neurons marked by the calcium binding protein parvalbumin, the transcription factor Etv1 and the voltage-gated sodium channel subunits Scn1a and Scn1b (Sharma et al., 2020; Renthal et al., 2020). In the human data, the small H15 group of NTRK3-positive cells had this expression pattern (Figure 1C, Figure 1—figure supplement 5) implying that proprioceptors have conserved transcriptomic markers in humans and mice. Similarly, small groups of both Aδ-low threshold mechanosensors (H13) and cool responsive neurons (H8) were identified by their characteristic expression profiles of functionally important transcripts (Figure 1—figure supplement 5). Thus, large groups of human and mouse DRG neurons appear to share basic transcriptomic signatures and functional potential, supporting our data as informative about neuronal diversity among human somatosensory neurons.

Despite these broad similarities between the putative peptidergic, proprioceptive, cooling sensitive, Aβ and Aδ classes of DRG neurons in mice and humans there were differences in expression of many genes. These include molecules that modulate cellular responses to internal signals (e.g., growth factor receptors), sensory stimuli and also the mediators they may release. For example, in humans, the H8 putative cool responsive neurons expressing TRPM8 were strongly positive for the BDNF-receptor NTRK2 but hardly expressed the neuropeptide TAC1 whereas in rodents the converse was true (Figure 1—figure supplement 5). Other genes that have been shown to control sensory responses in mice exhibit a different expression pattern in human DRG neurons. For instance, Tmem100 encodes a protein that in mice has been implicated as playing an important role in functional interactions between Trpv1 and Trpa1 and contributing to persistent pain (Weng et al., 2015). By contrast it was almost undetectable in the human sequencing data (Figure 2A). Similarly, we did not detect marked expression of the sphingosine-1-phosphate receptor S1PR3 (Figure 2A) that has been suggested as a target for treating both pain and itch based on mouse work (Hill et al., 2018). More strikingly, a small group of human neurons, H5, expressing TRPA1 were resolved in our clustering (Figure 2B), whereas in mouse nuclear sequencing data no direct counterpart was detected (Figure 1—figure supplement 4). Whole cell-based sequencing (Sharma et al., 2020) of mouse DRG neurons does identify a group of peptidergic nociceptors (called CGRP-gamma) with abundant Trpa1 expression, highlighting the caution needed when interpreting differences across species and the value of alternative approaches. Nonetheless, whereas mouse CGRP-gamma neurons strongly express Calca and Ntrk1, H5 cells are essentially CALCA (CGRP) and NTRK1 negative and instead are strongly NTRK2 positive (Figure 1, Figure 1—figure supplements 4 and 5) suggesting that they may respond differently to external stimuli as well as in their signaling properties and therefore may not be a direct equivalent of mouse CGRP-gamma neurons. Thus, the availability of human transcriptomic data should help focus translational work in model organisms on promising targets with conserved expression patterns in humans.

Figure 2 with 2 supplements see all
Human DRG neurons exhibit specialization that distinguishes them from mouse counterparts.

(A) Universal manifold (UMAP) representation of mouse and human dorsal root ganglia (DRG) neurons showing relative expression level (blue) of two genes that have been linked to pain sensation in mice. Note that both TMEM100 and S1PR3 are more sporadically expressed by the human somatosensory neurons. (B) Classes of DRG neurons that are selectively detected in humans are highlighted together with their expression of key genes. H9 neurons coexpress the cool and mechanosensory ion channels; for comparison cool sensitive neurons (H8) that correspond more closely with their rodent counterparts are also highlighted. (C) Expression profiles of several itch-related genes in the mouse and human DRG transcriptome. (D) Confocal image of a region from a human DRG that was labeled using multiplexed in situ hybridization (ISH) for OSMR, TAC1, and NEFH as indicated in the key. Almost all neurons were OSMR, TAC1, or NEFH positive (Figure 2—figure supplement 1). However, few neurons were strongly positive for more than one of these markers (see Figure 2—figure supplement 2 for individual channels). Note that autofluorescence in all channels from lipofuscin associated with many human neurons should not be confused with real signal (see Figure 2—figure supplement 2 for more detail). Also note that NEFH is typically expressed in larger diameter neurons than the other two markers. Scale bar = 100 µm.

Human DRG neurons without clear transcriptomic equivalents in mice

Analysis of the gene expression patterns of the different classes of human somatosensory neurons revealed several groups for which we could not discern direct counterparts in the mouse. One small but prominent group of human DRG neurons (H9) expresses TRPM8, PIEZO2, SCN10A, and SCN11A (Figure 1—figure supplements 4 and 5, Figure 2B) and clearly segregates from the putative cool sensing cells (H8) that express TRPM8, GPR26, NTM, and FOXP2 but are devoid of both the light touch receptor and the pain-related sodium channels (Figure 2B, Figure 1—figure supplement 5). In mice, Trpm8 expression is simpler with the cool sensing, menthol responsive ion channel just expressed in cells with this latter gene expression pattern (Figure 1—figure supplement 5). Interestingly, single fiber recordings have identified human neurons that respond to both cooling and gentle touch as might be expected for cells expressing both TRPM8 and PIEZO2 (33). Moreover, recent single-cell sequencing of macaque DRG neurons also identified two populations of cells that match H8 and H9 neurons (Kupari et al., 2021). Finally, H9 neurons resemble (but also have differences from) human mechanosensory neurons that were recently engineered by transcriptional programming of stem cells (Nickolls et al., 2020).

A second larger group of human neurons H12 is marked by NTRK3 and the voltage-gated ion channel SCN1A, but is only weakly positive for NEFH, expresses moderate levels of PIEZO2 (Figure 2B, Figure 1—figure supplement 5) and appears distinct from any potential mouse counterpart. The H12 gene expression pattern is most consistent with these cells functioning as a type of mechanosensor that has no direct equivalent in mice. Similarly, we designated H4 as c-nociceptors because of their expression of nociception-related SCN10A and NTRK1 and low level of NEFH (Figure 2B, Figure 1—figure supplement 5). These neurons expressed low levels of neuropeptides, but their overall gene expression patterns did not resemble any mouse counterparts including the nonpeptidergic nociceptors (see below).

The two remaining large groups of neurons in the human dataset H10 and H11 that have no clear mouse counterpart exhibit most similarity with mouse c-type nonpeptidergic neurons (Figure 2—figure supplement 2A). At a functional level both H10 and H11 express receptors that in mice have roles in detecting pruritogens. For example, these clusters were positive for the two subunits (IL31RA and OSMR) of the interleukin 31 receptor and the histamine receptor HRH1 (Figure 2C) that mediate mast cell-related scratching in mice (Solinski et al., 2019). They also express the itch-related neuropeptide NPPB (Figure 2C), nociception-related sodium channels SCN10A and SCN11A as well as TRPV1 (Figure 2—figure supplement 2A) but not appreciable NEFH or TAC1 (Figure 1C). Therefore, it is likely that these are groups of putative unmyelinated, nonpeptidergic nociceptors with potential for triggering human itch responses. H10 and H11 characterization is revisited later in its own section of the results.

The peptidergic nociceptors, myelinated Aβ and Aδ neurons, rarer human-specific cells, and the two nonpeptidergic nociceptor clusters H10 and H11 account for all the neurons in our analysis with H10 and H11 totaling approx. 20% of the neurons. In marked contrast, mouse nonpeptidergic, small diameter neurons are far more numerous than H10 and H11 accounting for 40% of the sensory neurons in mouse DRGs (Renthal et al., 2020) and divide into four highly stereotyped transcriptional groups (Figure 1—figure supplement 4). Two of these classes of mouse neurons (NP2 and NP3) trigger itch (Mishra and Hoon, 2013; Han et al., 2013), one (NP1, expressing Mrgprd) responds to noxious mechanical stimulation (von Buchholtz et al., 2021b). NP1 neurons may have a role in mechanonociception (Gatto et al., 2019) and have recently been associated with suppression of skin inflammation (Zhang et al., 2021), which was hypothesized as relevant for human health. The fourth class corresponds with low threshold mechanosensors (cLTMRs) that are thought to mediate affective touch (Gatto et al., 2019; McGlone et al., 2014). Given this difference between the transcriptomic map of human DRG neurons and their rodent counterparts, we next used independent ISH-based analysis to test basic predictions of the sequencing. If transcriptomic characterization of human DRG neurons is accurate then one clear expectation is that TAC1, NEFH, and OSMR should be expressed by distinct and only partially overlapping populations of human DRG neurons. If it is also comprehensive, that is, not missing equivalent classes to mouse neurons NP1, NP2, and cLTMRs that constitute almost a third of mouse DRG neurons and do not significantly express Nefh, Tac1, or Osmr, then we would anticipate that the same three markers should label the vast majority of neurons. Multigene ISH demonstrates that both these predictions are true for human DRG neurons (Figure 2D, Figure 2—figure supplements 1B and 2B) with essentially every cell labeled by one of these probes but with very few exhibiting strong coexpression. Although NEFH expression could be detected in some of the cells positive for the other markers (Figure 2D, Figure 2—figure supplement 1B), many TAC1- or OSMR-positive small diameter neurons were negative for this neurofilament subunit. Moreover, TAC1 and OSMR labeled almost completely separate sets of cells (Figure 2—figure supplement 1). Notably, in keeping with our assignments based on transcriptomic data, the largest diameter neurons were strongly positive for NEFH whereas TAC1 and OSMR primarily labeled smaller cells (Figure 2D). Finally, these three markers each labeled a large group of neurons (Figure 2—figure supplement 1).

Co-clustering human and mouse DRG-neuron snRNAseq data

As detailed above, the expression of genes that are thought to be important for functional and morphological features of somatosensory neurons reveal similarities between groups of human and mouse neurons. They also expose differences that likely reflect distinct somatosensory adaptations in the two species. We next used co-clustering methods to test whether the wider transcriptome could reveal additional information about the relationships between classes of human and rodent DRG neurons using the same mouse dataset (Renthal et al., 2020) that we analyzed above (Figure 1, Figure 1—figure supplement 4). We used the well-established approach developed by the Satija lab (Stuart et al., 2019) as it has been shown to perform well without forcing false matches. As predicted, several classes of human neurons grouped with corresponding mouse counterparts including H15 – proprioceptors, H14 – Aβ cells, H13 – AδLTMRs, H11 – NP3 (Nppb) neurons, and H3/H6 – Aδ nociceptors (Figure 3A). This analysis suggested that H10 the other cluster that gene expression indicated are also itch related most closely resembled NP1 (Mrgprd) neurons rather than any other human or mouse class of sensory neurons. The H12 cluster, which is human specific, grouped close to larger diameter mouse neurons, whereas other clusters of human cells appeared better aligned with smaller diameter nociceptors. However, all types of peptidergic small diameter nociceptors were less organized in the co-clustering and separated from their potential mouse counterparts despite their qualitatively similar expression of some of the best-known functional markers (Figure 1—figure supplement 4).

Figure 3 with 1 supplement see all
Co-clustering of human and mouse neurons support tentative assignments based on select genes.

(A) Universal manifold (UMAP) representation of the co-clustering of mouse and human neurons. Upper panel shows the mouse neurons colored by their identity when analyzed alone (Figure 1—figure supplement 1); lower panel shows human neurons colored by their identity when analyzed alone (Figure 1). Note that large diameter human neurons match their expected mouse counterparts reasonably well and that the two classes of neurons expressing itch-related transcripts H10 and H11 best match NP1 and NP3 neurons, respectively. (B) Heatmap showing the natural logarithm (see scale bar) of Kullback–Leibler divergences for the various human neuron classes when compared to each class of mouse cells as a reference distribution; potentially human-specific classes based on functional markers are marked by *; see also Figure 3—figure supplement 1.

UMAP plots (Figure 3A) provide a visual representation of similarity between cells with related transcriptomic properties. However, co-clustering methods do not provide a quantitative measure of the similarity between individual cells or clusters of cells. A given cell cluster can be viewed as a multivariate probability distribution in gene expression space. While not commonly employed in gene expression analysis, in probability and information theory, the similarity/dissimilarity between two probability distributions is most commonly measured by calculating their Kullback–Leibler (KL) divergence. Recent advances have allowed KL divergence to be estimated between two samples in continuous multivariate space (Perez-Cruz, 2008) as is observed in dimensionality reduced gene expression data of cell populations. Therefore, we made use of KL divergence estimation to quantitate the similarity between human DRG-neuron clusters and all their potential mouse counterparts (Figure 3B). As expected, clusters that co-segregate in the UMAP analysis showed greatest similarity but additional relationships not apparent from the visual representation of the co-clustering were also seen. For example, the small cluster of human ‘cool’ responsive neurons H8 showed greatest similarity to mouse Trpm8 cells and several groups of human cells (H1, H2, and H5) that gene expression predicted should be c-type peptidergic nociceptors, indeed best matched these cells (Figure 3B). Interestingly, no class of human neurons showed appreciable similarity to mouse cLTMRs. Among the groups of cells that had human-specific gene expression patterns, H5 best matched c-peptidergic nociceptors, H9 (the putative cool and mechanical responsive cells) showed only weak similarity to any mouse neuron class. H12, which we considered likely to be mechanosensors best matched mouse proprioceptors and H4 neurons appeared distantly related to several classes of nociceptor but without a clear match in mice. We also extended this analysis to mouse clustering where similar cell populations had not been combined (Figure 3—figure supplement 1). Although 19 clusters of mouse neurons were now analyzed, neither mouse cLTMRs nor the potentially human-specific classes H4, H9, or H12 better matched a cell types of the other species. One important caveat to this type of analysis remains that any functional conclusions based on shared transcriptomic features still need to be verified experimentally.

Transcriptomically related neurons are spatially grouped in the human dorsal root ganglion

From sequence analysis, we identified a range of potential markers to better explore the diversity of human DRG neurons using ISH. To maximize information, we chose a multiplexed approach that allows localization of up to 12 probes (Figure 4) revealing the different classes of sensory neurons identified in the transcriptomic data. For example, TRPM8 expressing neurons clearly segregate into two distinct types (Figure 4A, Figure 4—figure supplement 1). One set of cells (H8) share other transcriptomic properties with mouse cooling responsive cells: TRPM8, the cool and menthol receptor is not coexpressed with the ion channels SCN10A or PIEZO2 (Figure 4A), but unlike in mice these cells are NTRK2 positive (Figure 4—figure supplement 1). By contrast, other cells (H9) coexpress the pain and light touch related ion channels (SCN10A and PIEZO2) with TRPM8 (Figure 4A, Figure 4—figure supplement 1). Similarly, putative proprioceptive neurons (H15) were distinguished by their expression of NEFH, PIEZO2, and PVALB and lack of NTRK2 (Figure 4B, Figure 4—figure supplement 1). One surprise (Figure 4A, B) was that in small fields of view, several examples of all three of these rare neuron types could be identified in human DRGs. However, much of the rest of the ganglion was devoid of these cell types and instead the neurons there had distinct sets of markers. Therefore, it appears that transcriptomic classes of human DRG sensory neurons may not be stochastically distributed in the ganglion. Indeed, when we examined the distribution of nociceptors and myelinated neurons at lower magnification (using strong selective probes), broad clustering of similar types of neurons was apparent, quantifiable, and statistically significant (Figure 4C, D). We carried out a similar analysis in mouse DRG neurons using Nefh and Scn10a probes and found that there too, cells expressing these markers are not uniformly distributed (Figure 4—figure supplement 2) suggesting spatial clustering of related DRG neurons may be a common feature across species.

Figure 4 with 2 supplements see all
Transcriptomically related classes of human dorsal root ganglia (DRG) neurons are spatially clustered in the ganglion.

Confocal images of sections through a human DRG probed for expression of key markers using multiplexed in situ hybridization (ISH); see Figure 4—figure supplement 1 for the individual panels and additional probes. (A) Left panel shows a group of four H8-neurons (yellow arrows) that express TRPM8 (green) but not PIEZO2 (red) or SCN10A (blue). By contrast, right panel shows a different region of the ganglion where three H9-neurons coexpress these three transcripts (double arrowheads). (B) Other regions of the ganglia were dominated by larger diameter neurons. Putative proprioceptors, highlighted by double arrowheads, expressing PIEZO2 (green) and PVALB (red), but not NTRK2 (blue) were typically highly clustered in the ganglion. (C) Lower magnification images of complete sections stained for NEFH (green) and SCN10A (red) highlight the extensive co-clustering of large and small diameter neurons in different individuals. Scale bars = 100 µm in (A and B); 500 µm in (C). (D) The nearest n neighbors of NEFH and SCN10A single positive neurons in (C) were identified (for n = 1–40): green lines represent proportion (mean, solid line ± standard error of mean [SEM], shaded) of cells surrounding NEFH-positive neurons; red lines, proportion (mean, solid line ± SEM, shaded) of cells surrounding SCN10A-positive neurons. Left panel: proportion of surrounding cells that were NEFH positive. Right panel: proportion of surrounding cells that were SCN10A positive. Dashed black lines are the proportions expected for randomly distributed cells. Insets schematically show a central cell (highlighted by a star) and the surrounding neurons. Neighboring NEFH-positive cells are colored green and SCN10A-positive cells are colored red when these are being scored in the associated graph; gray cells are positive for the other marker. Clustering was statistically significant across the complete range (1–40 neighbors), p ≤ 6.96 × 10−42 (one-tailed Mann–Whitney U-test); n = 803 single positive cells confirming both short and long-range grouping of similar classes of human DRG neurons.

H10 and H11 are distinct but related types of human nonpeptidergic neurons

Perhaps the most intriguing classes of human somatosensory neurons revealed by our transcriptomic approach are the H10 and H11 classes that primarily share features with the mouse nonpeptidergic nociceptors NP1–3 (Figures 2 and 3, Figure 2—figure supplement 2, Figure 5—figure supplement 1). ISH showed that the H10 and H11 classes of neurons, identified by their expression of OSMR were small diameter neurons comparable in size to the TAC1-expressing peptidergic nociceptors (Figure 5A). Two qualitatively different ratios of SCN10A and OSMR were apparent in these cells (Figure 5A) hinting at their distinct identities. Our data (Figure 2, Figure 2—figure supplement 2, Figure 5—figure supplement 1A) show that H10 and H11 neurons express a number of genes that are known markers of mouse NP3 cells and functionally important for triggering pruritic responses (Gatto et al., 2019; Solinski et al., 2019). They are also distinguished from each other by expression of genes that likely play roles in itch and other aspects of somatosensation (Figure 5B, Figure 2—figure supplement 2, Figure 5—figure supplement 2). For example, although not prominently expressed, the human chloroquine responsive receptor MRGPRX1 (27) localized selectively to H10 neurons (Figure 5B) perhaps suggesting a relationship to mouse NP2 cells. By contrast, Janus kinase 1 (JAK1), a mediator of itch through various types of cytokine signaling, including through OSMR (Oetjen et al., 2017), and the neuropeptide SST are particularly strongly expressed in H11 cells (Figure 5B). Both these genes are prominent markers of NP3 pruriceptors in mouse (Figure 5B). However, not all known itch-related transcripts are expressed in H10 and H11 neurons and both classes of cells express genes that better define NP1 neurons in mice as well as other cell types (Figure 5B, Figure 2—figure supplement 2, Figure 5—figure supplement 2).

Figure 5 with 2 supplements see all
Two related classes of human nonpeptidergic small diameter neurons that may mediate itch.

(A) Confocal image of a section through a human dorsal root ganglia (DRG) probed for nociception-related genes using multiplexed in situ hybridization (ISH). In this view, many neurons expressing the itch-related transcript, OSMR (red), are grouped together (arrowheads); cyan arrowheads point to cells that express relatively higher levels of OSMR than SCN10A (green). Peptidergic nociceptors marked by expression of TAC1 (blue) and additional SCN10A-positive cells are also present in this region of the ganglion. (B) Universal manifold (UMAP) representation of mouse and human DRG neurons showing relative expression level (blue) of genes that distinguish H10 and H11 and mark-specific sets of mouse NP1–3 neurons. MRGPRX1 is the human chloroquine receptor and the functional equivalent of Mrgpra3, which in mice marks NP2 cells. Note that coexpression patterns of SST and JAK1 in H11 neurons resembles their expression in mouse NP3 pruriceptors but the ion channel TRPC3 which also marks these cells is primarily expressed in mouse NP1 neurons; see Figure 5—figure supplement 1 for additional breakdown of similarities and differences between nonpeptidergic neurons in mice and humans. (C) Confocal image of the group of candidate pruriceptors shown in (A) probed for expression of genes that distinguish H11 (SST, blue) from H10 cells (PIEZO2, green); note that neurons highlighted with cyan arrowheads have gene expression expected for H11 cells, whereas some H10 cells express lower levels of SST and also exhibit variation in the level of PIEZO2 expression. Scale bars = 100 µm; see Figure 5—figure supplement 1B for individual channels and for expression of additional markers.

H10 cells are also distinguished from H11 and mouse pruriceptors by their prominent expression of the stretch-gated ion channel PIEZO2 (Figure 1C, Figure 2—figure supplement 2). The coexpression of itch-related transcripts and this low threshold mechanosensor hint that H10 neurons may be responsible for the familiar human sensation known as mechanical itch. However, their relationship to NP1-neurons revealed by co-clustering mouse and human data (Figure 3) and their expression of markers for various other cell types (Figure 2—figure supplement 2, Figure 5—figure supplement 2) possibly including nonpeptidergic cLTMRs suggest that their role in somatosensation may not be limited to itch alone.

A problem with single-cell sequencing approaches is the sparse nature of the data making it difficult to disentangle expression level from proportional representation in any cluster. This means that except for the most highly expressed genes, there is inherent ambiguity in interpreting the expression patterns. ISH provides an independent and more analog assessment of expression level that can help resolve this issue. Multiplexed ISH showed that SST divides the OSMR-positive cells into two intermingled types (Figure 5C) in keeping with the sequence data (Figure 5B) and the relative expression patterns of SCN10A and OSMR (Figure 5A). Moreover, the prediction that PIEZO2-OSMR coexpression should mark SST-negative neurons was also largely borne out by ISH (Figure 5C, Figure 5—figure supplement 1). However, ISH also shows that some neurons expressing lower levels of SST are PIEZO2-positive and that some OSMR-positive H10 cells, contain only a very low level of the mechanosensory channel (Figure 5C, Figure 5—figure supplement 1). Therefore, H10 and H11 are by no means homogeneous populations and may not be as functionally distinct as snRNA clustering suggests.

Discussion

Transcriptomic analysis of DRG neurons confirms that mouse and human somatosensory neurons express many of the same genes (Ray et al., 2018). However, although gross similarity in the transcriptomic classification of these cells can be discerned in single-cell analyses (peptidergic versus nonpeptidergic; neurofilament rich, myelinated versus nonmyelinated), the patterns of coordinated gene expression across species are not well conserved and both species exhibit unique specializations. Recently, available transcriptomic data from the macaque further highlights the individuality of somatosensory neurons across species (Kupari et al., 2021). Surprisingly, in that study, despite major differences in gene expression between the monkey and mouse cell types, co-clustering approaches found an apparently close relationship between them including identification of large and distinct groups of NP1, NP2, and NP3 cells and putative cLTMRs (Kupari et al., 2021). However, perhaps because of fragility of large diameter neurons in cell-based droplet sequencing approaches, no candidate Aβ neurons and only very few Aδ cells were recovered (Kupari et al., 2021). By contrast our transcriptomic data and analysis combined with multiplexed ISH provide strong evidence for similarities between large diameter neurons between humans and mice but major differences in small diameter nonpeptidergic neurons as well as the existence of other human-specific cell types. Since in humans all nonpeptidergic neurons express TRPV1 whereas mouse NP1, NP2, and cLTMRs do not prominently express this gene, this fits well with recent data showing much broader expression of TRPV1 in humans than mice (Shiers et al., 2020).

At one level, this type of interspecies variation was unexpected given that there is similarity between the neuronal classes that comprise the mouse lumbar DRGs and trigeminal ganglia despite their very different types of innervation targets (Sharma et al., 2020; Nguyen et al., 2019). However, large changes in the receptive repertoire of other sensory systems have been observed and are thought to play a role in adaptation to specific ecological niches (Yarmolinsky et al., 2009). Thus, the evolution of DRG receptor cell diversity further highlights the importance of appropriate sensory input for fitness and survival of a species. What is unusual relative to other senses is that transcriptomic differences are not limited to just the receptor repertoire for sensing environmental stimuli but instead extend to genes involved in the development and maintenance of defined neuronal subtypes. It is possible that this reflects major differences between mouse and human skin including fur covering. From a translational viewpoint, these differences could explain some of the problems in replicating results from mouse-based therapies (Mogil, 2019; Yezierski and Hansson, 2018) in humans and the availability of the human data may help direct research toward new targets and even suggest precision medicine strategies (e.g., to treat cold pain).

Our analysis identified particularly surprising differences between small diameter nonpeptidergic neurons in mice and humans (H10 and H11 cell types). In mice, one distinctive subset of these cells are the cLTMRs that innervate hairy skin and are thought to be responsible for affective touch (McGlone et al., 2014). At a transcriptomic level, humans do not have a clearly identifiable correlate for these cells (Figure 3B) although careful microneurography has revealed human c-fibers that respond to stroking (Wessberg et al., 2003). We suspect that some of these stroking responsive cells may be the nonpeptidergic H10 neurons that express itch-related genes as well as high levels of PIEZO2. In keeping with this suggestion, a recent preprint of a spatial transcriptomic analysis of human somatosensory gene expression (Tavares-Ferreira et al., 2021) designates cLTMRs as a subset of cells resembling H10 neurons that appear to have lower expression of some pruriceptive markers. However, it is also possible that human cLTMRs are other PIEZO2-expressing neurons that are unique to humans (Figure 2B, Figure 1—figure supplement 3). H10 and H11 neuron classes have more clear similarity to the mouse NP1–3 neurons but also exhibit major differences to all three types of cells. For example, in mice NP1 cells express a large combination of diagnostic markers (Figure 5—figure supplement 2) including Mrgprd that we did not find in our sequencing of human ganglion neurons. Bulk sequencing studies have identified MRGPRD expression in human DRG (Price et al., 2016), but recent ISH localization studies suggest only very low-level expression of this transcript together with MRGPRX1 in somatosensory neurons (Klein et al., 2021). This would fit with our co-clustering that identifies the MRGPRX1-expressing H10 neurons as related to NP1 cells. However, many of the other NP1 markers have potential roles in signal detection and transduction but are not H10 selective (Figure 5—figure supplement 2). Moreover, in mice, Mrgpra3 (the functional equivalent of MRGPRX1) marks the distinct NP2 neurons.

Taken together, our results and analysis suggest that experiments in mice are likely to illustrate general principles that are important for sensory detection and perception in humans but also imply that specific details related both to genes and cell-type responses may differ. In future studies, the central projections and targets of human somatosensory neuron subtypes might provide independent approaches for inferring function. Similarly, using immunohistochemistry (IHC) to understand how these cell classes innervate the skin and other tissues may allow correlation of arborization patterns with microneurography results. Since microneurography can be complemented by microstimulation this could ultimately reveal the role of specific neuronal classes in sensory perception (Ackerley and Watkins, 2018).

Our data provide a searchable database for gene expression in human DRG neurons. However, there are some limitations to the data and interpretation. For example, neither the number of neurons sequenced, nor the depth of sequencing is as comprehensive as for mice (Sharma et al., 2020; Nguyen et al., 2019). This means that rare neuronal subtypes and the expression patterns of even some moderately expressed genes may not be clear therefore it is likely that future studies and larger samples will be needed to refine these issues. Nonetheless, highly multiplexed ISH (Figures 2 and 5) confirm the major findings both about cell types and also gene expression and therefore substantiate the overall value of the data. The nuclear-based sequencing approach used here has advantages in preventing gene expression changes during single-cell isolation and is also likely to be less biased than cell-based approaches in terms of representation of the different cell types (von Buchholtz et al., 2020). However, sn-RNA sequencing provides a somewhat distorted view of cellular gene expression, as has been described for sensory neurons in mice (Nguyen et al., 2019). Therefore, it will be important to confirm expression levels of specific genes using complementary approaches. Finally, any functional roles for neuronal classes identified here have been extrapolated from expression of markers and distant similarity to mouse counterparts. Given the extensive differences that we report, some of these conclusions may need to be revised once cell class can be linked to neuronal function in human subjects and/or using physiological tools in vitro with human tissues.

Materials and methods

Key resources table
Reagent type (species) or resourceDesignationSource or referenceIdentifiersAdditional information
Strain, strain
background
(Mus musculus)
C57BL/6NCrlCharles RiverStrain code: 027Male and female
Antibodyanti-NeuN (rabbit polyclonal)MilliporeCat#ABN78(1:4000 for nuclei isolation)(1:1000 for IHC)
Antibodyanti-Tubb3 (mouse monoclonal)ProteintechCat#66375-1-Ig(1:500)
Antibodyanti-rabbit Cy3-conjugated(donkey polyclonal)JacksonImmuno
Research
Cat#711-166-152(1:1000)
Antibodyanti-mouse FITC-conjugated(donkey polyclonal)JacksonImmuno
Research
cat#715-096-150(1:1000)
Sequence-based reagentChromium Next GEM Single Cell 3′ GEM, Library & Gel Bead Kit v3.1×10 GenomicsCat#1000128
Sequence-based reagentChromium Next GEM Chip G Single Cell Kit×10 GenomicsCat#1000127
Commercial assay or kitRNAscope HiPlex8 Detection reagentsAdvanced Cell
Diagnostics
Cat#324110
Commercial assay or kitRNAscope HiPlex12 Ancillary reagentsAdvanced Cell
Diagnostics
Cat#324120
Commercial assay or kitHuman RNAscope probes (HiPlex 12)Advanced Cell
Diagnostics
NEFH (cat# 448141); TRPM8 (cat# 543121); PIEZO2 (cat# 449951); SCN10A (cat# 406291); NTRK2 (cat# 402621); TAC1 (cat# 310711); OSMR (cat# 537121); SST (cat# 310591); TRPV1 (cat# 415381); PVALB (cat# 422181)
Commercial assay or kitRNAscope Fluorescent Multiplex assayAdvanced Cell
Diagnostics
Cat #320851
Commercial assay or kitMouse RNAscope probes(MultiPlex)Advanced Cell
Diagnostics
Scn10a (cat#426011); Nefh (cat#443671)
Commercial assay or kitHuman RNAscope probes(MultiPlex)Advanced Cell
Diagnostics
OSMR (cat#537121); SST (cat# 310591); HRH1 (cat#416501); NPPB (cat#448511);
Commercial assay or kitNeuroTrace GreenFisher ScientificCat#N21480(1:100)
Software, algorithmCellRanger×10 GenomicsVersion 2.1.1GRCh38.v25.
premRNA
Software, algorithmSeuratSatija labVersions 3-4.04https://satijalab.org/seurat/
Software, algorithmPythonpython.orgVersion 3.7.
Software, algorithmscipy.statsscipy.orgVersion 1.5.2
Software, algorithmscikit-learnscikit-learn.orgVersion 0.23.2
Software, algorithmRStudiohttps://www.rstudio.com/Version 1.4.1106
Software, algorithmRhttps://www.r-project.org/R version 4.1.1
Software, algorithmDoubletFinderhttps://github.com/chris-mcginnis-ucsf/DoubletFinderMcGinnis lab (McGinnis, 2021)
Software, algorithmImageJhttp://imagej.nih.gov/ijImageJ 1.53 c
Software, algorithmAdobe Photoshophttps://www.adobe.com/25.5.1 release
Software, algorithmkl_divergencehttps://gist.github.com/lars-von-buchholtz/636f542ce8d93d5a14ae52a6c538ced5636f542ce8d93d5a14ae52a6c538ced5
OtherRNA-laterThermo FisherCat# AM7021
OtherSpectrum Bessman tissue pulverizerFisher ScientificCat# 08-418-3
OtherDounce homogenizerFisher ScientificCat# 357538
Other40 μm cell strainerThermo FisherCat# 08-771-1
OtherLow bind microfuge tubesSorenson
BioScience
Cat# 11,700
OtherSUPERaseIn RNase inhibitorThermo FisherCat# AM26960.2 U/μl
Otheranti-rabbit IgG microbeadsMiltenyi biotecCat# 130-048-602
OtherLS columnMiltenyi biotecCat# 130-042-401
OtherMouse DRG datasetWoolf labGSE154659Renthal et al., 2020, Neuron

Study design

Request a detailed protocol

Transcriptomic analysis of human DRG neurons was carried out to establish similarities and differences between human somatosensory neurons and their counterparts in model organisms and to provide a resource. We chose a nuclear-based strategy because of its simplicity and quantitative nature relative to isolation of cells (Nguyen et al., 2019). All tissue was obtained from de-identified organ donors and was not preselected or otherwise restricted according to health conditions. We used DRGs from both male and female donors for the sequencing and ISH localization experiments. Randomization and blinding were not used because of the nature of our experiments. Similarly, before starting this study, we had no relevant information for setting sample size for snRNA sequencing from human DRG neurons. Therefore, we stopped data collection when we empirically determined that the cost of adding extra data outweighed the benefit of additional sequencing. In essence, numbers of sn-transcriptomes analyzed were limited by the availability of material and the difficulty of isolating human DRG nuclei with preservation of their transcriptome. We considered that the dataset would serve as a valuable and relatively comprehensive resource once including additional material from an individual preparation made only minor differences to the pattern of clustering we observed. Criteria for data exclusion followed standards in the field (see below) and sample sizes and numbers of replicates are also typical for this type of study and are described in the relevant experimental sections. Apart from the exclusions described for sn experiments, all data obtained were included in our study.

Isolation of human DRG nuclei

Request a detailed protocol

DRG recovery was reviewed by the University of Cincinnati IRB #00003152; Study ID: 2015-5302, title Human dorsal root ganglia and was exempted. Lumbar L4 and L5 DRGs were recovered from donors withing 90 min of cross-clamp (Valtcheva et al., 2016). For ISH and IHC, DRGs were immersion fixed in 4% paraformaldehyde in phosphate-buffered saline (PBS) overnight, cryoprotected in 30% sucrose and were frozen in Optimal Cutting Temperature compound (Tissue Tek). For RNA sequencing, human DRGs immediately were cut into 1–2 mm pieces and stored in RNA-later (Thermo Fisher, Cat# AM7021). Excess RNA-later was removed and the tissue was frozen on dry ice and stored at −80°C. Nuclei were isolated from each donor separately as described previously (Nguyen et al., 2019) with minor modification. Briefly, the tissues were homogenized with a Spectrum Bessman tissue pulverizer (Fisher Scientific, CAT# 08-418-3) in liquid nitrogen. The sample was then transferred to a Dounce homogenizer (Fisher Scientific, Cat# 357538) in 1 ml of freshly prepared ice-cold homogenization buffer (250 mM sucrose, 25 mM KCl, 5 mM MgCl2, 10 mM Tris, pH 8.0, 1 µM DTT, 0.1% Triton X-100 [vol/vol]). To lyse cells and preserve nuclei homogenization used 5 strokes with the ‘loose’ pestle (A) and 15 strokes with the ‘tight’ pestle (B). The homogenate was filtered through a 40 µm cell strainer (Thermo Fisher, cat# 08-771-1), was transferred to low bind microfuge tubes (Sorenson BioScience, cat# 11700) and centrifuged at 800 g for 8 min at 4°C. The supernatant was removed, the pellet gently resuspended in 1 ml of PBS with 1% bovine serum Albumin (BSA) and SUPERaseIn RNase Inhibitor (0.2 U/µl; Thermo Fisher, Cat#AM2696) and incubated on ice for 10 min.

Neuronal nuclei selection was performed by incubating the sample with a rabbit polyclonal anti-NeuN antibody (Millipore, cat#ABN78) at 1:4000 dilution with rotation at 4°C for 30 min. The sample was then washed with 1 ml of PBS with 1% BSA and SUPERaseIn RNase Inhibitor and centrifuged at 800 × g for 8 min at 4°C. The resulting pellet was resuspended in 80 µl of PBS, 0.5% BSA, and 2 mM EDTA. 20 µl of anti-rabbit IgG microbeads (Miltenyi biotec, cat# 130-048-602) were added to the sample followed by a 20-min incubation at 4°C. Nuclei with attached microbeads were isolated using an LS column (Miltenyi Biotec, cat# 130-042-401) according to the manufacturer’s instruction. The neuronal nuclei enriched eluate was centrifuged at 500 × g for 10 min, 4°C. The supernatant was discarded, and the pellet was resuspended in 1.5 ml of PBS with 1% BSA. To disrupt any clumped nuclei, the sample was homogenized on ice with an Ultra-Turrax homogenizer (setting 1) for 30 s. An aliquot was then stained with trypan blue and the nuclei were counted using a hemocytometer. The nuclei were pelleted at 800 × g, 8 min at 4°C and resuspended in an appropriate volume for ×10 chromium capture. A second count was performed to confirm nuclei concentration and for visual inspection of nuclei quality. Note that prior to sequencing no check was made that this selection procedure was unbiased but we have previously used the approach to purify mouse trigeminal ganglion neuron nuclei (Nguyen et al., 2019) and have demonstrated highly quantitative recovery of all neuronal types (von Buchholtz et al., 2020). IHC analysis (Figure 1—figure supplement 1) supports the use of NeuN enrichment as a means to purify human DRG neurons and together with ISH quantitation (Figure 2—figure supplement 1) substantiates the relatively unbiased purification of various classes of neuronal nuclei using this method.

Single nuclear capture, sequencing, and data analysis

Request a detailed protocol

10x-chromium capture and library generation were performed according to the manufacturer’s instructions using v3 chemistry kits. Next generation sequencing was performed using Illumina sequencers. 10x chromium data were mapped using CellRanger to a pre-mRNA modified human genome (GRCh38.v25.premRNA). Data analysis used the Seurat V3 packages developed by the Satija lab and followed standard procedures for co-clustering (Stuart et al., 2019). For sn-RNA sequencing experiments cell filtering was performed as follows: outliers were identified and removed based on the number of expressed genes (500–10,000 retained) and mitochondrial proportion (<10% retained). Normalization and variance stabilization used regularized negative binomial regression (sctransform). After initial co-clustering of data from the different preparations, non-neuronal cell clusters were identified by their gene expression profiles. Clusters not expressing high levels of neuronal or somatosensory genes like SNAP25, SCN9A, SCN10A, PIEZO2, NEFH, etc. but instead expressing elevated levels of markers of non-neuronal cells including PRP1, MBP, QKI, LPAR1, and APOE were tagged as non-neuronal and were removed to allow reclustering of ‘purified’ human DRG neurons. A total of 1837 human DRG neuronal nuclei were included in the analysis (Table 1). The mean number of genes detected per nucleus was 2839 (range 501–9652), with a standard deviation of 1917. Doublet detection was performed on the individual datasets using DoubletFinder (McGinnis et al., 2019). For the clustering shown in the main figures the small number of potential doublets (Figure 1—figure supplement 2E) were not removed; principal components (PCs) were determined from integrated assay data and PCs 1–16 were used both for UMAP display of the data and for determining clusters. The resolution for clustering used relatively low stringency (2.0) and closely related clusters without distinguishing markers were merged. Changes in display clustering parameters and in the cutoffs for data inclusion/exclusion as well as leaving out nuclei from any single preparation made differences in how the data were represented graphically and the number of clusters identified but not to the main conclusions (Figure 1—figure supplement 3). All the different transcriptomically related neuron types described here could still be readily discerned in UMAP analysis of expression data. Identification of markers used the FindMarkers and FindAllMarkers functions of Seurat using default settings and with markers limited by a minimal level of expression in the positive cluster set using (min.pct = 0–0.3) and the difference in expression between positive and negative clusters set using (min.dif.pct = 0–0.3).

Table 1
Sequencing data for the five individual DRG nuclear preparations.
Preparation#DRG1-F36DRG2-M36DRG3a-F34DRG3b-F34DRG4-F35DRG5-F55
Number of reads273,123,313121,333,026105,113,253101,427,23490,306,79394,069,716
Valid barcodes97.50%97.60%98.20%96.10%93.60%97.90%
Sequencing saturation93.20%93.10%39.70%82.10%38.70%55.90%
Q30 bases in barcode96.30%97.40%97.80%97.50%97.70%96.80%
Q30 bases in RNA read88.30%85.90%87.40%91.50%93.40%92.20%
Q30 bases in sample index96.00%96.80%97.40%96.60%96.60%94.80%
Q30 bases in UMI95.50%97.00%97.60%97.00%97.20%96.10%
Reads mapped confidently to genome86.50%81.20%71.20%16.90%10.80%18.10%
Reads mapped confidently to intergenic regions4.10%4.20%4.00%1.30%0.90%1.20%
Estimated number of cells5842736,180223872999
Mean reads per cell467,676444,44317,008454,830135,18994,163
Median genes per cell1917709786481872812
Total genes detected24,64617,55928,75917,33824,79825,027
Median UMI counts per cell29298629656631,293988
DRG cells in final object21215277080281342

For analysis of the mouse, a random subset of data from sn-RNA sequencing of DRGs from wild type mice were extracted from data deposited by the Woolf lab (Renthal et al., 2020) using the R-function: sample. The data were filtered according to gene count (400–12,000 retained) and mitochondrial DNA (<1% retained) leaving 6895 DRG nuclei that were clustered using standard methods (Stuart et al., 2019). For the data that are displayed in most figures, PCs 1–20 were used for UMAP display with resolution for clustering set at 3.5; closely related clusters were merged. For Figure 3—figure supplement 1, PCs 1–30 were used for UMAP display with resolution for clustering set at 1.6; related clusters were not merged. The expression patterns that are described for genes in mice can also be checked in the outstanding and easy to search single-cell analysis provided by Sharma et al., 2020.

Co-clustering of mouse and human data used methods described by the Satija lab (Stuart et al., 2019). Briefly, we capitalized gene names in both mouse and human DRG-neuron data and then limited the datasets to genes that were shared between them. Mouse and human data were individually normalized and subjected to variance stabilization (sctransform) and were then combined using default parameters for integration (dims = 15). After integration, data were scaled and 30 PCs were calculated, used for UMAP display and clustering. Clustering resolution was set at 0.5 and cells were color coded by transferring their identity in the original clustering of human or mouse data. We experimented using different numbers of mouse neurons and found broadly similar results when approx. 1800, 3500, or 7000 mouse nuclei were used. However, although the relationships shown in Figure 3 could be discerned using the lower numbers of mouse nuclei the mouse neurons were less organized (with some of the mouse clusters splitting). When we used substantially greater numbers of mouse neurons from the full Renthal dataset (Renthal et al., 2020), the mouse neurons dominated the clustering and only co-clustering of large diameter neurons across species was observed with the other human and mouse nuclei clustering separately from each other. In order to quantify the similarity/dissimilarity between a given mouse cluster and any human cluster, the KL divergence between their distributions in 30-dimensional continuous space (PCs 1–30 from the co-clustering) was estimated as described previously (Perez-Cruz, 2008) using R (https://gist.github.com/lars-von-buchholtz/636f542ce8d93d5a14ae52a6c538ced5; von Buchholtz, 2021a). The natural logarithm of the KL divergences for each mouse/human pairing was plotted as a heatmap in R.

ISH and IHC

Request a detailed protocol

Cryosections from human DRGs (from different donors than those used for sequencing) were cut at 20 µm and used for ISH with the RNAscope HiPlex Assay (Advanced Cell Diagnostics) following the manufacturer’s instructions. The following probes were used: NEFH (cat# 448141); TRPM8 (cat# 543121); PIEZO2 (cat# 449951); SCN10A (cat# 406291); NTRK2 (cat# 402621); TAC1 (cat# 310711); OSMR (cat# 537121); SST (cat# 310591); TRPV1 (cat# 415381); and PVALB (cat# 422181). We also used RNAscope Multiplex Fluorescent Assay (Advanced Cell Diagnostics) for localization of OSMR (cat# 537121); NPPB (cat# 448511); HRH1 (cat# 416501); and SST (cat# 310591).

Confocal microscopy (5 µm spaced optical sections for HiPlex and 1 µm spaced optical sections for Multiplex Assays) was performed with a Nikon C2 Eclipse Ti (Nikon) using a ×40 objective. All confocal images shown are collapsed (maximum projection) stacks. HiPlex images were aligned and adjusted for brightness and contrast in ImageJ as previously described (von Buchholtz et al., 2020). Diagnostic probe combinations were used on at least three sections from at least two different individuals with qualitatively similar results. Overall, we used sections from ganglia from five different donors, however, the signal intensity of all probes varied between the individual ganglia making identification of positive signal risky in some cases. Strongest signals were observed for sections from an 18-year-old male donor and 23- and 35-year-old female donors. All images displayed here, and our analysis including cell counts were only from sections of ganglia isolated from these three individuals. When individual channels are displayed in the supplements, the strongest autofluorescence signals have been selected and superimposed using photoshop to help focus attention on ISH signal.

IHC was carried out using 20 µm sections of human DRGs using standard procedures. Primary antibodies to NeuN (Millipore) and β3 tubulin (Proteintech) were used at 1:1000 and 1:500 dilution, respectively, and were visualized using appropriate secondary antibodies conjugated to fluorescent reporters (Jackson Immuno Research, 1:1000). Sections were counterstained with 4′,6-diamidino-2-phenylindole (DAPI) and/or green NeuroTrace (Fisher Scientific). Confocal images were acquired with a Nikon C2 Eclipse Ti (Nikon) using a ×40 objective with optical sections 5 µm apart. Images are maximum projection (collapsed) stacks of individual optical sections; consistency of NeuN staining was assessed using sections from two individuals; images were processed using Adobe Photoshop CC to adjust brightness, contrast, and set channel color for display.

Animal experiments were carried out in strict accordance with the US National Institutes of Health (NIH) guidelines for the care and use of laboratory animals and were approved by the National Institute of Dental and Craniofacial Research animal care and use committee (protocol #20-1041). Male and female mice were used for all experiments but were not analyzed separately. Mice were C57BL/6NCrl and were 6 weeks or older. Cryosections from fresh frozen mouse lumbar DRGs in OCT (Tissue-Tek) were cut at 10 µm and used for ISH with RNAscope Multiplex Fluorescent Assay (Advanced Cell Diagnostics) according to the manufacturer’s instructions. Confocal images were acquired with a Nikon C2 Eclipse Ti (Nikon) using a ×10 objective at 1 µm optical section. Images are maximum projection (collapsed) stacks of 10 individual optical sections; consistency of staining was assessed using multiple sections from at least three mice as is considered standard; images were processed using Adobe Photoshop CC to adjust brightness, contrast, and set channel color for display.

Spatial analysis of cell clusters

Request a detailed protocol

In order to quantify spatial clustering of cell types, neurons in ISH images (for humans, one male, one female and for mice several sections from different DRGs and animals) were manually outlined and annotated as NEFH-only or SCN10A-only; cells expressing both genes were not counted. Centroid coordinates of these cells and their distances were analyzed in Python 3.7. The nearest neighbors were identified based on Euclidean distance (Scikit-Learn package) and the percentage of NEFH and SCN10A cells in each neighborhood of size 1–40 cells was calculated. Statistical significance between NEFH- and SCN10A-surrounding neighborhoods was determined using a one-tailed Mann–Whitney U-test (Scipy Stats package).

Data availability

Sequence data is available in GEO, accession number GSE168243; a searchable version of the data will also be made available at https://seqseek.ninds.nih.gov/home.

The following data sets were generated
    1. Nguyen M
    2. Ryba N
    3. Davidson S
    (2021) NCBI Gene Expression Omnibus
    ID GSE168243. Single nucleus transcriptomic analysis of human dorsal root ganglion neurons.
The following previously published data sets were used
    1. Renthal W
    2. Yang L
    3. Tochitsky I
    4. Woolf C
    (2020) NCBI Gene Expression Omnibus
    ID GSE154659. Transcriptional reprogramming of distinct peripheral sensory neuron subtypes after axonal injury.

References

    1. Mogil JS
    (2019) The translatability of pain across species
    Philosophical Transactions of the Royal Society B 374:20190286.
    https://doi.org/10.1098/rstb.2019.0286
  1. Conference
    1. Perez-Cruz F
    (2008) 2008 IEEE International Symposium on Information Theory - ISIT
    Kullback-Leibler Divergence Estimation of Continuous Distributions. pp. 1666–1670.
    https://doi.org/10.1109/ISIT.2008.4595271

Decision letter

  1. Rebecca Seal
    Reviewing Editor; University of Pittsburgh School of Medicine, United States
  2. Catherine Dulac
    Senior Editor; Harvard University, United States

Our editorial process produces two outputs: i) public reviews designed to be posted alongside the preprint for the benefit of readers; ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.

Decision letter after peer review:

Thank you for submitting your article "Single nucleus transcriptomic analysis of human dorsal root ganglion neurons" for consideration by eLife. Your article has been reviewed by 4 peer reviewers, one of whom is a member of our Board of Reviewing Editors, and the evaluation has been overseen by Catherine Dulac as the Senior Editor. The reviewers have opted to remain anonymous.

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

This is clearly an important and timely study of the gene expression profiles of individual human primary sensory neurons. The authors have acute access to quickly harvest tissue from donors which is extremely beneficial. While the number of neurons and depth of sequencing is on the low side, the authors provide new insights into the molecular organization as well as a comparison to mouse which will be valuable for the field. To be suitable for publication in eLife, the reviewers request that the manuscript is first strengthened by providing additional data, analysis and clarifications as outlined below.

Essential revisions:

1. There is concern that additional evidence is needed to support two major claims of the paper: that there are mouse-specific and human-specific neuron subtypes. These are issues stemming from the bioinformatics and from comparing datasets that were not generated in parallel. The in situ hybridization experiments provide some support for species specificity, but expression of individual marker genes may evolve while the overall cell population remains consistent. To address these issues, the authors should:

a) Add doublet detection;

b) Attempt different levels of clustering (number of clusters) and show the results are robust;

c) Try different methods and parameters of co-clustering across species;

d) Try a method of comparing across species that doesn't depend on the exact clusters same clusters found in the analysis. For example, this might involve training a classifier in mouse to predict cell type from gene expression. Then, applying it to human and showing that clearly conserved cell types show up, but that none of the species-specific populations are correctly labeled.

2. The reviewers had concern about the NeuN method of isolation. Was the technique in any way validated given the large number of non-neuronal cells and could there also have been bias in neurons that were selected by using this method, as it was pointed out that the levels of NeuN can vary across neuronal subtypes? Evidence validating the sorting and evidence against bias should be included, including plotting the distribution of UMIs and genes as well as using more than one glial marker.

3. The number of nuclei that were sequenced and the depth of sequencing are low. Can the authors make a short statement indicating that although this work enables a better understanding of the subsets of human sensory neurons, the readers should keep in mind that this classification will be further refined as more studies or larger samples are added

4. The authors should perform in situ for NPPB and also quantify all of the in situ data.

5. The execution of the experiment looking at organization of NEHF and SCN10A positive cells in the human DRG is not entirely compelling based on the few images shown but the quantification is potentially better. Depending on how strongly the authors want to claim the "spatial clustering" model, the authors should perform additional quantification/modeling and also similar analysis with mouse DRG neurons.

Reviewer #3 (Recommendations for the authors):

The most obvious factor to strengthen this paper would be to increase the number of nuclei used in the clustering analysis. This may be complicated by experimental issues, such as the need to run samples in parallel. They can perhaps make a short statement indicating that although this work enables a better understanding of the subsets of human sensory neurons, that the readers should keep in mind that this classification will be further refined as more studies or larger samples are added.

One other note is that there is less interest in the multiple comparisons to the mouse classification, but more on defining the sensory profiles in human. Just describing how one class of sensory neurons has the possibility to be polymodal based on its expression of different receptors, and build on that in the Results section. The part about adaptation to ecological niches in the discussion is exactly what this paper can bring that is unique and exciting.

For instance, I find this very interesting that sensory neurons innervating hairy vs. non-hairy skin neurons get prioritized for human somatosensation compared to the mouse.

On another note: do the authors know whether they aren't biasing their analysis for neurons of large or small diameter, they shouldn't be (since there isn't really a size filtration step) but maybe it would be worthwhile for them to do a sort of population study where they look at the overall diameters of neurons and generate ratios of large to medium to small diameter neurons in the human DRG and compare it to mouse. Then maybe comparing the markers used for a FISH to generate ratios of how much those cells occur in their snRNASeq data to confirm that they aren't accidentally biasing their analysis.

Figure 1: they need to state the number of nuclei analyzed in the text, not just the figure legend.

Figure 1: they should show a table with the number of nuclei from each donor.

Lines 379-381: they acknowledge that they have higher neuron recovery, but the tradeoff is lower read depth – wouldn't read depth be of the utmost priority when making claims about new classes of sensory neurons in humans?

Figure 2D: The authors claim that most neurons are labelled, yet from the image it seems like there are major gaps (that look like places where we'd find large diameter neurons based on the size). Did they use DAPI or any other counterstain at all to illustrate what the number of neurons actually was?

The authors claim they do not detect Mrgprd. Can they comment on the accuracy of snRNAseq at detecting low level transcripts? Is it that other low-level transcripts are undetectable? Can they comment on whether this could be a limitation of the technique?

Can the authors discuss the Shiers study published in Bioarchives from the lab of Ted Price, where they found that TRPV1 was broadly expressed?

Overall, I find this manuscript to be very important for the field of pain research and will be accessed frequently by pain scientists.

Reviewer #4 (Recommendations for the authors):

The number of cells and the depth of the sequencing are on the low end. One potential consequence is that populations such as c-LTRMs, which were reported by mouse studies and by the macaque and other human study are missing. Thus, it likely that there are not enough data to generate a comprehensive representation of all of the cell types. A bioinformatic comparison to the macaque would help to clarify some issues with what is different between these two important studies. In any case, a table that compares the cell types identified in the macaque and the other human study with this one would be very helpful. Information about the number of nuclei and "depth of sequencing" comparison across studies within this table or figure would also be helpful. The number of neuronal nuclei and genes per nuclei need to be stated upfront in the results.

The absence of C-LTMRs in the data is surprising. In mouse, C-LTMRs share features and genes with non-peptidergic neurons, are cooling sensitive and mechanically sensitive though the cooling sensitivity in mice is not likely mediated by Trpm8. How does the macaque c-LTMR cluster compare to the human data? How do the C-LTMRs assigned by the Price lab study compare?

For the comparison to mouse, the authors state that a random subset of Renthal data was used. Perhaps a better description of this would be helpful particularly what was in the subset? It is unclear why the full set would enrich for larger neurons.

Authors list genes that are not found in human data but are in mouse as well as a population that is in human that may not be in mouse. Page 8. The authors should use RT-PCR to see if the gene isn't expressed or if there is efficient nuclear export issue or annotation issue.

The authors should quantify the in situ data.

Confusing about how H5 has a mouse counterpart with a few different genes (single cell comparison) but then H9 doesn't have a mouse counterpart. Were computational methods used to compare to the single cell mouse data?

The logic of sentence on lines 217-220 escapes me. If the snRNA.Seq representation was not comprehensive I assume that you could get the same result – all cells have one of these genes and they are generally non-overlapping.

Organization to the DRG would make sense as it is a common theme for sensory systems and for the somatosensory system. The data are in the supplementary where the images are not entirely convincing (3-D reconstruction of 3 or 4 DRG would be more convincing) but the analysis seems to suggest this could be true. The authors suggest this organization is present in human but not mouse. The distribution of myelinated and unmyelinated should be shown for mouse using the same approach.

Reviewer #5 (Recommendations for the authors):

– It's a bit odd that quality of the NeuN sorting wasn't experimentally validated, especially given the large numbers of non-neuronal nuclei for a sorted experiment. The lack of structure in the non-neuronal cells could indicate that these droplets were empty or contained highly damaged cells. Plotting the distribution of UMIs and genes, as well as more than one glial marker could help resolve this.

– Levels of NeuN itself can vary across neurons subtypes (https://celltypes.brain-map.org/rnaseq/human_m1_10x). Could these be responsible for the differences in cell type across species?

– Seurat has different options for clustering included, each of them with a number of parameter options. There is not enough detail in the methods to identify if it was appropriate. In particular, how were the number of clusters identified? Do they really reflect distinct cell types or are there continuous signals? Similarly, the parameters outlier removal aren't described. These are especially critical for post-mortem experiments that can contain damaged cells.

– The length of time between subject death and sample collection should be included.

– There should be computation procedure for identifying/removing doublets.

– The number of and distribution of UMIs per nucleus in the human dataset should be reported.

– What parameters and function were used for the co-clustering between human and mouse?

– The KL divergence calculated identifies similarities for cell classes and seems to perform reasonably. However, I am not aware of that test being commonly used to link cell populations across species. What is the justification for using it.

– The procedure for identifying markers within and across species is not well-described. In particular, differences in the representation of different cell types across species could lead to artifacts in which markers are identified. Furthermore, since much of manuscript is devoted to differences between mouse and human expression, a statistical procedure should be used to verify the ability of particular genes to label a population in one species versus another.

[Editors' note: further revisions were suggested prior to acceptance, as described below.]

Thank you for resubmitting your work entitled "Single nucleus transcriptomic analysis of human dorsal root ganglion neurons" for further consideration by eLife. Your revised article has been reviewed by 3 peer reviewers and the evaluation has been overseen by Catherine Dulac as the Senior Editor, and a Reviewing Editor.

The manuscript has been significantly improved. There is one remaining issue that still needs to be further clarified or the claim tempered; namely there is concern that the authors have not convincingly proven that species specific cell types exist. Although the authors have argued the limitations of the bioinformatic approaches, it is suggested that they try generating substantially smaller clusters of cells in mouse (~twice the number of clusters). In this case, if none of these mouse sub-populations have a match to the human clusters based on KL divergence, that would substantially strengthen the argument that the human population is not conserved in mouse.

https://doi.org/10.7554/eLife.71752.sa1

Author response

Essential revisions:

1. There is concern that additional evidence is needed to support two major claims of the paper: that there are mouse-specific and human-specific neuron subtypes. These are issues stemming from the bioinformatics and from comparing datasets that were not generated in parallel. The in situ hybridization experiments provide some support for species specificity, but expression of individual marker genes may evolve while the overall cell population remains consistent. To address these issues, the authors should:

a) Add doublet detection;

We added doublet detection Figure 1—figure supplement 2E.

b) Attempt different levels of clustering (number of clusters) and show the results are robust;

We added different clustering parameters to show robustness (Figure 1—figure supplement 3).

c) Try different methods and parameters of co-clustering across species;

Co-clustering of DRG data across species is a significant challenge because of the dramatic differences in gene expression even amongst the most conserved classes. As stated in the paper the Satija lab methodology was developed to find similarities across datasets without forcing matches. We also used CONOS: a program that provides significant flexibility for setting parameters. Using PCAs or PCCAs CONOS did not align mouse and human DRG data. With CCAs and parameters as used by Kupari et al., human and mouse DRG data appear to be closely related. However, these parameters inappropriately match human cortical excitatory and inhibitory neurons, DRG neurons with cortical excitatory neurons and even peripheral blood mononuclear cells with brain inhibitory neuron nuclei. Therefore, we have not included the analysis in our revised manuscript.

d) Try a method of comparing across species that doesn't depend on the exact clusters same clusters found in the analysis. For example, this might involve training a classifier in mouse to predict cell type from gene expression. Then, applying it to human and showing that clearly conserved cell types show up, but that none of the species-specific populations are correctly labeled.

It is not entirely clear what the reviewers are asking for; Kupari et al. trained a classifier with mouse data and then used it on macaque data (which seems to be what the reviewers want us to do), but this method forces a cell to “match” one of the mouse neuron types and has no control to ascertain validity. A serious concern is that all statistical evaluation of such an analysis assumes the basic tenet that testing data come from the same distribution as training data. This is clearly violated in across species comparisons. Moreover, the analysis does not evaluate closeness of relationships: the probability of a match is not a measure of similarity but of ambiguity between the limited prediction choices. For example, human cells in a cluster with gene expression features resembling three classes of neurons in mouse would distribute amongst the three clusters because of random differences at the single cell level. By contrast, cells in a cluster with very little similarity to any type of mouse neurons may be assigned to a single class because each cell needs to be matched to one of the mouse classes. We have not come up with an alternative that overcomes these problems. This is why we used Kullback-Leibler divergence to examine similarities between the robust clusters identified in mouse and human data. We now explain the appropriateness and rationale for using this method in the Results section (lines 234-241; line numbers refer to the word document with track changes all markup selected).

2. The reviewers had concern about the NeuN method of isolation. Was the technique in any way validated given the large number of non-neuronal cells and could there also have been bias in neurons that were selected by using this method, as it was pointed out that the levels of NeuN can vary across neuronal subtypes? Evidence validating the sorting and evidence against bias should be included, including plotting the distribution of UMIs and genes as well as using more than one glial marker.

We have used this approach in mice but had not tested if the NeuN method was biased in humans before use. However, we needed to deplete the roughly 100-fold excess of non-neuronal nuclei in human DRGs. We have now carried out immunohistochemistry using the NeuN antibody and demonstrate that the vast majority of neurons (but not other cells) are robustly stained in sections (new Figure 1—figure supplement 1). We believe this will be valuable information for readers and should assuage reviewers’ concerns.

3. The number of nuclei that were sequenced and the depth of sequencing are low. Can the authors make a short statement indicating that although this work enables a better understanding of the subsets of human sensory neurons, the readers should keep in mind that this classification will be further refined as more studies or larger samples are added

Statement added to discussion (lines 427-428)

4. The authors should perform in situ for NPPB and also quantify all of the in situ data.

NPPB expression is not a focus of this paper and is only mentioned once in passing in the results. Nonetheless, we obtained relevant probes and carried out this request. As expected from the sequence data and in agreement with previous ISH data published by other groups, our data demonstrate low-level expression of NPPB in H10 and H11 neurons (new Figure 5—figure supplement 1A).

At the reviewers’ request, we also include extensive quantitation of ISH data (new Figure 2—figure supplement 1). Despite recent trends to “quantitate” all ISH data, we are not keen on this level of analysis (unless it can be done by automated signal detection) because many factors can influence the scoring of positives. Extensive autofluorescence makes automated signal detection impractical for human DRGs, thus we concentrated on diagnostic genes.

5. The execution of the experiment looking at organization of NEHF and SCN10A positive cells in the human DRG is not entirely compelling based on the few images shown but the quantification is potentially better. Depending on how strongly the authors want to claim the "spatial clustering" model, the authors should perform additional quantification/modeling and also similar analysis with mouse DRG neurons.

This is not a major point of the paper but helps explain the anatomically clustered nature of rare but related cell-types including proprioceptors and “cool sensing” cells. Nearest neighbor analysis provides a strong argument for organization at the level of the large diameter non-nociceptive neurons (NEFH-only) and the smaller nociceptors (SCN10A-only): for up to 40 nearest neighbors of 803 single positive cells in 2 sections, the maximum p-value was < 7 x 10-42. We have moved this analysis to the main figure to highlight the quantitative measure. We carried out the same analysis in mice Figure 4—figure supplement 2 as requested by the referees. Our data demonstrate that a similar organization is also present there strengthening the conclusions by showing conservation across species.

Reviewer #3 (Recommendations for the authors):

[…] Figure 1: they need to state the number of nuclei analyzed in the text, not just the figure legend.

Figure 1: they should show a table with the number of nuclei from each donor.

We moved the number to the text and list the nuclei from each donor in the table.

Figure 2D: The authors claim that most neurons are labelled, yet from the image it seems like there are major gaps (that look like places where we'd find large diameter neurons based on the size). Did they use DAPI or any other counterstain at all to illustrate what the number of neurons actually was?

Neurons are well spaced in human DRGs and surrounded by nerve fibers and other structures; we provide an image in the supplement where the brightness/contrast has been adjusted to highlight this point. See also new Figure 1—figure supplement 1 for DAPI and other counterstaining.

The authors claim they do not detect Mrgprd. Can they comment on the accuracy of snRNAseq at detecting low level transcripts? Is it that other low-level transcripts are undetectable? Can they comment on whether this could be a limitation of the technique?

We already addressed the limitations of snRNAseq in the paper and explained that the low expression of MRGPRD in human DRG neurons was consistent with recently published ISH images for this transcript showing expression close to the detection limit of the technique. By contrast, in mouse Mrgprd-expression is very robust.

Can the authors discuss the Shiers study published in Bioarchives from the lab of Ted Price, where they found that TRPV1 was broadly expressed?

Our data and theirs are consistent; we have added a sentence to the discussion and reference to the now published Shiers paper (ref-43, lines 346-349).

Reviewer #4 (Recommendations for the authors):

The number of cells and the depth of the sequencing are on the low end. One potential consequence is that populations such as c-LTRMs, which were reported by mouse studies and by the macaque and other human study are missing. Thus, it likely that there are not enough data to generate a comprehensive representation of all of the cell types. A bioinformatic comparison to the macaque would help to clarify some issues with what is different between these two important studies. In any case, a table that compares the cell types identified in the macaque and the other human study with this one would be very helpful. Information about the number of nuclei and "depth of sequencing" comparison across studies within this table or figure would also be helpful. The number of neuronal nuclei and genes per nuclei need to be stated upfront in the results.

The absence of C-LTMRs in the data is surprising. In mouse, C-LTMRs share features and genes with non-peptidergic neurons, are cooling sensitive and mechanically sensitive though the cooling sensitivity in mice is not likely mediated by Trpm8. How does the macaque c-LTMR cluster compare to the human data? How do the C-LTMRs assigned by the Price lab study compare?

Note we do not think that humans have no cLTMRs merely that they don’t have an identifiable group of cells that correspond with the mouse neurons that are thought to play this role (see lines 391-400). We have added some discussion about the macaque study (lines 160-161 and 335-402) but the differences in approaches (e.g., cells vs nuclei) make extensive comparisons difficult. The Price lab study has not yet been published and the original bioRxiv version did not assign any neurons as cLTMRs. This was changed at a later time making it dangerous to extensively evaluate similarities and differences between our study and theirs since their interpretation may yet change. Therefore, we have removed most of the discussion we had about the spatial transcriptomic analysis but now point out (lines 397-400) that our speculation that H10 neurons may function as cLTMRs is consistent with their assignment of transcriptomically similar cells as cLTMRs.

For the comparison to mouse, the authors state that a random subset of Renthal data was used. Perhaps a better description of this would be helpful particularly what was in the subset? It is unclear why the full set would enrich for larger neurons.

We added more detail to the methods. The full Renthal dataset does not enrich for larger neurons but clusters in the same way as smaller subsets; we used the smaller group to prevent the mouse data overwhelming the human data in co-clustering as is now explained more carefully in the methods. The smaller subsets of mouse data (including just 1800 nuclei like our human dataset) are also (1) just as informative as the full dataset in terms of clustering and (2) more efficient to use computationally.

The logic of sentence on lines 217-220 escapes me. If the snRNA.Seq representation was not comprehensive I assume that you could get the same result – all cells have one of these genes and they are generally non-overlapping.

Thanks, we have rewritten the sentence to better explain our interest in these 3 genes.

Reviewer #5 (Recommendations for the authors):

[…] – The number of and distribution of UMIs per nucleus in the human dataset should be reported.

Added to Figure 1—figure supplement 2D.

– What parameters and function were used for the co-clustering between human and mouse?

– The procedure for identifying markers within and across species is not well-described. In particular, differences in the representation of different cell types across species could lead to artifacts in which markers are identified.

We have added more detail to the methods.

[Editors' note: further revisions were suggested prior to acceptance, as described below.]

The manuscript has been significantly improved. There is one remaining issue that still needs to be further clarified or the claim tempered; namely there is concern that the authors have not convincingly proven that species specific cell types exist. Although the authors have argued the limitations of the bioinformatic approaches, it is suggested that they try generating substantially smaller clusters of cells in mouse (~twice the number of clusters). In this case, if none of these mouse sub-populations have a match to the human clusters based on KL divergence, that would substantially strengthen the argument that the human population is not conserved in mouse.

We carried out the requested analysis (19 clusters rather than 10) and demonstrate that with increased granularity in the mouse data, the same relationships between human and mouse cell types are seen (new Figure 3—figure supplement 1). Importantly, neither subcluster of mouse cLTMRs is a good match for any type of human DRG-neurons. Similarly, the human neuronal classes that only poorly matched the mouse clusters were not better matched to any of the subdivided mouse cell types. Nonetheless we also toned down the abstract.

https://doi.org/10.7554/eLife.71752.sa2

Article and author information

Author details

  1. Minh Q Nguyen

    National Institute of Dental and Craniofacial Research, National Institutes of Health, Bethesda, United States
    Contribution
    Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Writing – original draft, Writing – review and editing
    Competing interests
    No competing interests declared
  2. Lars J von Buchholtz

    National Institute of Dental and Craniofacial Research, National Institutes of Health, Bethesda, United States
    Contribution
    Data curation, Formal analysis, Methodology, Software, Writing – original draft, Writing – review and editing
    Competing interests
    No competing interests declared
  3. Ashlie N Reker

    Department of Anesthesiology, College of Medicine, University of Cincinnati, Cincinnati, United States
    Contribution
    Investigation
    Competing interests
    No competing interests declared
  4. Nicholas JP Ryba

    National Institute of Dental and Craniofacial Research, National Institutes of Health, Bethesda, United States
    Contribution
    Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Supervision, Writing – original draft, Writing – review and editing
    For correspondence
    nick.ryba@nih.gov
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-2060-8393
  5. Steve Davidson

    Department of Anesthesiology, College of Medicine, University of Cincinnati, Cincinnati, United States
    Contribution
    Conceptualization, Data curation, Funding acquisition, Investigation, Methodology, Writing – original draft, Writing – review and editing
    For correspondence
    davidsst@ucmail.uc.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-2944-1144

Funding

National Institute of Dental and Craniofacial Research (ZIC DE 000561)

  • Nicholas JP Ryba

National Institute of Neurological Disorders and Stroke (R01NS107356)

  • Steve Davidson

National Institutes of Health (RFNS113881)

  • Steve Davidson

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We thank the Genomics and Computational Biology Core (National Institute on Deafness and Other Communication Disorders) for sequencing; this work utilized the computational resources of the NIH HPC Biowulf cluster (http://hpc.nih.gov) and was supported in part by the Intramural Program of the National Institutes of Health, National Institute of Dental and Craniofacial Research. We are indebted to the LifeCenter, Cincinnati and the donor families for their generosity. We also thank Drs. Mark Hoon, Claire Le Pichon, and Alexander Chesler and members of our groups for valuable suggestions.

Ethics

Animal experiments were carried out in strict accordance with the US National Institutes of Health (NIH) guidelines for the care and use of laboratory animals and were approved by the National Institute of Dental and Craniofacial Research animal care and use committee (protocol #20-1041).

Senior Editor

  1. Catherine Dulac, Harvard University, United States

Reviewing Editor

  1. Rebecca Seal, University of Pittsburgh School of Medicine, United States

Publication history

  1. Received: June 28, 2021
  2. Preprint posted: July 4, 2021 (view preprint)
  3. Accepted: November 6, 2021
  4. Version of Record published: November 26, 2021 (version 1)

Copyright

This is an open-access article, free of all copyright, and may be freely reproduced, distributed, transmitted, modified, built upon, or otherwise used by anyone for any lawful purpose. The work is made available under the Creative Commons CC0 public domain dedication.

Metrics

  • 1,075
    Page views
  • 113
    Downloads
  • 0
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)

Further reading

    1. Neuroscience
    Rawan AlSubaie et al.
    Research Article Updated

    Projections from the basal amygdala (BA) to the ventral hippocampus (vH) are proposed to provide information about the rewarding or threatening nature of learned associations to support appropriate goal-directed and anxiety-like behaviour. Such behaviour occurs via the differential activity of multiple, parallel populations of pyramidal neurons in vH that project to distinct downstream targets, but the nature of BA input and how it connects with these populations is unclear. Using channelrhodopsin-2-assisted circuit mapping in mice, we show that BA input to vH consists of both excitatory and inhibitory projections. Excitatory input specifically targets BA- and nucleus accumbens-projecting vH neurons and avoids prefrontal cortex-projecting vH neurons, while inhibitory input preferentially targets BA-projecting neurons. Through this specific connectivity, BA inhibitory projections gate place-value associations by controlling the activity of nucleus accumbens-projecting vH neurons. Our results define a parallel excitatory and inhibitory projection from BA to vH that can support goal-directed behaviour.

    1. Cell Biology
    2. Neuroscience
    Angela Kim et al.
    Research Article Updated

    Insulin-induced hypoglycemia is a major treatment barrier in type-1 diabetes (T1D). Accordingly, it is important that we understand the mechanisms regulating the circulating levels of glucagon. Varying glucose over the range of concentrations that occur physiologically between the fed and fuel-deprived states (8 to 4 mM) has no significant effect on glucagon secretion in the perfused mouse pancreas or in isolated mouse islets (in vitro), and yet associates with dramatic increases in plasma glucagon. The identity of the systemic factor(s) that elevates circulating glucagon remains unknown. Here, we show that arginine-vasopressin (AVP), secreted from the posterior pituitary, stimulates glucagon secretion. Alpha-cells express high levels of the vasopressin 1b receptor (V1bR) gene (Avpr1b). Activation of AVP neurons in vivo increased circulating copeptin (the C-terminal segment of the AVP precursor peptide) and increased blood glucose; effects blocked by pharmacological antagonism of either the glucagon receptor or V1bR. AVP also mediates the stimulatory effects of hypoglycemia produced by exogenous insulin and 2-deoxy-D-glucose on glucagon secretion. We show that the A1/C1 neurons of the medulla oblongata drive AVP neuron activation in response to insulin-induced hypoglycemia. AVP injection increased cytoplasmic Ca2+ in alpha-cells (implanted into the anterior chamber of the eye) and glucagon release. Hypoglycemia also increases circulating levels of AVP/copeptin in humans and this hormone stimulates glucagon secretion from human islets. In patients with T1D, hypoglycemia failed to increase both copeptin and glucagon. These findings suggest that AVP is a physiological systemic regulator of glucagon secretion and that this mechanism becomes impaired in T1D.