Transcriptional profiling at whole population and single cell levels reveals somatosensory neuron molecular diversity

Abstract

The somatosensory nervous system is critical for the organism's ability to respond to mechanical, thermal, and nociceptive stimuli. Somatosensory neurons are functionally and anatomically diverse but their molecular profiles are not well-defined. Here, we used transcriptional profiling to analyze the detailed molecular signatures of dorsal root ganglion (DRG) sensory neurons. We used two mouse reporter lines and surface IB4 labeling to purify three major non-overlapping classes of neurons: 1) IB4(+)SNS-Cre/TdTomato(+), 2) IB4(-)SNS-Cre/TdTomato(+), and 3) Parv-Cre/TdTomato(+) cells, encompassing the majority of nociceptive, pruriceptive, and proprioceptive neurons. These neurons displayed distinct expression patterns of ion channels, transcription factors, and GPCRs. Highly parallel qRT-PCR analysis of 334 single neurons selected by membership of the three populations demonstrated further diversity, with unbiased clustering analysis identifying six distinct subgroups. These data significantly increase our knowledge of the molecular identities of known DRG populations and uncover potentially novel subsets, revealing the complexity and diversity of those neurons underlying somatosensation.

Keywords: DRG; evolutionary biology; genomics; mouse; neuroscience; nociception; peripheral nervous system; proprioception; somatosensation; transcriptome.

Figure 1.. Fluorescent characterization of SNS-Cre/TdTomato and Parv-Cre/TdTomato DRG populations.

(A) SNS-Cre/TdTomato and Parv-Cre/TdTomato lumbar DRG sections imaged for TdTomato (red), IB4-FITC, anti-CGRP, or anti-Parvalbumin (green). Scale bars, 50 μm. (B–C) Proportions of IB4⁺, CGRP⁺, NF200⁺, Parvalbumin⁺ populations expressing SNS-Cre/TdTomato or Parv-Cre/TdTomato, and converse TdTomato proportions expressing each co-stained marker (mean ± s.e.m., n = 8–20 fields from 3 animals). (D) Venn diagram depicting distinct DRG populations as labeled by Isolectin B4, NF200, and TdTomato populations. (E) For transcriptional profiling, three non-overlapping DRG populations were FACS purified: IB4⁺SNS-Cre/TdTomato⁺, IB4⁻SNS-Cre/TdTomato⁺, and Parv-Cre/TdTomato⁺ cells. DOI: http://dx.doi.org/10.7554/eLife.04660.003

Figure 1—figure supplement 1.. SNS-Cre/TdTomato and Parv-Cre/TdTomato DRG and spinal cord characterization.

(A) SNS-Cre and Parv-Cre mice were bred with Rosa26-TdTomato mice to generate lineage reporter progeny. (B) Confocal microscopy images of whole mount L4 DRG from TdTomato progeny. Scale bars, 50 μm. (C) Lumbar spinal cord sections were stained with Isolectin B4-FITC (green) and anti-CGRP (blue). SNS-Cre/TdTomato fibers overlapped with CGRP and IB4 staining in dorsal horn laminae I–II. By contrast, Parv-Cre/TdTomato fibers extended to lamina III, Clark Nucleus (C.N.) and ventral horns. (D) Lumbar sections show SNS-Cre/TdTomato fibers in lamina II (colocalized with IB4), but not lamina III stained by anti-PKC-γ. Parv-Cre/TdTomato does not innervate superficial laminae. Scale bars, 100 μm. DOI: http://dx.doi.org/10.7554/eLife.04660.004

Figure 2.. Electrophysiological properties of SNS-Cre/TdTomato and Parv-Cre/TdTomato neurons.

Whole cell current clamp recordings were conducted on SNS-Cre/TdTomato and Parv-Cre/TdTomato neurons in response to 200 pA injection. (A) Representative action potential waveforms before and after application of 500 nM TTX. (B–C) Statistical comparisons of action potential (AP) half-widths and capacitances between sensory populations (SNS-Cre/TdT⁺, n = 13; Parv-Cre/TdT⁺, n = 9; p-values by Student's t test). DOI: http://dx.doi.org/10.7554/eLife.04660.005

Figure 3.. FACS purification of distinct somatosensory neuron populations.

(A) Mouse DRG cells were stained with DAPI and subjected to flow cytometry. After gating on large cells by forward and side scatter (R1), dead cells were excluded by gating on the DAPI⁻ events; Next, TdTomato (hi) events were purified. Following purification, fluorescence and DIC microscopy show that the majority of sorted neurons are TdTomato⁺ (images on right). (B) Representative FACS plots of Parv-Cre/TdTomato⁺ and SNS-Cre/TdTomato⁺ DRG populations. Right, quantification of proportions of DAPI⁻ events in the DRG constituting each neuron population (n = 5 SNS-Cre/TdTomato mice, n = 4 Parv-Cre/TdTomato mice; p-values, Student's t test; Error bars, mean ± s.e.m.). (C) Representative FACS plot shows relative percentages of IB4-FITC surface stained and IB4⁻ neuronal populations among the total SNS-Cre/TdTomato (hi) gate. DOI: http://dx.doi.org/10.7554/eLife.04660.006

Figure 3—figure supplement 1.. Flow cytometric sorting and analysis of TdTomato⁺ neurons.

(A) By FACS analysis, TdTomato labeled both ‘high’ and ‘low’ fluorescence populations (see gates). Purified high-expressing populations corresponded to neuronal cell bodies, while the lower fluorescence consisted of fluorescent axonal debris, as shown by microscopy images post-sorting (right). (B) TdTomato neurons purified and plated onto glass slides. After 24 hr, post-sorted SNS-Cre/TdT⁺ neurons showed neurite outgrowth and relatively pure populations compared to unsorted SNS-Cre/TdT⁺ neurons. (C) Representative FACS plot overlay of light scattering properties for Parv-Cre/TdT⁺ vs SNS-Cre/TdT⁺ populations. Comparison of forward and side scatter properties on left (SNS-Cre/TdT, n = 4; Parv-Cre/TdT, n = 4; error bars, s.e.m). DOI: http://dx.doi.org/10.7554/eLife.04660.007

Figure 3—figure supplement 2.. Transcriptome analysis of purified neuronal samples relative to whole DRG tissues.

(A) Individual expression profile comparisons of sorted neuron samples (Red and green show numbers of trancripts >twofold differential expression). (B) Plots of absolute RMA normalized transcript levels for myelin associated, nociceptor associated, and proprioceptor associated genes in FACS purified SNS-Cre/TdT⁺ and Parv-Cre/TdT⁺ samples vs whole DRG samples. p-values by One-way ANOVA: ***p < 0.001. (C) Fold-change vs fold-change comparison of sorted neurons vs whole DRG datasets (red transcripts are >twofold enriched in whole DRG; blue transcripts are >twofold enriched in both sorted subsets). DOI: http://dx.doi.org/10.7554/eLife.04660.008

Figure 4.. Hierarchical clustering and principal components analysis of transcriptomes.

(A) Hierarchical clustering of sorted neuron molecular profiles (top 15% probesets by coefficient of variation), showing distinct groups of transcripts enriched in IB4⁺SNS-Cre/TdT⁺, IB4⁻SNS-Cre/TdT⁺, and Parv-Cre/TdT⁺ neuron populations. (B) Principal component analysis shows distinct transcriptome segregation for the purified populations along three principal components axes. DOI: http://dx.doi.org/10.7554/eLife.04660.010

Figure 5.. Functional somatosensory mediators show clustered gene expression across purified DRG populations.

Heat-map showing relative transcript levels for known somatosensory mediators plotted across IB4⁺SNS-Cre/TdTomato⁺, IB4⁻SNS-Cre/TdTomato⁺, and Parv-Cre/TdTomato⁺ purified neuron transcriptomes (rows show individual samples; columns are specific transcripts). Genes were grouped based on known roles linked to thermosensation/nociception, pruriception, tactile function, neurotrophic receptors, and proprioception. DOI: http://dx.doi.org/10.7554/eLife.04660.011

Figure 6.. Heat-map distribution of voltage-gated and TRP channels across neuronal subsets.

Expression patterns of different sub-types of voltage-gated ion channels and transient receptor potential (TRP) channels were hierarchically clustered and analyzed across IB4⁺SNS-Cre/TdT⁺, IB4⁻SNS-Cre/TdT⁺ and Parv-Cre/TdT⁺ purified neuron samples (columns are individual samples, heat-maps). (A) Sodium channel levels, (B) calcium channel levels, (C) potassium channel levels (top 60 differentially expressed transcripts by CoV), (D) chloride channel levels, and (E) TRP channel levels are plotted as heat-maps. For A–E, plotted transcripts show minimum expression >100 in at least one neuronal subgroup. DOI: http://dx.doi.org/10.7554/eLife.04660.012

Figure 7.. Heat-map distribution of ligand-gated ion channels, G-protein coupled receptors, and transcription factors across neuronal subsets.

(A) Expression patterns of ligand-gated ion channels, including glutamatergic, chlorinergic, HCN, P2X channels, were analyzed by hierarchical clustering (columns are individual samples). (B) Differentially expressed G-protein coupled receptors (GPCRs) were clustered and plotted across sensory subsets (Top 60 by CoV are shown). (C) Differentially expressed transcription factors were clustered and plotted across sensory subsets as a heat-map (Top 60 by CoV are shown). For A–C, plotted transcripts show minimum expression >100 in at least one neuronal subgroup. DOI: http://dx.doi.org/10.7554/eLife.04660.013

Figure 8.. Differential volcano plot analysis of SNS-Cre/TdTomato vs Parv-Cre/TdTomato transcriptomes.

(A) Pairwise comparison of SNS-Cre/TdT⁺ vs Parv-Cre/TdT⁺ profiles showing differentially expressed (DE) transcripts as a volcano plot (blue transcripts, Parv-Cre/TdT enriched; red, SNS-Cre/TdT enriched, twofold, p < 0.05). (B) Most enriched Gene ontology (GO) categories and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways in SNS-Cre/TdT vs Parv-Cre/TdT enriched transcripts, plotted as heat-map of −log (p-value). (C) Volcano plots depicting (i) calcium channels, (ii) potassium channels, and (iii) TRP channels expression differences between populations. Individual transcripts highlighted (red, SNS-Cre/TdT⁺ enriched; green, Parv-Cre/TdT⁺ enriched; blue, not significantly different: twofold, p < 0.01). DOI: http://dx.doi.org/10.7554/eLife.04660.014

Figure 9.. Differential volcano plot analysis of IB4⁺ and IB4⁻ SNS-Cre/TdTomato subset transcriptomes.

(A) Pairwise comparison of IB4⁺SNS-Cre/TdT⁺ vs IB4⁻SNS-Cre/TdT⁺ neuronal profiles show differentially expressed (DE) genes by volcano plot (blue, IB4⁺ enriched; red, IB4⁻enriched, twofold, p < 0.05). (B) Top Gene ontology (GO) categories of biological processes (BP) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways for IB4⁺SNS-Cre/TdT⁺ and IB4⁻SNS-Cre/TdT⁺ enriched transcripts, plotted as heat-maps of −log (p-value). (C) Volcano plots showing differential expression of (i) ion channels, (ii) cell adhesion molecules, and (iii) G-protein coupled receptors between neuronal populations (red, IB4⁺ enriched transcripts; green, IB4⁻ enriched; blue, not significantly different: twofold, p < 0.01). DOI: http://dx.doi.org/10.7554/eLife.04660.015

Figure 10.. Analysis of most enriched marker expression by IB4⁺, IB4⁻ SNS-Cre/TdTomato and Parv-Cre/TdTomato⁺ populations.

(A–C) Fold-change/fold-change comparisons illustrate most differentially enriched genes in each subset (highlighted in color are threefold and twofold enriched numbers). (D) Heat-maps showing relative expression of the top 40 transcripts enriched in each of the three neuronal subsets (>threefold), ranked by product of fold-change differences. DOI: http://dx.doi.org/10.7554/eLife.04660.016

Figure 10—figure supplement 1.. Fluidigm analysis of 100 and 10 cell-samples.

FACS sorted 100 cell or 10 cell samples consisting of IB4⁺SNS-Cre/TdT⁺, IB4⁻SNS-Cre/TdT⁺, and Parv-Cre/TdT⁺ neurons were analyzed by Fluidigm for 80 different transcript levels chosen based on microarray results, and normalized to Gapdh expression. Hierarchical clustering of transcript levels is shown for 100 cell and 10 cell groups as heat-maps. DOI: http://dx.doi.org/10.7554/eLife.04660.017

Figure 11.. Single cell transcript levels show log-scale distribution across neuronal populations.

Normalized transcript levels in single cells determined by parallel qRT-PCR are plotted on a log-scale comparing IB4⁺SNS-Cre/TdT⁺, IB4⁻SNS-Cre/TdT⁺, and Parv-Cre/TdT⁺ cells. (A) Nociceptor related transcript levels (Trpv1, Trpa1, Mrgprd, P2rx3, Nppb, Ptgir), (B) Proprioception related transcript levels (Pvalb, Runx3, Cdh12). Individual neurons are shown as dots in plots. DOI: http://dx.doi.org/10.7554/eLife.04660.019

Figure 12.. Hierarchical clustering analysis of single cell qRT-PCR data reveals distinct neuronal subgroups.

Heat-map of 334 single neurons and 80 genes after spearman-rank hierarchical analysis of RT-PCR data (relative gene expression normalized to gapdh). Each column represents a single sorted cell, and each transcript is shown per row. Clustering analysis finds seven distinct subgroups (I, II, III, IV, V, VI, VII). Characteristic transcript expression patterns that delineate each somatosensory subset are written below. DOI: http://dx.doi.org/10.7554/eLife.04660.020

Figure 12—figure supplement 1.. Expression of neuronal-associated transcripts across purified single cell samples by qRT-PCR.

Heat-map showing expression levels of neuron-associated transcripts across single cells (from 1 to 334) by qRT-PCR. DOI: http://dx.doi.org/10.7554/eLife.04660.021

Figure 12—figure supplement 2.. Transcript expression levels for characteristic marker genes in single cell neuron Group I and Group VII.

Plotted are normalized transcript levels of Group I and Group VII transcripts, ordered from highest to lowest expression (i.e., Grik1 to Wnt2b for Group I, Pvalb to Cdh12 for Group VII). DOI: http://dx.doi.org/10.7554/eLife.04660.022

Figure 13.. Single cell subgroups distribute differentially across originally purified populations.

(A) Principal Components Analysis of single cell transcriptional data shows distinct segregation of Groups I, V, and VII neurons. (B) Proportions of each neuronal subgroup relative to original labeled IB4⁺SNS-Cre/TdTomato⁺, IB4⁻SNS-Cre/TdTomato⁺, and Parv-Cre/TdTomato⁺ neurons. DOI: http://dx.doi.org/10.7554/eLife.04660.023

Figure 14.. Focused analysis of single cell heterogeneity and transcript enrichment in neuronal subgroups.

(A) Relative expression levels of subgroup specific transcripts in single cells across each neuronal subgroup (each bar = 1 cell). Group I (Lpar3, Mrgprd), group VI (Il31ra, Nppb), and group VII markers (Gpcr5b) show subset enrichment and highly heterogeneous expression at the single cell level. (B–C) Nearest neighbor analysis by pearson correlation of Mrgprd and Pvalb transcript levels to all 80 probes across the single cell expression dataset was generated. Correlation levels go from left to right. DOI: http://dx.doi.org/10.7554/eLife.04660.024

Figure 14—figure supplement 1.. Defining the transcriptional characteristics of Group I, II, and IV neurons.

Transcript levels for selected genes that define the characteristics of specific neuronal subgroups Group I, II, and IV neurons were plotted across all 334 individual neurons. (A) Group 1 neurons were found with high levels of P2rx3, Lpar3. (B) Group II neurons show high levels of Ntrk1 and Kcnv1. (C) Group IV are characterized by Trpv1 expression but lack of Scn10a expression. DOI: http://dx.doi.org/10.7554/eLife.04660.025

Figure 14—figure supplement 2.. Expression plots of nociceptor-associated transcripts across single cell transcriptional data.

Transcript levels for nociceptor associated genes (A) Trpa1, (B) Scn11a, and (C) Aqp1 were plotted across all individual neurons. DOI: http://dx.doi.org/10.7554/eLife.04660.026

Figure 15.. DRG subgroups I, VI, and VII characteristics defined by double RNA in situ hybridization.

(A) Double RNA in situ hybridization in SNS-Cre/TdTomato and Parv-Cre/TdTomato lumbar DRG sections for TdTomato (red) with Lpar3, Il31ra, or Gpcr5b (green), which are Group I, VI, and VII markers respectively. Lpar3 and IL31ra expression colocalize with SNS-Cre/TdTomato but not Parv-TdTomato, while Gpcr5b colocalizes with Parv-Cre/TdTomato but not SNS-Cre/TdTomato. (B) Double in situ hybridization in lumbar DRG sections for group VI marker IL31ra vs Group I marker Lpar3, Group VI marker Gpcr5b, or Group VI marker Nppb. Il31ra and Nppb in shown in a distinct subset of DRG neurons. Scale bars, 100 μm. DOI: http://dx.doi.org/10.7554/eLife.04660.028

Figure 15—figure supplement 1.. Immunofluorescence characteristics of DRG subgroup V.

(A) Expression plot shows enrichment of Th expression in Group V neurons. (B) SNS-Cre/TdTomato lumbar DRG sections were imaged for TdTomato (red), anti-TH (blue), and IB4-FITC (green). (C) Quantification of neuronal proportions TH+ neurons that are IB4⁻SNS-Cre/TdT⁺, IB4⁻SNS-Cre/TdT⁺, or Parv-Cre/TdT⁺ neurons expressing TH. Statistical analysis by Student's t test (n = 8–10 fields from 3 mice each). Scale bars, 100 μm. DOI: http://dx.doi.org/10.7554/eLife.04660.029

Figure 15—figure supplement 2.. Group I marker Prkcq is in a distinct subset of DRG neurons.

(A) Transcript levels for Prkcq plotted across all individual neuron subgroups. (B) Double in situ hybridization (ISH) of lumbar DRG sections for TdTomato (red) and for Lpar (green) shows that Prkcq⁺ neurons showed SNS-Cre/TdTomato expression whereas they were did not express SNS-Cre/TdTomato. Scale bars, 100 μm. (C) Double ISH of lumbar DRG sections shows that Prkcq does not colocalize with Group VI marker IL31ra. DOI: http://dx.doi.org/10.7554/eLife.04660.030