Abstract
Viscera receive innervation from sensory ganglia located adjacent to multiple levels of the brainstem and spinal cord. Here we examined whether molecular profiling could be used to identify functional clusters of colon afferents from thoracolumbar (TL), lumbosacral (LS), and nodose ganglia (NG) in male and female mice. Profiling of TL and LS bladder afferents was also performed. Visceral afferents were back-labeled using retrograde tracers injected into proximal and distal regions of colon or bladder, followed by single-cell qRT-PCR and analysis via an automated hierarchical clustering method. Genes were chosen for assay (32 for bladder; 48 for colon) based on their established role in stimulus detection, regulation of sensitivity/function, or neuroimmune interaction. A total of 132 colon afferents (from NG, TL, and LS ganglia) and 128 bladder afferents (from TL and LS ganglia) were analyzed. Retrograde labeling from the colon showed that NG and TL afferents innervate proximal and distal regions of the colon, whereas 98% of LS afferents only project to distal regions. There were clusters of colon and bladder afferents, defined by mRNA profiling, that localized to either TL or LS ganglia. Mixed TL/LS clustering also was found. In addition, transcriptionally, NG colon afferents were almost completely segregated from colon TL and LS neurons. Furthermore, colon and bladder afferents expressed genes at similar levels, although different gene combinations defined the clusters. These results indicate that genes implicated in both homeostatic regulation and conscious sensations are found at all anatomic levels, suggesting that afferents from different portions of the neuraxis have overlapping functions.
SIGNIFICANCE STATEMENT Visceral organs are innervated by sensory neurons whose cell bodies are located in multiple ganglia associated with the brainstem and spinal cord. For the colon, this overlapping innervation is proposed to facilitate visceral sensation and homeostasis, where sensation and pain are mediated by spinal afferents and fear and anxiety (the affective aspects of visceral pain) are the domain of nodose afferents. The transcriptomic analysis performed here reveals that genes implicated in both homeostatic regulation and pain are found in afferents across all ganglia types, suggesting that conscious sensation and homeostatic regulation are the result of convergence, and not segregation, of sensory input.
Introduction
Visceral organs receive sensory innervation from primary sensory neurons arising from multiple levels of the neuraxis. Depending on the organ, innervation can originate from vagal afferents [including afferents in the jugular ganglia and nodose ganglia (NG)], thoracolumbar (TL) spinal afferents [arising from lower thoracic and upper lumbar dorsal root ganglia (DRG)], or lumbosacral (LS) spinal afferents (arising from lower lumbar and upper sacral DRG). Our hypotheses to account for why viscera requires monitoring by multiple populations of afferents include the following: (1) that different levels of innervation are involved in different qualitative aspects of stimulus detection, including pain; (2) that because the first synapse is on CNS neurons associated with sympathetic or parasympathetic circuits, different afferent populations provide the basis for integration of autonomic/homeostatic functions; or (3) that the different levels play complementary roles in immune modulation.
With regard to shaping the nature of the visceral sensory experience, it had been hypothesized that spinal afferents, but not vagal afferents, relay nociceptive sensations (e.g., stabbing, burning, cramping), whereas vagal afferents underlie affective aspects of visceral pain, including depression, anxiety, nausea, and fear (Altschuler et al., 1993; Berthoud and Neuhuber, 2000; Grundy, 2002; Sengupta, 2009). However, recent evidence has shown that vagal afferents may have a significant role in nociception and pain (Lee et al., 2003; Kang et al., 2004; Yu et al., 2005; Bielefeldt et al., 2006; Canning et al., 2006; Nassenstein et al., 2008; Taylor-Clark et al., 2009). Even spinal afferents from different levels may have specific roles; colonic afferents arising from the TL level have been reported to be important for inflammatory pain but not for acute pain, whereas LS colonic afferents are proposed to be involved in both acute and inflammatory pain (Traub et al., 1994; Traub, 2000; Wang et al., 2005).
In terms of homeostasis, spinal visceral afferents that innervate colon and bladder have cell bodies in both TL and LS DRG projecting to organs they innervate along sympathetic (TL) or parasympathetic (LS) splanchnic nerves and terminate broadly in the spinal cord, including in the region of sympathetic (TL) or parasympathetic (LS) preganglionic neurons (Harrington et al., 2019; Smith-Edwards et al., 2019). In contrast, vagal afferents that innervate the colon (and possibly the bladder; Herrity et al., 2014) have cell bodies in the nodose ganglia, run in the vagus nerve with preganglionic parasympathetic fibers, and terminate in the brainstem nucleus tractus solitarius, which in turn projects to nuclei in the brainstem, thalamus, and cortex (Norgren, 1978; van der Kooy et al., 1984). Thus, redundant sensory input to the CNS could be informing different components of the autonomic nervous system.
Afferent/immune interactions are also critical for gut function and homeostasis (i.e., the vagus nerve has a crucial role in the anti-inflammatory reflex; Tracey, 2002), with vagal afferent response to cytokine release as the first step in this reflex. Primary afferents (vagal and spinal) release multiple peptides including calcitonin gene-related peptide (CGRP) and substance P, which bind receptors on immune cells, allowing direct afferent and immune communication (Brain and Williams, 1985; Caceres et al., 2009; Altmayr et al., 2010; Riol-Blanco et al., 2014; Cohen et al., 2019). Ablation or silencing of primary afferents (vagal and spinal) also can modulate the immune response across multiple organs, including the lungs, small intestine, and colon (Engel et al., 2011; Baral et al., 2018; Lai et al., 2020).
In this study we examined how genes required for sensation, homeostasis and neuroimmune regulation are divided among different visceral afferents. Single-cell qRT-PCR analysis combined with automated hierarchical clustering (AHC) of afferents innervating the bladder and colon were conducted. Results show that location matters; that genes expressed by an afferent and where that afferent terminates (especially for the colon) is determined by the ganglionic location. Furthermore, the clustering of unique genes to afferents from different ganglia suggests that all afferents contribute to homeostatic and immune regulation, as well as signal all aspects of visceral sensation.
Materials and Methods
Animals.
Experiments were conducted on male and female adult (6–12 week of age) C57BL/6 mice from The Jackson Laboratory (catalog #000664). Animals were group housed with a 12 h light/dark cycle and ad libitum access to food and water. All procedures were approved by the Institutional Animal Care and Use Committee at the University of Pittsburgh and were conducted in accordance with The Association for Assessment and Accreditation of Laboratory Animal Care International-approved practices.
Back-labeling of neurons.
Mice were anesthetized with isoflurane (2%), and a laparotomy was performed to access the pelvic viscera. Fluorescent retrograde dyes [cholera toxin subunit β (CTB) or wheat germ agglutinin (WGA); Thermo Fisher Scientific] were injected into colon or bladder as previously described (Christianson et al., 2007). Briefly, using an insulin syringe (31 gauge needle), 5–10 μl of Alexa Fluor 555-conjugated CTB was injected in 2–3 μl aliquots beneath the serosal layer of the distal colon at the level of the base of the bladder. Alexa Fluor 488-conjugated CTB was injected in the proximal colon just distal to the cecum. To ensure that dye did not spread between areas, colons were removed and visualized under a fluorescence microscope. In all cases, there was at least a 20 mm area of middle colon that had no visible fluorescence. In a separate cohort of mice, 5-10 μl of Alexa Fluor 488-conjugated CTB or WGA was injected in 2–3 μl aliquots beneath the serosal layer in the urinary bladder body. WGA was used to examine whether CTB was labeling a specific subset of afferents (n = 3, three females). We found that the number of cells expressing a given gene was virtually the same (except for some of the rarely detected genes like Mrgprd and Trpm8) regardless of which tracer was used (Table 1). Therefore, CTB was used in subsequent experiments. In a subset of these animals, 555-conjugated CTB was injected into both proximal and distal colon, combined with 488-CTB injection into the bladder to label convergent bladder/colon afferents. All incisions were then sutured, and animals were allowed to recover before returning to their home cage.
Dissociation of ganglia and single neuron isolation.
Three to 5 d after back-labeling, mice were anesthetized with isoflurane and transcardially perfused with cold Ca2+/Mg2+-free HBSS. The NG, TL DRG (T10–L1), and LS DRG (L5–S1) were dissected into cold HBSS and enzymatically treated with cysteine, papain, collagenase type II, and dispase type II to facilitate isolation by mechanical trituration. Cells were plated on poly-d-lysine/laminin-coated coverslips in 35mmx10mm Petri dishes. Coverslips were flooded with DMEM F12 containing 10% fetal bovine serum and 1% penicillin/streptomycin 2 h later. Cells were picked up from each culture using glass capillaries (World Precision Instruments) held by a three-axis micromanipulator. Typically, there were 10–20 back-labeled cells on each coverslip. To avoid selection bias based on cell size, only isolated single neurons (those not adjacent to another labeled cell) were used. Six to 12 coverslips were plated for each mouse, and cells were sampled from multiple coverslips. To monitor for size selection bias, the relative size of each cell was recorded at the time of collection. The great majority of back-labeled cells were small (diameter, 10–20 µm) or medium (diameter, 20–30 µm). Occasionally, large cells were seen and collected in a manner proportional to their relative abundance. Collected cells were transferred to tubes containing 3 μl of lysis buffer (MessageBOOSTER Kit, Epicentre), and stored at –80°.
Single-cell amplification and qPCR.
The methods for single-cell qRT-PCR and clustering, and validation of these techniques have been previously described (Adelman et al., 2019). Briefly, the RNA collected from each cell was reverse transcribed and amplified using T7 linear amplification (MessageBOOSTER Kit for Cell Lysate, Epicentre), cleaned with RNA Cleaner & Concentrator-5 columns (Zymo Research), and assessed using qPCR with optimized primers and SsoAdvanced SYBR Green Master Mix (BIO-RAD). Cycle time (Ct) values were determined via regression. Quantification threshold for PCR was defined as the point at which there was a 95% replication rate (35 Ct; Reiter et al., 2011). GAPDH threshold was thus defined as 25 Ct to ensure detection of transcripts at least a thousand-fold less common than GAPDH.
Primer design and validation.
Unique forward and reverse primer sequences were chosen for each gene within 500 bases of the 3′ end. Stock solutions of cDNA were generated by extracting RNA from the whole DRG and both 10 and 160 pg aliquots of the RNA were amplified using the same procedure described above for single cells. Serial dilutions of these aliquots were used to calculate primer efficiencies over the range of RNA concentrations observed in single cells. Expression level was determined relative to GAPDH and corrected for these primer efficiencies (Pfaffl, 2001).
Automated hierarchical clustering.
All back-labeled neurons were clustered using the unweighted pair group method with averaging (UPGMA) using expression information obtained from single-cell PCR in MATLAB (MathWorks). The preprocessing for this data analysis consists of taking the ΔCt values and replacing the samples that failed to generate a value for a given gene with the limit of detection for that gene.
Tissue preparation and estimation of the number of back-labeled neurons.
Three to 5 d after back-labeling, mice were anesthetized with isoflurane and transcardially perfused with cold isotonic saline. The NG, and TL and LS DRG were removed and fixed with 4% paraformaldehyde solution for 30 min. For a negative control, the L3 DRG was also removed. After postfixation, the ganglia were moved to 25% sucrose for overnight cryoprotection. After cryoprotection, ganglia were embedded in optimal cutting temperature compound (Thermo Fisher Scientific) and cryosectioned at 14 µm, allowed to dry, and immediately stored at −80°C. Each ganglion was placed on six slides with serial sections alternating between slides. Ganglia sections were analyzed using fluorescent microscopy at 20×. Only two slides from the six were analyzed with 42 µm between slides; visual confirmation was used to avoid analyzing the same cells more than once. At least six sections from each ganglion were analyzed per animal. For each animal, the number of 555-conjugated CTB-positive, 488-conjugated CTB-positive cells, and double-labeled cells were counted. The percentage of CTB-positive cells at each level of the neuraxis was found for each mouse and averaged within groups.
Calcium imaging protocols.
Colon neurons were back-labeled and prepared for culture as described above, except that the proximal and distal colon were both labeled with CTB-555. Cells were incubated at 37°C and imaged within 36 h, as previously described (Malin et al., 2006). Before imaging, cells were incubated for 30 min at 37°C with 1 μl Invitrogen fura-2 AM (2 μm; Thermo Fisher Scientific) and 2 μl of 20% Pluronic F-127 dissolved in DMSO (AnaSpec) and diluted in HBSS containing 10 mg/ml BSA (Sigma-Aldrich). Coverslips were mounted on an inverted Olympus IX71 microscope stage with HBSS buffer flowing at 5 ml/min, controlled by a Perfusion Fast-Step System (AutoMate Scientific). Perfusate temperature was maintained at 30°C using a heated stage and an inline heating system (Warner Instruments). Chemicals were delivered with a rapid-switching local perfusion system. Firmly attached, CTB-positive neurons were identified using a 555 nm filter and chosen as regions of interest using Simple HCI software (Compix Imaging Systems). Unlabeled adjacent cells were also identified and imaged. All fields were first tested with a brief application (4 s) of 50 mm K+ (high K+) to ensure that cells were responsive. Following a 5 min recovery period, agonists were applied in a randomized order with at least a 5 min recovery period between agonists. Stock solutions of capsaicin (Sigma-Aldrich) and mustard oil (Sigma-Aldrich) were made in 1-methyl-2-pyrrolidinone, and α,β-methylene ATP (Sigma-Aldrich) was made in water. The stock solution of mouse interferon-α A (PBL Assay Science); recombinant mouse interleukin-4 (PeproTech) and recombinant mouse interferon-γ and recombinant human interkeukin-8 (R&D Systems) was prepared in PBS containing 0.1% BSA, as previously described (Wang et al., 2017). Final concentrations applied to cells were 1 μm capsaicin, 30 μm α,β-methylene ATP and 50 μm mustard oil for 4 s; and 500 ng/ml interleukin-4, interleukin-8, and interferon-γ, and 1000 IU/ml interferon-α for 90 s. Absorbance data at 340 and 380 nm were collected at 1 frame/s. Responses were measured as the ratio of 340/380 nm excitation and 510 nm emission over baseline (ΔF/F0; DG4, Sutter Instrument). Peak responses were included in the analysis if the response was 5 SDs above baseline. The prevalence of responsive colon afferents was determined as a percentage of total healthy (high K+-responsive) CTB-positive cells. Any cell with a significantly diminished fura-2 signal over the duration of the experiment or a signal that did not recover to baseline before the second agonist application was not included in the analysis.
Experimental design and statistical analysis.
For back-labeling experiments, the mean values of the number of back-labeled cells were compared between level (NG, TL, LS) and target (proximal, distal, dual) using a two-way ANOVA with Tukey's post hoc analysis. A two-tailed paired t test was used to test for significance of the difference in the extent of spread (length) of CTB injections in the proximal and distal colon. For calcium imaging experiments, mean peak values of the calcium transients were compared between level (NG, TL, LS) and target (proximal, distal, dual) using a mixed-effects model with Tukey's post hoc analysis. Statistical tests were performed in Excel (Microsoft) and GraphPad Prism (GraphPad Software). Data are expressed as the mean ± SEM, where n represents mice used, unless indicated otherwise. The number of animals and statistical values are reported in the Results section; significance is defined as p < 0.05.
Results
Distribution of back-labeled neurons from injection sites in distal and proximal colon
In the mouse, the spinal innervation of the distal colon and bladder (including dually projecting colon and bladder neurons) have been described (Robinson et al., 2004; Christianson et al., 2007), whereas innervation of the proximal colon has not been fully characterized. To fill this knowledge gap, the spinal and nodose sensory innervation of the proximal colon was investigated and compared with innervation of the distal colon.
Injection of CTB into the proximal and distal colon of individual mice labeled, on average, a total of 659 ± 75 CTB-positive neurons within the nodose and spinal ganglia (n = 6, 4 males, 2 females; Table 2, Fig. 1A). There were nearly equal percentages of neurons arising from nodose (32.29 ± 3.73%, p > 0.05), TL (30.29 ± 4.02%), and LS (37.43 ± 5.77%) ganglia. Injections into the distal colon produced more CTB-labeled cells (451 ± 84) compared with the proximal colon (181 ± 54; F(3,15) = 22.19, p = 0.0256; Table 2), resulting in 65.42 ± 8.12% of CTB-positive afferents originating from distal colon, 27.87 ± 7.55% originating from proximal, colon and 6.70 ± 1.02% being dual labeled (i.e., innervating both proximal and distal colon). This difference in the number of cells innervating the proximal and distal colons occurs because the majority of LS afferents project to the distal colon; of the 258 ± 56 back-labeled cells in LS ganglia, 252 ± 59 were labeled by injections into the distal colon (this included 6 ± 1 dual-labeled afferents; Table 2). Thus, LS colon afferents almost exclusively innervate the distal colon. For the proximal colon, the largest afferent contribution came from nodose and TL ganglia [58.29 ± 7.20% (NG), 38.22 ± 6.00% (TL ganglia), and 3.49 ± 1.63% (LS ganglia); F(4,20) = 25.21, p < 0.0001 for NG vs LS ganglia; p = 0.0002 for TL vs LS ganglia; Fig. 1D]. The size of the proximal and distal injections was not a factor as the distances over which fluorescence could be detected in the proximal colon (13.0 ± 1.9 mm) and distal colon (14.0 ± 1.1 mm; t(7) = 0.5739, p = 0.58) were equivalent (n = 8, 5 males, 3 females; Fig. 1E). In addition, there was a 20–30 mm gap between the leading edge of each injection. That TL afferents projected equally to proximal and distal colon, while LS afferents were more restricted to the distal colon is consistent with the idea that neural components that regulate colon function change along the length of the colon (Li et al., 2019).
Automated hierarchical clustering based on mRNA expression and afferent location reveal distinct neuronal clusters
To classify afferent subtypes we used single-cell qRT-PCR to analyze the transcripts of 28 genes previously used in a study for classification of identified cutaneous afferents (Adelman et al., 2019). These genes were shown to identify clusters that largely replicated the clusters produced by bulk, single-cell RNA sequencing (RNAseq) of DRG neurons unidentified with respect to innervation target or functional phenotype (Usoskin et al., 2015; Zeisel et al., 2018). Importantly Adelman et al. (2019) combined AHC of these 28 genes with functional identification of neuron subtype with the use of an ex vivo skin-nerve physiological preparation. This approach allowed correlation of the molecular phenotype of individual cells with physiological response properties of intact afferents (e.g., conduction velocity, response to mechanical and thermal stimuli). Using this approach, these 28 genes produced clusters that corresponded to functional afferent subtypes [e.g., CMH (C-mechanoheat), CH (C-heat), and HTMR (high-threshold mechanoreceptors)]. For analysis of colon and bladder afferents, four additional genes encoding opioid receptors and the immune checkpoint inhibitor programmed cell death 1 ligand 1 [Cd274 (also called programmed death-ligand 1 or PDL-1] were added. For colon afferents only, 16 genes associated with immune function were also examined. Based on previous RNAseq analysis, these immune genes are highly expressed in both the NG (Wang et al., 2017) and lumbar DRG (Usoskin et al., 2015).
The protocols used here were previously shown to produce a linear correlation between starting mRNA transcript level and the amount of mRNA detected following reverse transcription and amplification (Adelman et al., 2019). Thus, for each gene we can determine the relative level of expression for each cell analyzed (but cannot compare the level of expression between mRNA transcript species because of differential primer efficiency). Expression level information was converted into the heat maps shown in Figures 2 and 3. In addition, descriptors of expression level in each cluster are reported as “undetected” (every cell in a group had no detectable level of transcript), “low” (when the average expression level for cells in a cluster was in the lowest 30th percentile), “moderate” (between 31st and 70th percentiles) or “high” (between 71st and 100th percentiles).
Clusters were named by combining the approaches of Zeisel et al. (2018) and Hockley et al. (2019). Abbreviations used for clusters are as follows: PEP for peptidergic [clusters containing cells with expression of mRNA coding for substance P (Tac1) and calcitonin gene-related peptide (Calca)]; NP for nonpeptidergic (clusters of cells not expressing neuropeptide related mRNAs); NF for neurofilament (clusters of cells with expression of the high-molecular-weight neurofilament gene (Nefh), a marker of myelinated afferents; Lawson and Waddell, 1991; Adelman et al., 2019); tl for clusters containing cells originating predominately (>85%) from thoracolumbar ganglia; ls for clusters of cells from lumbosacral ganglia; ng for clusters of cells from nodose ganglia; and m for mixed clusters containing cells from TL and LS ganglia. Clusters are named based on whether they are ng, ls, tl, or m plus whether they are PEP, NP, NF, or some combination (e.g., NF-PEP). Because it is not possible to interrogate the unbiased algorithms used in this study, the contribution of each gene to the decision-making process is unknown. Therefore, genes that are highlighted in the discussion of each cluster were chosen based on our inspection of the heatmaps and the identification of ones most highly expressed in a cluster and ones that are either absent or expressed at low levels.
Colon afferents
Using the previously validated AHC methods of Adelman et al. (2019), a total of 132 colon afferents from NG (n = 43), TL ganglia (n = 48), and LS ganglia (n = 41) separated into 13 distinct clusters based on mRNA level (n = 7 animals, 4 males, 3 females; Figs. 2, 4A). Sixteen cells did not fit into any of these clusters and are not represented on the heatmap in Figure 2.
As the AHC is unbiased, it is interesting that the first branch point in the cell dendrogram separates nodose afferents from spinal afferents (Figs. 2, 4A). Major distinguishing characteristics of nodose afferents were the high levels of purinergic receptors (P2rx2, P2rx3) and Cd274, and the low expression of calcitonin gene-related polypeptide a (Calca) and preprotachykinin-1 (Tac1). Nodose afferents were further divided into four clusters based on genes involved in growth factor receptor signaling, chemokine receptors, and GPCR expression.
For AHC of spinal afferents, the first branch point separated cells based on levels of the heavy neurofilament (Nefh), a marker of myelination. The Nefh branch point gave rise to five Nefh-low and four Nefh-high final clusters. In the Nefh-low group, two clusters came almost exclusively from TL ganglia and one came from LS ganglia, whereas two other clusters contained a mixture of TL and LS neurons. For the Nefh-high clusters, two contained only LS afferents and the other two contained a mix of TL and LS neurons. Thus, even at the initial level of analysis, the AHC revealed multiple mRNA-based clusters that were unique to afferents arising from anatomically distinct neuron populations.
Nodose afferent clusters
Four nodose clusters with nearly equal numbers of cells (8–11 in each cluster) were identified (Figs. 2, 4A). Each cluster had both proximal colon and distal colon projecting cells. The ngNPa cluster has high levels of interleukin-8 receptor type 2 (Cxcr2), prostaglandin I2 receptor (Ptgir), and interferon-γ receptor 2 (Ifngr2). The ngNPb cluster also had high Cxrc2, interleukin-1 receptor type 1 (Il1r1), and cholinergic receptor nicotinic α 3 subunit (Chrna3). The ngNF cluster was a distinct nodose cluster with high expression of Nefh, glial cell line-derived neurotrophic factor family receptor α 2 (Gfra2; one of three coreceptors for Ret), sodium voltage-gated channel α subunit 10 (Scn10a), purinergic receptor P2ry1 (a Gq-coupled GPCR detected only with uniform, moderate expression in this one cluster), and opioid receptor κ 1 (Oprk1). The nodose ngNF group was also unique in the uniform and low level of transient receptor potential cation channel family v member 1 (Trpv1) and Chrna3. The final nodose cluster, ngNPc, had high expression of Scn10a and glial cell line-derived neurotrophic factor family receptor α 1 (Gfra1, along with moderate expression of its coreceptor Ret) and was the only nodose cluster with a moderate level of transient receptor potential cation channel family a member 1 (Trpa1).
Spinal afferents
Thoracolumbar clusters
Of the two TL clusters (tlPEPa and tlPEPb), tlPEPa, overall the largest cluster, had high levels of Calcα and Ret. This cluster also had moderate expression of glial cell line-derived neurotrophic factor family receptor α 3 (Gfra3; which was absent in nodose clusters), Tac1, tropomyosin receptor kinase a (Trka), Trpa1, Trpv1, and interleukin-4 receptor α (Il4ra; Figs. 2, 4A). The tlPEPa cluster neurons project primarily to the proximal colon. As anatomic origin (proximal vs distal) was not a variable in the AHC algorithm, this clustering indicates that transcripts of proximal colon afferents are unique from those of the distal colon. The tlPEPb cluster, the second smallest of the spinal groups, had moderate Il4ra, low Trpa1, and high Trka, relative to the tlPEPa cluster. This cluster also had high interleukin-1 receptor type 1 (Il1r1) and moderate transient receptor potential cation channel family m member 8 (Trpm8), purinergic receptor P2ry2, and Oprk1.
Lumbosacral clusters
Within the four identified LS clusters, 20 of 23 of the back-labeled afferents were from the distal colon (two were dually labeled), indicating that the distal colon receives functionally distinct afferent input relative to the proximal colon. All but the lsPEP cluster had high or moderate Nefh, piezo type mechanosensitive ion channel component 2 (Piezo2), Ifngr2, and Gfra1. The sodium voltage-gated channel α subunit 9 (Scn9a) was expressed in all LS neurons (including the lsPEP cluster) as was the acid sensing ion channel 2 (Asic2). The lsNF-PEP cluster was distinguished by high levels of the immune genes interferon-γ receptor 1 (Ifngr1), Ifngr2, Cxcr2, Il1rl, Ptgir, and interferon-α receptor 2 (Ifnar2). This cluster also had high Chrna3 and moderate Th, P2rx2, Tac1, and Calca. The lsNFa and lsNFb clusters were the most similar of the LS clusters, distinguished by relatively high levels of Ifngr1 in the lsNFa cluster and Ifngr2 in the lsNFb cluster. The lsNFb cluster also had moderate tropomyosin receptor kinase b (Trkb), which was absent from other spinal clusters except the lsNF-PEP cluster. As mentioned above, the lsPEP cluster had transcripts not detected in other LS clusters and was unique in having the high levels of Trpv1 and tyrosine hydroxylase (Th).
Mixed clusters
There were three mixed clusters of spinal afferents that were defined by neurons within both TL and LS ganglia (Figs. 2, 4A). The mPEPa and mPEPb clusters had high or moderate Calcα, Gfra3, Ret, and Trpv1 transcripts. The mPEPa cluster also had moderate levels of the immune genes Il4ra, Ifngr2, and fractalkine receptor (Cx3cr1). The mPEPb cluster was the only cluster with high transient receptor potential cation channel family m member 3 (Trpm3) and Oprk1, a gene with minimal expression in other spinal clusters (although present at high levels in one nodose cluster).
The other mixed cluster (mNF) was similar to the LS clusters with respect to the relative levels of Nefh, Piezo2, Scn9a, Gfra1, and Asic2. This cluster also contained the immune-related genes Ifnar1, Ifngr2, Cd274, Cxcr2, Ifngr1, Ptgir, and Cx3cr1 as well as P2rx3, acid-sensing ion channel 1 (Asic1), acid-sensing ion channel 3 (Asic3), and tropomyosin receptor kinase c (Trkc). This was the only cluster (spinal or NG) that had no detectable level of Chrna3.
Comparison of gene clustering of vagal and spinal colon afferents
A significant finding of the current study is the nearly exclusive manner in which colon NG neurons segregate from colon spinal afferents. Only one colon DRG neuron appeared in the four nodose clusters, and only one nodose neuron sorted with neurons in the spinal clusters. The major contributor to this segregated pattern was the uniformly high levels of P2rx2 and P2rx3 in NG neurons. Similar high levels of the prostacyclin receptor (Ptgir) and Cd274 [AKA PDL-1 (programmed cell death ligand 1)] were in the NG. Ptgir and Cd274 have important roles in immune responses, consistent with the role of nodose afferents in neuroimmune regulation (Sharkey and Mawe, 2002; Tracey, 2002; Baral et al., 2018). The individual NG clusters were also unique. No NG cluster had characteristics similar to the peptidergic clusters at TL or LS levels. There was a conspicuous absence or low level of Trka, Calca, Gfra3, and Tac1 in NG neurons, consistent with previous findings (Wang et al., 2017; Trancikova et al., 2018; Mazzone et al., 2020).
Bladder afferent clusters
For the bladder, 32 genes were analyzed from 128 back-labeled afferents from LS (80 afferents) and TL (48 afferents) levels (n = 11 animals, 3 males, 8 females; Figs. 3, 4B). Neurons sorted into seven distinct clusters with nine cells not fitting into any cluster (not shown in the heat map). Because some afferents can project to both bladder and colon (Keast and De Groat, 1992; Christianson et al., 2007), in some experiments CTB conjugated to different fluorophores was injected into bladder and colon. From these animals (n = 5, 2 males, 3 females), 14 convergent bladder/colon afferents were analyzed. Interestingly, these convergent afferents did not cluster into a separated group but were spread across the different clusters.
Unlike colon afferents, the first branch point in clustering of bladder afferents was defined by Gfra3 and Gfra1. Neurons in these groups had either high levels of Gfra1 or Gfra3, but never both. Nefh, the gene associated with the first branch point for colon afferents, was also highly segregated, being expressed at moderate to high levels in all Gfra1 clusters, but moderately in only one Gfra3 cluster. Similar to colon afferents, there were distinct TL and LS clusters as well as mixed clusters.
Thoracolumbar cluster
Only one bladder afferent cluster was composed exclusively of TL afferents (tlPEP; Figs. 3, 4B). The tlPEP cluster had high Calcα, Cd274, Scn10a, Scn11a, and Trka, and moderate levels of Nefh, P2rx3, and Piezo2. This cluster was also distinctive in having both Gfra1 and Gfra3, although at low to moderate levels. The low coexpression of Gfra1 and Gfra3 is likely the reason for positioning of this cluster at the border of the major branch point.
Lumbosacral clusters
Bladder afferents formed three LS exclusive clusters. The lsPEP cluster had low Gfra3 and almost no GFRa1 (undetectable in 19 of 28 cluster and moderate to low in the rest), high Th, and moderate Calca, Tac1, Ret, and Tprv1. The lsNP and lsNF clusters had no cells with detectable Gfra3, but moderate Piezo2. The lsNP and lsNF clusters can be separated by the high levels of Asic2, Gfra1, and Nefh in the lsNF cluster. The lsNP cluster also had high Trka (moderate in the NP vs low in the NF).
Mixed clusters
Bladder afferents formed three mixed clusters, mPEP, mNF, and mNF-PEP. The mPEP cluster was the largest and had the highest level of Gfra3 of any bladder cluster. mPEP also had high Tac1 and Ret, and moderate Calca, Trpa1, Trpv1, Th, and Scn11a. The mixed mNF and mNF-PEP clusters were small and similar to the LS clusters based on Gfra1 and Nefh, but also had high levels of several other genes, including Scn10a, Scn9a, Asic3, and Piezo2. The mNF cluster can be differentiated from the mNF-PEP cluster by high Cd274, P2rx3, Asic2, Trpm3, and Trkc. In contrast, the mNF-PEP cluster had high Calca, Trka, Asic1, and P2ry2.
Comparison of gene clustering for bladder and colon afferents
Comparison of high-level transcripts across colon and bladder clusters shows that many are present in afferents that innervate the two organs. However, different combinations of transcripts define the clusters (Fig. 4A,B). Dendrograms that illustrate the segregation of the 32 genes examined in spinal afferents for both organs (Fig. 5A,B) show some significant differences; the first major branch point for colon is associated with a peptidergic phenotype (e.g., Calca, Tac1, Trka, Trpa1, Trpv1, Gfra3) versus a nonpeptidergic phenotype (e.g., Asic1, Trkb, Gfra2, Mrgprd, Mrgpra3, Trpm8, Piezo2). The clustering in bladder shows a similar first branch point; however, three high-level and closely clustered genes (Oprk1, Cd274, and P2rx3) are shifted from the nonpeptidergic to the peptidergic group. In addition, in colon, Trpv1 nearest neighbors are peptidergic-associated genes (Gfra3 and Trpa1), whereas in the bladder Trpv1 was separated from these genes by three branch points and clustered with Trpm3. Given that the overall percentage of cells expressing all 32 genes was not dramatically different for bladder and colon afferents, these differences suggest that at least some bladder afferents have novel properties not present in colon afferents.
Calcium imaging and validation of gene function in colon afferents
Because RNA expression may not correlate with synthesis of functional protein (Adelman et al., 2019), we performed in vitro calcium imaging of back-labeled colon afferents to examine the relationship between mRNA expression and function. Afferents from different levels of the neuraxis were analyzed. Three agonists (capsaicin, mustard oil, and α,β-methylene ATP) known to activate colon afferents (Christianson et al., 2010; Shinoda et al., 2010) via membrane-bound receptors that were included in this analysis were used. We also examined the response properties of identified afferents to ligands of immune-related receptors (interferon-α, interferon-γ, interleukin-4, and interleukin-8) previously shown to produce calcium transients in primary afferents (Oetjen et al., 2017; Wang et al., 2017). The percentage of responders was compared with the percentage of cells that exhibited detectable (low, medium, or high) levels of gene transcripts (Table 3).
Calcium imaging was performed on 131 colon afferents [56 NG afferents (n = 6, 4 males, 2 females), 35 TL afferents (n = 8, 2 males, 6 females), and 40 LS afferents (n = 8, 4 males, 4 females; Fig. 6A)]. Many TL afferents were positive for Trpa1 (91.67% vs 48.84% for NG afferents and 43.90% for LS afferents), and this correlated with a high percentage of mustard oil-responsive neurons (63.64%, compared with 36.58% for NG neurons and 47.37% for LS neurons). Capsaicin activated 74.28% of NG neurons, 77.78% of TL neurons, and 95.45% of LS neurons, consistent with the relatively broad Trpv1 expression across all levels (Table 3, Fig. 3). Adelman et al. (2019) have shown for cutaneous afferents that neurons responsive to capsaicin had a ΔCt of ≤9 for Trpv1 (i.e., a ΔCt of 9 was the expression threshold). In this study, this level of Trpv1 was met in ∼80% of the colon afferents, approximately matching the percentage of neurons with capsaicin responses (Table 3).
The P2rX2/3- and P2rX3-selective agonist α,β-methylene ATP activated 96.15% of NG neurons, 34.78% of TL neurons, and 34.48% of LS neurons. For the NG neurons, the percentage of α,β-methylene ATP responders (96.15%) matched the uniformly high level of P2rx3 and P2rx2 detected in NG afferents. In contrast, for TL and LS neurons, there were fewer responders to α,β-methylene ATP (34.78% and 34.48%, respectively) relative to the percentage of neurons expressing detectable levels of P2rx3 and P2rx2 (Table 3). Adelman et al. (2019) found that for cutaneous afferents expressing the P2rx3 receptor (P2rx2 was not examined), a ΔCt of 6 was the threshold for the presence of functional receptors. This value was measured in only 20.45% of TL neurons and 14.28% of LS neurons, consistent with the low percentage of afferents responding to α,β-methylene ATP. In terms of the amplitude of the calcium transients in response to α,β-methylene ATP, for neurons that did respond, the peak responses of NG afferents was significantly higher than for LS and TL afferents (F(12,210)= 3.065; p = 0.0204 for TL; p < 0.0001 for LS; Fig. 6B). This suggests that for P2X receptors, there is a direct relationship between the transcript level and the response to ligand.
With the exception of Il4ra in NG and LS neurons, all immune-related receptors assayed were detected in nearly all cells. Unlike other ligands examined, the percentage of neurons activated by the immune receptor ligands was significantly lower than the percentage of cells positive for these receptor genes (except for l4ra in the NG). In NG neurons, calcium responses to interferon-α, interferon-γ, interleukin-8, and interleukin-4 occurred in approximately half (37–58%) of the cells expressing the corresponding receptor mRNA. A lower percentage of TL and LS neurons exhibited calcium transients in response to these ligands (24–47% for LS neurons and 15–34% for TL neurons). These results indicate that mRNA alone cannot be used as an indicator of functional activity for these receptors. This may be because these immune receptors are not ionotropic receptors and often act as receptor complexes that activate complex intracellular signaling cascades (Bezbradica and Medzhitov, 2009). For cells where immune ligands did cause a calcium transient, the amplitude of these transients were similar to those for capsaicin, mustard oil, and α,β-methylene ATP (Fig. 6B).
Discussion
Identification of unique protein or mRNA markers that can be used as surrogates for functional phenotyping has long been a goal of primary afferent biology. Here we used single-cell mRNA analyses of up to 42 genes to determine whether a unique molecular signature for colon and bladder primary afferents mapped to their origin (i.e., NG, or TL or LS ganglia). Using unbiased AHC, we found that for both colon and bladder afferents, most clusters had neurons primarily from a single level of the neuraxis (>85% of cells from NG, or TL or LS ganglia) than from mixed levels (Fig. 7). For the 13 colon clusters, 4 originated in the NG (ngNPa, ngNPb, ngNPc, ngNF), 2 from TL ganglia (tlPEPa and tlPEPb), and 4 from LS ganglia (lsPEP, lsNFa, lsNFb, and lsNF-PEP). Of the seven bladder clusters, one came from TL ganglia (tlPEP) and three came from LS ganglia (lsPEP, lsNP, and lsNF; Fig. 7). For colon, there were also clusters that primarily innervated the proximal (tlPEPa) or distal (mPEPa, lsPEP, lsNFa, and lsNFb) colon. Genes in our panel were chosen based on a validated role in stimulus detection (e.g., TRP channels, piezo2, P2X receptors, ASICs, cytokine receptors) or responsiveness (e.g., opioid receptors, neuropeptides, sodium channels), suggesting that function-specific components of sensory information arriving in the CNS is segregated by sensory ganglia at different levels.
The mRNA analysis from the current study revealed some key differences between colon and bladder afferents compared with studies examining unidentified (i.e., not back-labeled from a specific target) DRG afferents (Usoskin et al., 2015; Zeisel et al., 2018). Almost all spinal colon afferents expressed some level of Calcα, whereas previous studies of unidentified afferents reported clusters with no Calcα expression. In our study, >80% of colon and bladder afferents express Th, whereas in previous studies of unidentified afferents Th expression was restricted to putative C-low-threshold mechanoreceptors. However, studies of unidentified sacral DRG afferents did replicate the findings here with both studies showing high levels of Trpv1, Chrna3, and Tac1 (Smith-Anttila et al., 2020). This may reflect the relatively high proportion of visceral afferents in these ganglia.
Unlike RNAseq analysis, our study of visceral afferents focused on a small number of curated genes. The clusters identified are similar to those detected in a previous RNAseq analysis of back-labeled distal colon afferents from both TL and LS levels (Hockley et al., 2019). Using single-cell consensus clustering (Kiselev et al., 2017) of 314 cells in the work by Hockley et al. (2019), seven clusters were identified based on the presence of peptidergic (Calcα, Tac1, Trpa1) or nonpeptidergic (Mrgprd, P2rx3, Gfrα2) markers, and the expression of Nefh (which was coexpressed with Piezo2). Hockley et al. (2019) also identified mixed clusters arising from a combination of TL and LS afferents and clusters that were exclusively from LS afferents (defined as “pelvic nerve” afferents). The present study identified clusters that arise exclusively from TL afferents that were not identified by Hockley et al. (2019). However, these TL clusters have a majority of afferents originating from the proximal colon, which were not included in the study by Hockley et al. (2019).
The NG clusters were unique in that, unlike TL or LS clusters, NG clusters projected virtually equally to proximal and distal colons (45.13% proximal vs 45.44% distal, compared with a 38.25% proximal vs 52.94% distal for TL ganglia and 2.54% proximal vs 95.35% distal for LS ganglia). This suggests that unlike clusters arising from DRG afferents, innervation of different NG clusters is required throughout the length of the colon. This could reflect a role for NG afferents in neuroimmune monitoring, which would be in line with the high percentage of NG neurons that exhibited calcium transients in response to cytokine application.
RNAseq studies have been performed on NG afferents and compared with jugular afferents (Wang et al., 2017; Kupari et al., 2019; Mazzone et al., 2020). In the mouse, jugular afferent somata are located in a ganglion fused to the placode-derived NG and, like DRG neurons, are neural crest derived (Baker and Schlosser, 2005). These studies found that, like the DRG neurons examined here, jugular neurons are molecularly distinct from NG afferents. However, cluster analysis found several more unique NG clusters than identified in the present study [18 nodose clusters with a major breakdown between “mechanosensitive” groups (1–11) and “nociceptive” clusters (12–18); Kupari et al., 2019]. A likely reason for this difference is that we used back-labeled neurons from colon, whereas in the study by Kupari et al. (2019) clusters from all the thoracic and abdominal organs innervated by the NGs were represented.
A major theory to explain the need for overlapping sensory innervation of the viscera, especially for those organs that receive nodose and spinal input like the colon, is that different pathways relay unique aspects of the sensory experience. For example, in the case of pain arising from abdominal organs (including the colon), spinal afferents have been proposed as the transducer of noxious pain (e.g., sharp, burning), whereas nodose afferents transmit affective aspects of visceral pain (e.g., fear, anxiety, nausea; Berthoud and Neuhuber, 2000; Grundy, 2002; Sengupta, 2009). These distinctions are supported by patient reports of sensations produced by whole vagal nerve stimulation (Sackeim et al., 2001; Ben-Menachem, 2002), but are contradicted by the similarities in firing properties of nodose and spinal afferents innervating the gut (Sengupta et al., 1990; Ozaki and Gebhart, 2001; Yu et al., 2005; Bielefeldt et al., 2006) and the anecdotal reports of noxious gut pain in patients with spinal transections (i.e., patients that lack ascending input from spinal afferents; Yung and Groah, 2001; Finnerup et al., 2008; Levinthal and Bielefeldt, 2012).
Given that in most cases activation of visceral afferents does not produce conscious sensations (Critchley and Harrison, 2013), the data presented here could be used to argue that the driving force behind the evolution of the complex pattern of sensory innervation arises out of the requirements of homeostatic regulation and not because the different afferent populations are “labeled-lines” carrying segregated qualities of visceral sensation. First, the projections of NG and TL clusters overlap anatomically in a manner that suggests these molecularly distinct afferent populations are making unique contributions to homeostatic function; both NG and TL afferents innervate the entire colon, with approximately equal numbers of afferents, providing direct input to their associated parasympathetic (NGs) and sympathetic (TL ganglia) CNS circuits in the spinal cord and brainstem. Second, while it is clear that neurons at NG, and TL and LS ganglia levels have unique transcriptomic profiles with respect to receptor expression (e.g., TRP and purinergic channels, cytokine receptors) that should “tune” them to different aspects of the internal milieu, all clusters, regardless of origin, express receptors required for monitoring of homeostasis-relevant stimuli. Third, many organs, in particular the small and large bowel, contain large numbers of innate immune cells, and their interactions with their local microbiome affect the entire body (Cryan and O'Mahony, 2011; Maynard et al., 2012). The near ubiquitous expression of immune-related genes in all clusters is consistent with the idea that these afferents play a role in the integration of information required for immunologic homeostasis. Moreover, the efferent function of visceral afferents [i.e., their ability to release molecules (e.g., peptides likes CGRP and substance P)] that can modulate immune cells has recently been shown to play a central role in somatic and visceral tissue responses to pathogenic insult (Baral et al., 2018; Lai et al., 2020; Saloman et al., 2020).
The corollary to this hypothesis is that different aspects of conscious visceral sensation (e.g., both noxious and affective aspects of visceral pain) is the result of the integration of information gathered from all sensory ganglia and then processed by different (but likely overlapping) CNS circuits. Five colon afferent fiber types have been identified (mucosal, muscular, serosal, mesenteric, and muscular/mucosal; Brierley et al., 2004, 2005), as well as four classes of bladder afferents (serosal, muscular, muscular/urothelial, and urothelial; Zagorodnyuk et al., 2006; Xu and Gebhart, 2008). For colon, both TL and LS neurons have mucosal, muscular, and serosal afferents, whereas mesenteric afferents are specific to TL afferents and muscular/mucosal afferents are specific to LS afferents. For the bladder, the TL level is composed of primarily serosal and muscular afferents, whereas the LS afferents contain all four subtypes. This suggests that for both bladder and colon, some TL and LS afferents detect similar mechanical stimuli, whereas others code for unique aspects of these stimuli. That neurons from all levels of the neuraxis (including the NGs) have a role in conscious sensation including pain is supported by the data presented here; all of the colon clusters express proteins shown to play important roles and/or are “required” for nociception in both somatic and visceral afferents (Delafoy et al., 2006; Bielefeldt and Davis, 2008; McIlwrath et al., 2009; Kiyatkin et al., 2013). It seems unlikely that noxious pain would be detected by only a subset of these afferents from one or two levels with the remainder transmitting signals specific to non-noxious qualities of pain (i.e., affective aspects). The more parsimonious explanation is that afferents in all clusters, at all levels, are collaborating to provide a comprehensive report of peripheral conditions, and occasionally this information reaches qualitative and quantitative thresholds required for conscious visceral sensation, including pain.
Footnotes
The authors declare no competing financial interests.
This research was supported by National Institute of Neurological Disorders and Stroke Grants NS-023725 and NS-096705 (H.R.K.); National Institute of Arthritis and Musculoskeletal and Skin Diseases Grant AR-069951 (K.M.A.); National Institute of Diabetes and Digestive and Kidney Diseases Grants DK-124955 (B.M.D.) and OT2-OD-023859 (Marthe J. Howard Principal Designer, B.M.D., Principal Investigator). We thank Mr. Christopher Sullivan for expert technical support and mouse husbandry.
- Correspondence should be addressed to Brian M. Davis at bmd1{at}pitt.edu