The mammalian somatosensory system responds to mechanical, thermal, and chemical stimuli to provide valuable sensory information about both the environment and the internal physiological state.19,47,53 The sensory neurons of the trigeminal and dorsal root ganglia (DRG) are the primary somatosensory receptors, innervating the periphery and transmitting signals to the central nervous system. Recently single cell (sc) and single nuclear (sn) sequencing have defined an array of somatosensory neuronal classes and confirmed extensive similarity between these ganglia.7,13,24,31,32,42,48 Some of these classes correspond to cells that have been previously targeted genetically as they are defined by individual functional markers.3,11,27 Such experiments support the idea that genetically defined neural classes have distinct and specific roles.17,37 However, other transcriptomic classes are only distinguished by complex patterns of gene expression.7,13,24,31,32,42,48 Moreover, these diagnostic genes tend to have less clear function in sensory detection raising questions as to the significance of transcriptomic class and their selectivity in somatosensation. A problem with sc-sequencing is that cell isolation destroys information about the in vivo system. Consequently, it has been difficult to explore how neural class is related to anatomical features of the cells including their projection patterns.
Single cell sequencing of other neuronal populations has also vastly expanded the view of their diversity by exposing a complex array of new transcriptomic classes.46,51 Just as in the somatosensory system, mapping the anatomical organization of these molecularly defined neural classes is a fundamental challenge. To address these inherent challenges, several groups have developed related in situ hybridization (ISH)-based methods6,30,41,49 that are effective at determining the cellular expression of many genes. These methods appear very powerful but require investment in probes and equipment as well as high-resolution imaging, precise alignment, and processing of massive data sets, limiting their widespread use. On the other hand, similar methods can be applied to a much smaller number of transcripts simplifying both methodology and analysis. For example, Spatial Genomic Analysis was developed to map cells in developing neural crest25 by localizing just 35 genes. Similarly, an sc-quantitative polymerase chain reaction method has been used to classify somatosensory neurons using 28 genes.1 Therefore, we reasoned that an even smaller number of highly expressed transcripts selectively marking partially overlapping sets of scRNA-defined cell classes might mean that multicolor ISH, low-resolution imaging, and just a few repeat cycles could be used to economically project transcriptomic classification of cells to anatomical sections as a means to test and extend predictions from scRNA sequencing.
Here, we developed and refined this ISH-based approach to probe how neural diversity in the somatosensory system might orchestrate sensory discrimination and elicit select behavioral responses. To simplify analysis of the data and to make this type of strategy applicable to similar problems in other systems,46,51 we developed a powerful U-Net machine learning algorithm that automatically segments cells and assigns their classes. We then used this platform to characterize critical features of trigeminal sensory neurons including determining ion-channel expression profiles and examining peripheral targeting.
2. Material and methods
2.1. Probe selection and validation
Single cell RNA sequence data32 analyzed using the Seurat package developed in the Satija lab40 were used to identify probes of multiplexed ISH. Potential probes were selected from genes that were prominent in the principle components used for clustering and were tested for their power at distinguishing trigeminal neural classes by analyzing their expression profiles in the sc-data. Criteria for probe selection included sharp boundaries in expression between classes, high level of expression in positive cells, and expression in a large proportion of cells that made up a positive class. The 15 probes chosen provided redundancy for identification of most neural classes (Fig. 1A) but only a small number of markers with potential to segregate some groups, for example, C1 from C2 or C9 from C10 could be identified.
2.2. Mice and retrograde tracing
Animal experiments were performed in strict accordance with the US National Institutes of Health guidelines for the care and use of laboratory animals and were approved by the NIDCR or NINDS ACUCs. C57BL/6N and FVB/N mice were purchased from Harlan. Mas-related G-protein–coupled receptor member A3 (Mrgpra3)-Cre mice17 were crossed to Ai14 tdTomato reporter mice (Jackson Labs). Both male and female mice were used for experiments.
Wheat germ agglutinin (WGA) coupled to Alexa594 (Invitrogen, Carlsbad, CA) was injected into various target fields to retrogradely label trigeminal neurons from their terminal processes. All labeling procedures used WGA-Alexa594 at 1 mg/mL; mice were anesthetized with isoflurane. To label the eye, the cornea was first abraded with a beveled 30-g syringe needle; 0.5- to 1-μL WGA-Alexa594 was applied. To label the meninges, a midline incision into the skin allowed exposure of the skull, six 1-mm holes were drilled into the dorsal skull without perforating the meninges. Wells were built around each craniotomy with dental cement; 2 × 0.5-μL WGA-Alexa594 was applied to each well. Craniotomies were sealed with coverslips and dental cement, and the skin incision was closed and sutured. In all cases, animals were euthanized 16 to 20 hours after WGA-Alexa594 application and trigeminal ganglia dissected.
WGA-Alexa594 fluorescence was imaged in ISH-imaging buffer using 594-nm laser stimulation and a 9-nm wavelength scan, which helped separate signal from broader autofluorescence of the tissue. Slides were washed 3 times in 2 × saline sodium citrate (SSC) and processed further for ISH.
2.3. In situ hybridization
Hybridization chain reaction (HCR) version 3.08 was used for all ISH. Buffers and probes against the mouse genes Cd34, Sstr2, Trpv1, Etv1, Tmem233, S100b, Synpr, Tmem45b, Calca, Trpa1, Mrgprd, Trpm8, Nppb, Fxyd2, Nefh, Piezo2, Tubb3, sodium channel protein type 1 subunit alpha (Scn1a), sodium channel protein type 8 subunit alpha (Scn8a), sodium channel protein type 10 subunit alpha (Scn10a), and sodium channel protein type 11 subunit alpha (Scn11a) and against Cre recombinase and tdTomato (tdT) were purchased from Molecular Instruments (Los Angeles, CA) as a mix of 10 or more oligonucleotide pairs per probe.
Trigeminal ganglia dissected from 2- to 6-month-old C57BL/6N or FVB/N mice were embedded in optimal cutting temperature medium and frozen at −80°C. Approximately 20-μm sections were cut along a horizontal axis through the ganglion and mounted on microscope slides. Sections were dried (1 hour on a 37°C hotplate), were fixed in on ice for 15 minutes 4% paraformaldehyde (Electron Microscopy Sciences, Hatfield, PA) in phosphate-buffered saline (PBS), and then were washed 6 times in PBS. Slides were dehydrated in an ethanol series (50%, 70%, 100%, and 100% for 5 minutes each) and stored for up to 1 week. This preparation resulted in excellent retention of sections through multiple rounds of hybridization and detection of probes but was not specifically optimized for preservation of signal.
Slides were dried completely, washed 3 times in PBS, and prehybridized for 30 minutes at 37°C in HCR hybridization solution (Molecular Instruments). Hybridization and amplification were performed as described previously8 with minor modifications. Hybridization typically used 5 probes with adapters B1 to B5 at a concentration of 4 nM per probe in Coverwell incubation chambers (Grace Biolabs, Bend, OR) and was carried out for 48 to 72 hours at 37°C. Washing and preparation of amplifier hairpins for adapters B1 to B5 conjugated to Alexa488, Alexa514, Alexa561, Alexa594, and Alexa647 was as described by the manufacturer. Amplification was performed in a Coverwell chamber for 12 to 16 hours at room temperature in the dark. Slides were washed 2 times for 15 minutes in 5 × SSC containing 0.1% Tween 20 and coverslipped in Imaging Buffer (3-U/mL pyranose oxidase, 0.8% D-glucose, 2 × SSC, 10-mM Tris HCl pH 7.4).
After imaging, sections were rinsed in 2 × SSC, and DNA probes and amplifiers were removed by incubation in 250-U/mL RNase-free DNase (Roche Diagnostics, Indianapolis, IN) for 90 minutes at room temperature. Sections were washed 6 times (5 minutes in 2 × SSC) and prehybridized for the next round of hybridization with the next 5 probes. This procedure was repeated for up to 4 rounds of hybridization.
2.4. Confocal imaging and signal unmixing
Imaging was performed using a Zeiss LSM 880 confocal microscope with spectral detector using a 10×/0.45NA air objective with the pinhole opened to 132.9 μm to allow capture of the whole section and provide enough signal intensity at this comparatively low magnification. Alexa647 was stimulated with a 633 laser and emission was collected from 644 to 752 nm. Alexa546 and Alexa594 were stimulated with 561 and 594 lasers, respectively, and wavelength scans with 9-nm windows were collected to allow unmixing signals from different fluorophores. To resolve Alexa488 from Alexa514, sequential wavelength scans with 3-nm windows were performed while stimulating with a 488-nm laser. Spectral unmixing was performed using Zeiss ZEN software.
Image series from consecutive hybridizations were combined by manually aligning transmitted light images from each imaging session using the TurboReg plugin in ImageJ/Fiji.
2.5. Unsupervised clustering of manually outlined cells
Data processing including the supervised vs unsupervised steps is schematically described in Figure S1 (available at http://links.lww.com/PAIN/B13). As a starting point, cell regions of interest (ROIs) in 4 complete ganglion sections were manually outlined in ImageJ (based on the individual ISH images), and mean fluorescence for each cell and probe was extracted using custom ImageJ macroscripts. Further processing was performed using Python scripts: For subsequent clustering, local background in a ring around the cell was subtracted and the resulting values were standardized to mean 0 and variance 1 for each probe within a whole tissue section typically containing 1500 to 2000 cells.
The dimensionality of the data was then reduced from the 13 relevant probes (Cd34, Trpv1, Etv1, Tmem233, S100b, Synpr, Tmem45b, Calca, Trpa1, Mrgprd, Trpm8, Nppb, and Fxyd2) to a 2-dimensional representation by t-distributed stochastic neighbor embedding (tSNE). Cells clusters were identified using a density-based algorithm (https://github.com/alexandreday/fast_density_clustering) and mapped to the clusters derived from sc-sequencing based on their average expression of marker genes (see Fig. S1 for details, available at http://links.lww.com/PAIN/B13).
2.6. Automatic detection of cell classes using supervised Deep Learning
To automatically detect cell classes in tissue section images, we performed semantic segmentation by designing a custom Fully Convolutional Neural Network that follows the popular U-Net architecture.38 Several network designs and hyperparameters such as numbers of layers, learning rate, dropout rate, final activation function, and loss function were evaluated before settling on the final model.
The multichannel ISH input image data were first tiled into 256 × 256 images, standardized to mean 0 and variance 1 for each probe channel and divided into training (96 images) and validation data (48 images). The desired output masks were generated as multihot encoded images (11 channels for 10 classes and background) from the hand-annotated ROIs and the cell class assignment from unsupervised clustering. To improve training performance, the available input data and their corresponding cell class masks were augmented by rotation and flipping. The final model closely mirrors the original U-Net architecture and contains 5 downward convolutional blocks with 2 layers each and 4 corresponding symmetric blocks for upconvolution. Image resolution was reduced to half in each consecutive block, and the number of channels was doubled starting from 16 in the first block. The output was a 256 × 256 × 11 matrix calculated by a softmax activation layer. The categorical cross-entropy loss function was optimized with an Adam optimizer at a learning rate of 0.0005 to convergence on the validation loss (∼120 epochs).
To achieve good class separation with the minimal number of hybridizations, we trained a model that takes 8 channels (Tmem233, S100b, Calca, Trpa1, Mrgprd, Trpm8, Nppb, and Fxyd2) as an input. This model was validated and used to predict cells and their classes in WGA tracing experiments. To determine cell class identity of manually outlined WGA-labeled cells, the mean probability for each cell class was calculated for the ROI and the most probable cell class was assigned to the cell. Cells with a maximum class probability of less than 0.1 were kept unassigned. Rare retrogradely labeled cell classes (<5% from any target site) were manually inspected and corrected on a cell-by-cell basis.
To automatically segment cells, probability predictions for each class were first blurred with a Gaussian filter (3 × 3 kernel) and discretized by choosing the most probable class for each pixel. After morphological opening of the foreground classes with a 6 × 6 circular kernel, cells were segmented by a distance transform watershed algorithm.
Expression level of sodium channels was determined automatically by measuring signal intensity in segmented and classified cells. For comparison purposes, violin plots of channel expression in sc- and sn-data31,32 were generated. This representation is particularly useful for analyzing sparse transcriptomic data because it eliminates extremes (primarily from the zero values) thereby providing a valuable representation of gene expression across classes. The ISH data are not of this form, with no true zero value; we therefore calculated the relative standardized expression levels (mean = 0 and variance = 1) and used a histogram representation to approximate the violin plot of the transcriptomic data.
All computations were performed either in Python using Keras/Tensorflow, sklearn, Numpy, Pandas, skimage, opencv, and scipy libraries or in custom ImageJ macroscripts. All code is available in Github (https://github.com/lars-von-buchholtz/InSituClassification).
3.1. Development of an in situ hybridization approach for classifying trigeminal neurons
scRNA sequencing of trigeminal neurons defined about a dozen neuronal classes that closely matched classes identified by transcriptomic analysis of DRG neurons.32 More recent sn-analysis31 distinguished most of the same cell populations but further subdivided classes of putative large diameter neurons, which were much more prominently represented in this study. However, some of the diversity amongst nociceptors was not observed in the nuclei-based sequencing. Thus, only 10 of the 13 classes directly corresponded between studies. Here, we developed a simple, robust approach for mapping trigeminal neural classes to cells in tissue sections to expose their anatomical features including their distribution and projection targets as well as test predictions from these transcriptomic studies.
As a starting point, we performed quantitative ISH using an array of 15 genes that we predicted should, in combination, identify all the sc-sequence defined classes. The chosen marker genes featured prominently in the dimensional reduction (principal components) used in the scRNA-sequence cluster analysis. In the sc-data, each of these genes is present in a distinct subset of trigeminal neural classes. Importantly, the chosen transcripts were generally very good markers of the clusters where they were present, that is, because they were highly expressed genes, we were confident of their presence in the vast majority of cells, whereas they were essentially absent from other types of neuron (Fig. 1A). Thus, the combinatorial expression pattern of these transcripts would be expected to allow trigeminal neural classes to be reliably distinguished, providing a framework for classifying cells in situ and a straightforward approach for testing predictions from the transcriptomic analyses.
Using the well-characterized HCR protocol,8 we performed 3 rounds of 5 probe ISH and used spectral unmixing to resolve specific signals from the 15 diagnostic probes (Fig. 1B). The HCR approach has the sensitivity and resolution to detect single transcripts when sections are imaged at high resolution.8,41 Importantly, however, for studying genes that are highly expressed in subsets of neurons, we could use lower magnification imaging, saving time, and reducing the complexity of data processing as well as maximizing the number of cells analyzed. Such integrated HCR signals have been shown to accurately quantify the relative expression level of a transcript between cells8 and thus served to quantitate the diagnostic genes at a cellular level. It should be noted, however, that higher resolution imaging is also possible and may be useful for genes with lower expression level.
Fourteen of the 15 probes clearly differentiated positive from negative cells, but the probe to Sstr2 was not diagnostic at this magnification (Fig. 1B), thus clusters C9 and C10 cannot be resolved. These clusters also failed to segregate in sn-analysis,31 and few other distinguishing markers could be identified from the sc-data. Therefore, it is likely that C9 and C10 are broadly similar types of cells. In the future, higher resolution imaging and use of additional marker genes as well as increasing the gene coverage of the Sstr2-probe set may help resolve these 2 classes of neurons, which can be distinguished in transcriptomic analysis of DRG neurons.24,42,48
In situ hybridization images from the 14 diagnostic probes could be aligned by simple translation and rotation (Fig. S2A, http://links.lww.com/PAIN/B13). After alignment, the ISH data (Fig. 1B) are equivalent to 91 different double-label experiments (Movie S1, available at http://links.lww.com/PAIN/B14) or 364 combinations of 3 markers applied to a single section. At the most basic analytical level, almost all the pairwise expression patterns that were predicted (Fig. 1A) were detected by ISH (Movie S1, available at http://links.lww.com/PAIN/B14). However, the expected division of Trpm8-expressing cells into groups expressing or not expressing Nefh (Fig. 1A) was graded rather than binary (positive vs negative, see Fig. S3A, http://links.lww.com/PAIN/B13). Although we might be able to divide 2 classes of Trpm8-neurons by thresholding Nefh-expression, this seemed arbitrary and subject to error. Thus, the Nefh-based ISH approach does not reliably distinguish classes C1 and C2. Interestingly, although this split was predicted by the cell-based analysis,32 it was not seen in nuclei-based sequencing.31 A second distinguishing marker that was predicted to divide the C1 and C2 clusters in the sc-analysis32 encodes the neuropeptide galanin (Gal). We therefore tested whether Gal expression subdivides the Trpm8-positive neurons into 2 groups using ISH (Fig. S3B, http://links.lww.com/PAIN/B13). Our results indicated that Gal is not coexpressed with Trpm8 (Fig. S3B, http://links.lww.com/PAIN/B13), again matching sn-data. Therefore, the ISH data indicate that variation amongst Trpm8 cells is more continuous that we had originally predicted32 supporting more recent results,31 which indicate that these neurons although not homogeneous, in fact comprise a single transcriptomic class with graded expression differences between the cells.
3.2. Classifying trigeminal neurons in manually segmented sections defines transcriptomic class
Trigeminal sections from 4 mice were analyzed by multiplexed ISH, and neurons were manually segmented to measure the cellular expression of marker genes. For 3 of these animals, the nondiagnostic probes for Nefh and Sstr2 were replaced with probes for other genes including neuronal β-tubulin (Tubb3, a pan-neuronal marker), providing a means to identify all somatosensory neurons in the sections. After alignment of images, all neurons in the sections were manually segmented using the signal in the 15 different ISH channels as a guide. The cellular expression levels of the 13 diagnostic markers in these manually segmented cells were determined based on signal intensity (see methods for details). The resulting 13-dimensional expression data were subjected to tSNE followed by density-based cluster identification to classify neurons (Fig. 2A). The relative expression of each marker gene at a cellular level in the classified neurons is shown (Fig. 2B) highlighting a close match with predictions (Fig. 1B) for 10 trigeminal neuronal classes. Importantly, these classes of neurons match those that were detected in both the cell and nuclear-based analyses.31 Nonetheless, cells with an expression profile corresponding to an 11th class C7 that was only identified in the sc-study,32 expressing Trpv1 and Tmem233 but not Tmem45a (Fig. 1B), were also found using ISH (Fig. S3C, http://links.lww.com/PAIN/B13). However, these neurons were rare and did not separate from C9/10 cells in the cluster analysis (Fig. 2).
Quantitative analysis of ISH based clustering revealed a high degree of consistency in representation of neural classes between sections from different mice (Fig. S4, http://links.lww.com/PAIN/B13). Therefore, we are confident that distribution of neural classes (Fig. 2C) is representative of the ganglion as a whole. Notably, the combined ISH patterns of the 13 diagnostic probes and the expression of the pan-neuronal marker Tubb3 were almost identical (Fig. S5A, http://links.lww.com/PAIN/B13), confirming that the vast majority of trigeminal neurons were included in our classification. Thus, the multigene hybridization approach demonstrates the existence of at least 11 of the 13 types of trigeminal neuron predicted by scRNA sequencing32 and reveals their true representation in the ganglion.
Although primarily designed to identify the transcriptomic class of neurons in anatomical sections, the redundancy inherent in the probing strategy (Fig. 1A) meant that it could also be used to search for rare unexpected patterns of coexpression. By contrast, to sc-transcriptomic data where dropout and contamination generates ambiguity, the ISH approach allows an explicit check by visual examination. For instance, a small proportion of C6 neurons expressed the heat- and capsaicin-activated ion-channel Trpv1 (Fig. S5C, http://links.lww.com/PAIN/B13). This expression profile likely identifies this group of neurons as the Aδ-thermal nociceptors that have been reported to trigger rapid, Trpv1-dependent, withdrawal from noxious heat.28 Importantly, these Trpv1, Calca, S100b-positive cells that were not directly identified in sc-analyses7,13,24,31,32,42,48 are predicted to play a distinct functional role from other C6 neurons.
3.3. A machine learning–based approach for segmenting neurons and assigning cell class
Although the manual segmentation approach (Fig. 2) allows for classification of cells in tissue sections, it is slow and defining the cell boundaries is often somewhat subjective. Therefore, we reasoned that automation would greatly increase throughput, decrease bias, and thus provide a versatile platform for answering questions about gene expression and neuronal connectivity in the trigeminal ganglion.
Ideally, we wanted to develop a procedure that could take unprocessed aligned ISH images, define cell boundaries, and assign transcriptomic class without relying on additional input about cell shape. Unfortunately, traditional methods used in automated segmentation35 were not suitable for the ISH images where signal had variable intensity across positive cells and was often overlapping between neighboring cells as well as regionally localized with perinuclear vs cytoplasmic staining for different probes (Fig. 1B). However, since each ISH image contains both information about cell structure as well as gene expression, we predicted that the combined images might define both cell boundaries and neuronal transcriptomic class. Therefore, instead of segmenting individual cells, we set out to develop a custom semisupervised Deep Learning algorithm using a fully convolutional Neural Net based on the U-Net architecture.38 This approach uses the aligned individual ISH images as input and calculates the probability of each pixel being one of the 10 major classes (or non-neuronal, ie, not one of these classes; see Fig. S1 for a schematic description of all data-processing, available at http://links.lww.com/PAIN/B13). The Neural Net was trained on the manually outlined (supervised) cells that were assigned a class based on the unsupervised clustering described above (see Methods for details). Although we began by using the full data set, we realized that just a subset of the 13 genes (Trpm8, S100b, Fxyd2, Calca, Trpa1, Nppb, Tmem233, and Mrgprd, see Fig. 1A, Fig. S1, available at http://links.lww.com/PAIN/B13) should be sufficient to define the 10 neural classes that were distinguished in the manually segmented cells (Fig. 2) and here report results based on this 8-gene panel.
Just as we had hoped, the Neural Net transformed the ISH data of this 8-gene probe-set into images where cells and their classes were immediately recognizable (Fig. 3A). Indeed, cells were now easily segmented (Fig. 3B) using a standard watershed approach. Some of the ISH data, including the images shown (Fig. 3), had been manually segmented (revealing the “ground truth”) but had not been used to train the U-Net algorithm. The reliable segmentation and classification of neurons in these test sections of the trigeminal ganglion (Fig. 3B and Table 1) validate the automated Deep Learning approach; at a pixel level, the prediction displayed an accuracy of approx. 95%. One minor concern is that the neural net approach led to some distortion of cell size and shape that was most pronounced in areas where neurons were densely packed and ISH signals partially overlapped. This makes it difficult to assess how cell diameter varies across the various classes but does not alter conclusions about cellular identity (Table 1). Importantly, as demonstrated by comparing the combined predictions with Tubb3 expression (Fig. S5B, http://links.lww.com/PAIN/B13), we found that almost all neurons were classified by the U-Net algorithm. Equally importantly, non-neuronal regions were also accurately reported. Notably, minor modifications to the code would make this U-Net based strategy generally applicable for analyzing other types of complex image-based data sets. Thus, we envision that analogous multiplexed ISH or antibody-based approaches coupled with automatic classification of cells could be easily adapted to other sc data sets and serve as a simple platform for mapping cell class to anatomical sections.
3.4. Expression of voltage-gated sodium channels in trigeminal neurons classes
sc-transcriptomic analyses make use of many variable genes in clustering data and thus can tolerate the sparse nature of the underlying data sets, the consequent dropout of specific transcripts as well as occasional capture of false positives. However, the expression of any individual gene will be distorted by these artefacts. By contrast, ISH more accurately reflects the cellular level of mRNAs, and thus, the multiplexed approach should provide a platform for measuring cellular expression of genes across cell classes. Here, we illustrate the power of this approach by testing the expression profile of several voltage-gated sodium channels that are known to be differentially expressed in somatosensory neurons52 and play distinct roles in sensory detection.10,34,54
The patterns of expression of Scn10a (Nav1.8, a tetrodotoxin [Ttx] insensitive channel) and Scn1a (a Ttx sensitive channel Nav1.1) were very similar across the major neural classes in the sc- and sn-analyses (Fig. 4A). Expression of Scn10a in U-Net classified trigeminal neurons very closely matched the overall predictions from the sc-data (Fig. 4B). Notably using the ISH approach, each class of sensory neurons appeared quite homogeneous in its expression pattern closely following a unimodal distribution, but there were major differences in the mean expression level of the classes (Fig. 4B). By contrast, the sequence data were more ambiguous. This was particularly apparent for C3 or C7 to 10 cells where discontinuous violin plots (Fig. 4A) could be interpreted as only a subset of cells expressing this gene, or alternately as dropout of a moderately expressed transcript because of the sparse nature of the sequence information. The ISH indicates that the majority of neurons in these classes express Scn10a but that the average level of expression is low (Fig. 4B). Importantly, this sodium channel, which plays a crucial role in detecting painful input, is absent from the C1/2 cooling sensing neurons as well as the large diameter C4 and C5 putative mechanosensors but is present albeit at very different levels in all other classes.
The classes of neurons expressing Scn1a in the 3 data sets (sc- and sn-transcriptomic and ISH based) were also highly concordant (Fig. 4), emphasizing the effectiveness of the minimal U-Net approach at projecting the transcriptomic trigeminal neuronal class onto anatomical sections. Again, just as for Scn10a, the expression of Scn1a appeared homogeneous in most (8 of the 10) neuronal classes (Fig. 4B). This contrasted starkly with transcriptomic data where even the classes with the highest overall expression showed major dropout presumably because of the low to moderate expression of this gene resulted in false negatives in the transcriptomic data. Thus, these ISH gene expression data (Fig. 4B) not only validate sc-predictions but also extend them. C1/2, C4, and C5 neurons all express this sodium channel at a relatively high level with less expression detected in C3 neurons. C7/9/10, C11, C12, and C13 neurons were negative for this channel, while expression in C6 and C8 cells appeared heterogeneous with some neurons in these classes being negative and others strongly positive.
Interestingly, other voltage-gated sodium channels (Scn11a and Scn8a) were less consistent between the 2 transcriptomic data sets (Fig. 4A), probably reflecting their expression level, the depth of sequencing in the respective studies, instability of larger diameter neurons for the cell-based analysis, and perhaps differences in how these mRNAs are distributed in trigeminal neurons or regulated during cell isolation. By contrast, ISH analysis of Scn11a was diagnostic with each class exhibiting quite a homogeneous pattern of expression (Fig. 4B). The classes of neurons expressing this gene closely matched the expression profile of Scn10a except that Scn11a (a second Ttx-resistant channel, ie, selectively expressed in sensory neurons and involved in pain sensation) was not detected in the majority of C6 neurons that are predicted to represent Aδ-nociceptors.2,14,42 Our ISH analysis revealed that Scn8a is primarily expressed in large diameter neurons with highest levels in C4, C5, and C6 cells (Fig. 4B), is also detectable in C1/2, C3 and C12 neurons, and is essentially absent from other trigeminal neuronal classes. Taken together, these data further validate the minimal ISH approach for classifying trigeminal neurons in anatomical sections and reveal its power for quantitative expression analysis by exposing the stereotyped, strongly class-related, differential expression of Ttx-sensitive and -insensitive sodium channels in the trigeminal neuronal classes.
3.5. Classifying trigeminal neurons innervating specialized sensory targets
The trigeminal nerves innervate diverse tissues and structures in the head and neck. Multiplex ISH-based neural classification provides an opportunity to probe the transcriptomic identity of neurons innervating quite different types of specialized sensory environments. Here, we illustrate this using a retrograde tracer (fluorescently labeled WGA) to mark neurons projecting to the surface of the eye or the meninges (Fig. 5A). After one round of imaging to identify the labeled neurons, the cells projecting to these targets were classified using the 8-probe U-Net approach (Fig. 5B). Our data, using a comprehensive and unbiased approach extend previous studies that have used candidate transgenes to label subsets of neurons innervating tissues.5,18,36
An immediate surprise was that tracer applied to the surface of the eye labeled most classes of sensory neurons (Fig. 5C). However, the representation of these classes was dramatically different from the ganglion as a whole (Fig. 5C). The proportion of C8 cells (Trpv1/Trpa1 expressing nociceptors) innervating the eye was increased more than 10-fold. Based on their expression profile, including prominent expression of Trpa1, we expect C8 neurons to respond to irritants including the lachrymatory agent from onions, a volatile electrophile that activates Trpa1.26 C8 neurons also express Calca, consistent with previous studies demonstrating that CGRP-neurons innervate the cornea.5 Other peptidergic c-fibers C7/9/10 (Trpv1-positive, Trpa1-negative nociceptors) were not significantly different from their distribution in the whole ganglion, whereas C6 (CGRP-expressing Aδ-nociceptors) was slightly underrepresented. A second class of cells that were particularly prominent and overrepresented in neurons targeting the eye were the C1/2 (cool responsive) cells (Fig. 5C), supporting previous studies showing the cornea to be strongly innervated by Trpm8-expressing cells.18,36 By contrast, other types of sensory neurons rarely innervated the eye including C3 (putative c-fiber low threshold mechanoreceptors, LTMRs) and the itch related class (C11). Since we did not exclusively target the cornea with our injections, some of the labeled neurons may innervate regions around the eye including the conjunctiva, which, unlike the cornea is known to be rich in Mrgpra3 (C12) fibers.18 Major mechanosensory classes C4 (A-type LTMRs) and C13 (Mrgprd-expressing mechanonociceptors) were labeled at about 50% the expected frequency based on their distribution in the ganglion. Given the overall prominence of C4 and C13 neurons in the ganglion, this modest reduction still means that many putative mechanosensors target the eye. Since fast conducting mechanosensory neurons often make specialized connections in hairy and glabrous skin and exhibit select functions in gentle touch, it will be interesting to determine the terminal structure of these neurons and to assess what contribution they provide in corneal sensation.
Our results revealed that the meninges are also targeted by a broad mix of neurons but with a very different representation of classes (Fig. 5C). Specifically, C1/2 neurons, which prominently innervate the eye, rarely targeted the meninges and C8 neurons were completely absent. Remarkably the 2 classes of itch neurons (C11 and C12) prominently targeted this internal sensory environment (Fig. 5C). Initially, we were concerned that problems with the model and class assignment might be responsible for suggesting itch neurons project to a sensory environment that cannot be scratched. However, manual examination of positive neurons revealed that they were not misassigned by the model, with prominent expression of Nppb, a very selective marker for C1132,48 in WGA-labeled cells of this class (Fig. 6A). Since the classification of C12 cells in the eight-probe model relies on overlapping expression and absence of certain markers, we performed additional retrograde labeling experiments to specifically confirm that C12 (Mrgpra3-positive) neurons also project to the meninges. To do this, we used a well-characterized Mrgpra3-Cre transgenic mouse line17 to drive reporter gene expression and demonstrated that Cre and the reporter gene were coexpressed in neurons retrogradely labeled from the meninges (Fig. 6b). Thus, both classes of sensory neurons that selectively trigger itch responses when stimulated through the skin also target the meninges.
The biggest differences in the innervation of the eye vs the meninges were the overrepresentation of C1/2 and C8 cells targeting the eye vs C11 and C12 cells projecting to the meninges. Intriguingly, the cell bodies of these 4 neuronal classes were differentially distributed within the trigeminal ganglion with this spatial segregation of classes conserved between animals (Fig. S6A, http://links.lww.com/PAIN/B13). Similar segregation of these classes extended to areas that do not target the eye and meninges (Fig. S6B, http://links.lww.com/PAIN/B13), perhaps hinting that other target tissues share the same differential innervation patterns. However, although it is attractive to speculate that these differences in spatial representation reflect differences in peripheral targeting, the eye and meninges were represented by broad and partially overlapping fields. Thus, it is also possible that developmental processes (rather than projection per se) result in this topographic segregation of neuronal classes in the ganglion.
C3 and C4 neurons (c- and A-type LTMRs) were underrepresented in meningeal projections as were C13 neurons (Mrgprd-expressing, nonpeptidergic nociceptors) but still account for about a third of all neurons targeting this tissue. It will be interesting to determine the types of stimuli these putative mechanosensors respond to, examine their sensory termini in the meninges, and assess if they play a specific role in types of migraine. Notably, the Calca-positive C6 and C7/9/10 (A- and c-type) peptidergic nociceptors were more than 2-fold enriched amongst the meningeal-labeled cells. Interestingly, we also discovered that all the meningeal C12 neurons (Fig. 6) expressed Calca, although this gene is normally present in just a fraction of the C12 cells in the whole ganglion (Fig. S7, http://links.lww.com/PAIN/B13). Thus, almost half of the trigeminal neurons innervating the meninges are positive for Calca and consequently express the neuropeptide CGRP, which is well known to play an important role in certain types of migraine headache.15,43
Here, we set out to develop a simple but robust platform for assigning transcriptomic class to trigeminal ganglion neurons in anatomical sections. Just 8 probes were required to identify all the major classes of neurons that are consistently defined by cell- and nuclear-based scRNA sequence analyses.31,32 Moreover, the hybridization images from these 8 probes were all that was needed for automated and rapid segmentation of cells and class assignment using a novel U-Net based approach that we developed. We demonstrate the power of this simple platform by determining the expression of functionally relevant sodium channels in the different classes of neuron and analyzing the neural classes innervating the eye and meninges. This approach provides a much more consistent view of gene expression, at a cellular level, than sc-transcriptomic analysis (Fig. 4) and is particularly valuable for moderately expressed genes where dropout in sc-sequencing is a major drawback. It also provides a rapid and unbiased means for linking anatomy to gene expression (Fig. 5).
A surprise from our study of tissue innervation was that transcriptomic class did not appear to be the sole predictor of a neuron's sensory role. For example, in the skin, C11 and C12 neurons selectively respond to pruritogens and trigger itch.4,23 We found that these cells also innervate the meninges where they must drive different perceptual and behavioral responses. One possibility is that these types of neurons have common sentinel-type roles, for example, as detectors of immune cell activity44 or irritants in the different tissues. In this scenario, their differential central connectivity would then govern the consequences of their activation. Interestingly, it was recently shown that the lung is also innervated by C12 neurons,16 thus the meninges may not be unique in “repurposing” itch-related cells. Analogously, the innervation of the surface of the eye by C1/2 cool sensing neurons37 is most likely not related to the cornea being an important structure for sensing environmental temperature. Instead, the cooling effect of evaporation is probably the main sensory cue, acting as a selective driver for lacrimal-gland stimulation.36 Analogously, corneal C8 neurons (Trpv1/Trpa1 positive) probably play a specialized protective role by driving tear production in response to irritants, something familiar to most people as a rather painless but nonetheless dramatic side effect of chopping onions. Thus, a sensory neuron's gene-expression profile provides just one dimension of its identity rather than a direct indication of its sensory role. Our results strongly suggest that the significantly divergent patterns of connectivity amongst transcriptomically related neurons influence both their sensitivity and output. These conclusions about the functional flexibility of well-defined transcriptomic cell types are likely to be more widely applicable and not specific to peripheral sensory neurons.
Innervation of the meninges is particularly relevant in the context of headache and migraine. Thus, it was interesting that half the neurons targeting the meninges were from Calca (CGRP)-positive classes in keeping with the role of this neuropeptide in certain types of migraine.15,43 Perhaps surprisingly, other cell types expressing trigeminal genes with reported roles in driving headache, for example, Trpm812 (in C1/2 neurons) and Trpa121 (in C8 cells) minimally targeted the meninges (Fig. 5). For C1/2 neurons, this is in line with previous reports of only very sparse innervation of the dura but does suggest that this gene that has been genetically linked to migraine in humans12 may exert its effects at other target sites. Similarly, recent results suggest Trpa1 agonists act at nonmeningeal sites to induce migraine.22,45 For example, selective ablation of Trpv1-expressing trigeminal neurons projecting to the nasal cavity blocked Trpa1-agonist–induced changes in meningeal blood flow22 that are believed important in migraine induction. Moreover, Trpa1 was recently reported to be expressed in cerebral artery endothelial cells and to drive changes in blood flow to the brain in response to Trpa1 activation completely independent of the trigeminal system.45 Since our data demonstrate that C8 neurons do not innervate the meninges, these mechanisms22,45 likely dominate the migraine-inducing effects of compounds such as mustard oil that potently activate Trpa1.
The concept of neural class has been completely altered by recent large-scale scRNA sequencing efforts that have exposed new levels of transcriptomic diversity.46,51 Although there is evidence that at least some of these groupings delineate differences between cells, it is unclear how homogeneous each of these classes is.33 Just as importantly, many of these new transcriptomic classes share extensive similarity with each other and subdivide known types of neurons without obvious functional implication. Mapping class to cells in tissues will likely help address these issues and will provide important insight into neural diversity. Several groups have reported powerful approaches to identify many different genetic targets in anatomical sections that could greatly extend the power of scRNA sequencing.6,30,41,49 However, these methods are generally complex, generate extremely large amounts of data, and consequently represent a major undertaking that may limit their use. The basic approach that we describe here for classification of trigeminal neurons represents a compromise between these technically challenging methods and more standard ISH that is often anecdotal and limited to just identifying predicted coexpression patterns. Even if full characterization of the myriad new transcriptomic classes in tissues as complex as the cortex would require localization of more genes than the trigeminal ganglion, the approach we have taken is easily scalable. Moreover, careful selection of probes should make it possible to specifically focus on the most relevant cell types and to simplify the experimental platform. We envision that combining this type of ISH approach with viral tracing50 will be particularly useful for exposing rules guiding long-range connectivity. Similarly, determining the transcriptomic class of functionally labeled neurons9,29,39 could be highly informative.
Here, we demonstrate how defining transcriptomic class of trigeminal neurons in situ reveals important aspects of their role in sensation. This approach should be easy to extend to DRG neurons simply by adding one additional gene (eg, parvalbumin, Pvalb) to classify proprioceptors. In the future, this type of neuronal classification will be especially useful for examining changes in gene expression induced by pathological conditions that trigger chronic pain.20 Fully annotated code is available in Github (https://github.com/lars-von-buchholtz/InSituClassification).
Conflict of interest statement
The authors have no conflicts of interest to declare.
Appendix A. Supplemental digital content
Supplemental digital content associated with this article can be found online at http://links.lww.com/PAIN/B13 and http://links.lww.com/PAIN/B14.
Supplemental video content
A video abstract associated with this article can be found at, available at http://links.lww.com/PAIN/B15.
The authors are grateful to Dr. Xinzhong Dong for providing Mrgpra3-Cre reporter mice, to the NINDS Light Imaging Facility, and to other members of our labs for help and advice.
Supported by the Intramural program of the National Institutes of Health, National Institute of Dental and Craniofacial Research (N.J.P.R.) and National Center for Integrative and Complementary Health (A.T.C.) and included funding from Department of Defense in the Center for Neuroscience and Regenerative Medicine (A.T.C.).
Author contributions: study design: L.J. von Buchholtz, A.T. Chesler, and N.J.P. Ryba; experimentation: L.J. von Buchholtz, J.J. Emrick, and R.M. Lam; all software development: L.J. von Buchholtz; data analysis: L.J. von Buchholtz, J.J. Emrick, A.T. Chesler, and N.J.P. Ryba; manuscript draft: L.J. von Buchholtz, A.T. Chesler, and N.J.P. Ryba; final manuscript: L.J. von Buchholtz, J.J. Emrick, R.M. Lam, A.T. Chesler, and N.J.P. Ryba.
Fully annotated code is available in Github: https://github.com/lars-von-buchholtz/InSituClassification.
. Adelman PC, Baumbauer KM, Friedman R, Shah M, Wright M, Young E, Jankowski MP, Albers KM, Koerber HR. Single-cell q-PCR derived expression profiles of identified sensory neurons. Mol Pain 2019;15:174480691988449.
. Arcourt A, Gorham L, Dhandapani R, Prato V, Taberner FJ, Wende H, Gangadharan V, Birchmeier C, Heppenstall PA, Lechner SG. Touch receptor-derived sensory information alleviates acute pain signaling and fine-tunes nociceptive reflex coordination. Neuron 2017;93:179–93.
. Bautista DM, Siemens J, Glazer JM, Tsuruda PR, Basbaum AI, Stucky CL, Jordt SE, Julius D. The menthol receptor TRPM8 is the principal detector of environmental cold. Nature 2007;448:204–8.
. Bautista DM, Wilson SR, Hoon MA. Why we scratch an itch: the molecules, cells and circuits of itch. Nat Neurosci 2014;17:175–82.
. Bouheraoua N, Fouquet S, Marcos-Almaraz MT, Karagogeos D, Laroche L, Chedotal A. Genetic analysis of the organization, development, and plasticity of corneal innervation in mice. J Neurosci 2019;39:1150–68.
. Chen KH, Boettiger AN, Moffitt JR, Wang S, Zhuang X. RNA imaging. Spatially resolved, highly multiplexed RNA profiling in single cells. Science 2015;348:aaa6090.
. Chiu IM, Barrett LB, Williams EK, Strochlic DE, Lee S, Weyer AD, Lou S, Bryman GS, Roberson DP, Ghasemlou N, Piccoli C, Ahat E, Wang V, Cobos EJ, Stucky CL, Ma Q, Liberles SD, Woolf CJ. Transcriptional profiling at whole population and single cell levels reveals somatosensory neuron molecular diversity. eLife 2014;3:e04660.
. Choi HMT, Schwarzkopf M, Fornace ME, Acharya A, Artavanis G, Stegmaier J, Cunha A, Pierce NA. Third-generation in situ hybridization chain reaction: multiplexed, quantitative, sensitive, versatile, robust. Development 2018;145:dev165753.
. DeNardo LA, Liu CD, Allen WE, Adams EL, Friedmann D, Fu L, Guenthner CJ, Tessier-Lavigne M, Luo L. Temporal evolution of cortical ensembles promoting remote memory retrieval. Nat Neurosci 2019;22:460–9.
. Dib-Hajj SD, Black JA, Waxman SG. NaV1.9: a sodium channel linked to human pain. Nat Rev Neurosci 2015;16:511–19.
. Dong X, Han S-k, Zylka MJ, Simon MI, Anderson DJ. A diverse family of GPCRs expressed in specific subsets of nociceptive sensory neurons. Cell 2001;106:619–32.
. Dussor G, Cao Y-Q. TRPM8 and migraine. Headache. J Head Face Pain 2016;56:1406–17.
. Gatto G, Smith KM, Ross SE, Goulding M. Neuronal diversity in the somatosensory system: bridging the gap between cell type and function. Curr Opin Neurobiol 2019;56:167–74.
. Ghitani N, Barik A, Szczot M, Thompson JH, Li C, Le Pichon CE, Krashes MJ, Chesler AT. Specialized mechanosensory nociceptors mediating rapid responses to hair pull. Neuron 2017;95:944–54.e944.
. Goadsby PJ, Reuter U, Hallstrom Y, Broessner G, Bonner JH, Zhang F, Sapra S, Picard H, Mikol DD, Lenz RA. A controlled trial of erenumab for episodic migraine. New Engl J Med 2017;377:2123–32.
. Han L, Limjunyawong N, Ru F, Li Z, Hall OJ, Steele H, Zhu Y, Wilson J, Mitzner W, Kollarik M, Undem BJ, Canning BJ, Dong X. Mrgprs on vagal sensory neurons contribute to bronchoconstriction and airway hyper-responsiveness. Nat Neurosci 2018;21:324–8.
. Han L, Ma C, Liu Q, Weng HJ, Cui Y, Tang Z, Kim Y, Nie H, Qu L, Patel KN, Li Z, McNeil B, He S, Guan Y, Xiao B, Lamotte RH, Dong X. A subpopulation of nociceptors specifically linked to itch. Nat Neurosci 2013;16:174–82.
. Huang CC, Yang W, Guo C, Jiang H, Li F, Xiao M, Davidson S, Yu G, Duan B, Huang T, Huang AJW, Liu Q. Anatomical and functional dichotomy of ocular itch and pain. Nat Med 2018;24:1268–76.
. Julius D. TRP channels and pain. Annu Rev Cel Dev Biol 2013;29:355–84.
. Khoutorsky A, Price TJ. Translational control mechanisms in persistent pain. Trends Neurosci 2018;41:100–14.
. Kunkler PE, Ballard CJ, Oxford GS, Hurley JH. TRPA1 receptors mediate environmental irritant-induced meningeal vasodilatation. PAIN 2011;152:38–44.
. Kunkler PE, Ballard CJ, Pellman JJ, Zhang L, Oxford GS, Hurley JH. Intraganglionic signaling as a novel nasal-meningeal pathway for TRPA1-dependent trigeminovascular activation by inhaled environmental irritants. PLoS One 2014;9:e103086.
. LaMotte RH, Dong X, Ringkamp M. Sensory neurons and circuits mediating itch. Nat Rev Neurosci 2014;15:19–31.
. Li CL, Li KC, Wu D, Chen Y, Luo H, Zhao JR, Wang SS, Sun MM, Lu YJ, Zhong YQ, Hu XY, Hou R, Zhou BB, Bao L, Xiao HS, Zhang X. Somatosensory neuron types identified by high-coverage single-cell RNA-sequencing and functional heterogeneity. Cell Res 2016;26:83–102.
. Lignell A, Kerosuo L, Streichan SJ, Cai L, Bronner ME. Identification of a neural crest stem cell niche by Spatial Genomic Analysis. Nat Commun 2017;8:1830.
. Macpherson LJ, Dubin AE, Evans MJ, Marr F, Schultz PG, Cravatt BF, Patapoutian A. Noxious compounds activate TRPA1 ion channels through covalent modification of cysteines. Nature 2007;445:541–5.
. Mishra SK, Hoon MA. The cells and circuitry for itch responses in mice. Science 2013;340:968–71.
. Mitchell K, Lebovitz EE, Keller JM, Mannes AJ, Nemenov MI, Iadarola MJ. Nociception and inflammatory hyperalgesia evaluated in rodents using infrared laser stimulation after Trpv1 gene knockout or resiniferatoxin lesion. PAIN 2014;155:733–45.
. Moeyaert B, Holt G, Madangopal R, Perez-Alvarez A, Fearey BC, Trojanowski NF, Ledderose J, Zolnik TA, Das A, Patel D, Brown TA, Sachdev RNS, Eickholt BJ, Larkum ME, Turrigiano GG, Dana H, Gee CE, Oertner TG, Hope BT, Schreiter ER. Improved methods for marking active neuron populations. Nat Commun 2018;9:4440.
. Moffitt JR, Bambah-Mukku D, Eichhorn SW, Vaughn E, Shekhar K, Perez JD, Rubinstein ND, Hao J, Regev A, Dulac C, Zhuang X. Molecular, spatial, and functional single-cell profiling of the hypothalamic preoptic region. Science 2018;362:eaau5324.
. Nguyen MQ, Le Pichon CE, Ryba N. Stereotyped transcriptomic transformation of somatosensory neurons in response to injury. eLife 2019;8:e49679.
. Nguyen MQ, Wu Y, Bonilla LS, von Buchholtz LJ, Ryba NJP. Diversity amongst trigeminal
neurons revealed by high throughput single cell sequencing. PLoS One 2017;12:e0185543.
. Northcutt AJ, Kick DR, Otopalik AG, Goetz BM, Harris RM, Santin JM, Hofmann HA, Marder E, Schulz DJ. Molecular profiling of single neurons of known identity in two ganglia from the crab Cancer borealis. Proc Natl Acad Sci U S A 2019;116:26980–90.
. Osteen JD, Herzig V, Gilchrist J, Emrick JJ, Zhang C, Wang X, Castro J, Garcia-Caraballo S, Grundy L, Rychkov GY, Weyer AD, Dekan Z, Undheim EA, Alewood P, Stucky CL, Brierley SM, Basbaum AI, Bosmans F, King GF, Julius D. Selective spider toxins reveal a role for the Nav1.1 channel in mechanical pain. Nature 2016;534:494–9.
. Pal NR, Pal SK. A review on image segmentation techniques. Pattern Recognition 1993;26:1277–94.
. Parra A, Madrid R, Echevarria D, del Olmo S, Morenilla-Palao C, Acosta MC, Gallar J, Dhaka A, Viana F, Belmonte C. Ocular surface wetness is regulated by TRPM8-dependent cold thermoreceptors of the cornea. Nat Med 2010;16:1396–9.
. Pogorzala LA, Mishra SK, Hoon MA. The cellular code for mammalian thermosensation. J Neurosci 2013;33:5533–41.
. Ronneberger O, Fischer P, Brox T. U-net: convolutional networks for biomedical image segmentation. In: Navab N, Hornegger J, Wells W, Frangi A (eds). Medical Image Computing and Computer-Assisted Intervention – MICCAI 2015. MICCAI 2015. Lecture Notes in Computer Science, vol 9351.
. Sathyamurthy A, Johnson KR, Matson KJE, Dobrott CI, Li L, Ryba AR, Bergman TB, Kelly MC, Kelley MW, Levine AJ. Massively parallel single nucleus transcriptional profiling defines spinal cord neurons and their activity during behavior. Cell Rep 2018;22:2216–25.
. Satija R, Farrell JA, Gennert D, Schier AF, Regev A. Spatial reconstruction of single-cell gene expression
data. Nat Biotechnol 2015;33:495–502.
. Shah S, Lubeck E, Zhou W, Cai L. seqFISH accurately detects transcripts in single cells and reveals robust spatial organization in the Hippocampus. Neuron 2017;94:752–8.e1.
. Sharma N, Flaherty K, Lezgiyeva K, Wagner DE, Klein AM, Ginty DD. The emergence of transcriptional identity in somatosensory neurons. Nature 2020;577:392–8.
. Silberstein SD, Dodick DW, Bigal ME, Yeung PP, Goadsby PJ, Blankenbiller T, Grozinski-Wolff M, Yang R, Ma Y, Aycardi E. Fremanezumab for the preventive treatment of chronic migraine. New Engl J Med 2017;377:2113–22.
. Solinski HJ, Kriegbaum MC, Tseng PY, Earnest TW, Gu X, Barik A, Chesler AT, Hoon MA. Nppb neurons are sensors of mast cell-induced itch. Cell Rep 2019;26:3561–73.e4.
. Sullivan MN, Gonzales AL, Pires PW, Bruhl A, Leo MD, Li W, Oulidi A, Boop FA, Feng Y, Jaggar JH, Welsh DG, Earley S. Localized TRPA1 channel Ca2+signals stimulated by reactive oxygen species promote cerebral artery dilation. Sci Signaling 2015;8:ra2.
. Tasic B, Yao Z, Graybuck LT, Smith KA, Nguyen TN, Bertagnolli D, Goldy J, Garren E, Economo MN, Viswanathan S, Penn O, Bakken T, Menon V, Miller J, Fong O, Hirokawa KE, Lathia K, Rimorin C, Tieu M, Larsen R, Casper T, Barkan E, Kroll M, Parry S, Shapovalova NV, Hirschstein D, Pendergraft J, Sullivan HA, Kim TK, Szafer A, Dee N, Groblewski P, Wickersham I, Cetin A, Harris JA, Levi BP, Sunkin SM, Madisen L, Daigle TL, Looger L, Bernard A, Phillips J, Lein E, Hawrylycz M, Svoboda K, Jones AR, Koch C, Zeng H. Shared and distinct transcriptomic cell types across neocortical areas. Nature 2018;563:72–8.
. Umans BD, Liberles SD. Neural sensing of organ volume. Trends Neurosci 2018;41:911–24.
. Usoskin D, Furlan A, Islam S, Abdo H, Lonnerberg P, Lou D, Hjerling-Leffler J, Haeggstrom J, Kharchenko O, Kharchenko PV, Linnarsson S, Ernfors P. Unbiased classification of sensory neuron types by large-scale single-cell RNA sequencing. Nat Neurosci 2015;18:145–53.
. Wang X, Allen WE, Wright MA, Sylwestrak EL, Samusik N, Vesuna S, Evans K, Liu C, Ramakrishnan C, Liu J, Nolan GP, Bava FA, Deisseroth K. Three-dimensional intact-tissue sequencing of single-cell transcriptional states. Science 2018;361:eaat5691.
. Winnubst J, Bas E, Ferreira TA, Wu Z, Economo MN, Edson P, Arthur BJ, Bruns C, Rokicki K, Schauder D, Olbris DJ, Murphy SD, Ackerman DG, Arshadi C, Baldwin P, Blake R, Elsayed A, Hasan M, Ramirez D, Dos Santos B, Weldon M, Zafar A, Dudman JT, Gerfen CR, Hantman AW, Korff W, Sternson SM, Spruston N, Svoboda K, Chandrashekar J. Reconstruction of 1,000 projection neurons reveals new cell types and organization of long-range connectivity in the mouse brain. Cell 2019;179:268–81.e213.
. Zeisel A, Hochgerner H, Lonnerberg P, Johnsson A, Memic F, van der Zwan J, Haring M, Braun E, Borm LE, La Manno G, Codeluppi S, Furlan A, Lee K, Skene N, Harris KD, Hjerling-Leffler J, Arenas E, Ernfors P, Marklund U, Linnarsson S. Molecular architecture of the mouse nervous system. Cell 2018;174:999–1014.e1022.
. Zheng Y, Liu P, Bai L, Trimmer JS, Bean BP, Ginty DD. Deep sequencing of somatosensory neurons reveals molecular determinants of intrinsic physiological properties. Neuron 2019;103:598–616.e597.
. Zimmerman A, Bai L, Ginty DD. The gentle touch receptors of mammalian skin. Science 2014;346:950–4.
. Zimmermann K, Leffler A, Babes A, Cendan CM, Carr RW, Kobayashi J, Nau C, Wood JN, Reeh PW. Sensory neuron sodium channel Nav1.8 is essential for pain at low temperatures. Nature 2007;447:855–8.