Experiments were conducted with adult (8–12 weeks) female C57BL/6J mice (Jackson Laboratory, Bangor, ME, USA). Several chronic pain conditions have a higher prevalence in females, and numerous studies have reported higher pain prevalence in the SCI population among female patients (Cardenas et al., 2004). Women also report greater frequency, severity, and longer-lasting pain, as well as neuropathic pain below the level of injury, in comparison to men (Cardenas et al., 2004). The majority of research examines SCI in male rodents and this study will add to what is known in the literature by focusing on female mice (Cardenas et al., 2004). Naïve animals were group-housed; sham and spinal cord injured animals were individually caged. All animals were maintained on a 12:12-h light-dark cycle with a temperature-controlled environment, and given food and water ad libitum. All treatments and testing were approved by the University of Connecticut Health Center Institutional Animal Care and Use Committee.

Animals were anesthetized by inhalation of isoflurane and a 1.0-cm dorsal midline skin sterile incision was made over T8–T11, as per Ma et al. (2001). Connective and muscle tissue were removed to expose the bone from T9–T10, and a laminectomy was performed at the T10 vertebral level. SCI was produced by compression of the vertical plane of the spinal cord using a vessel clip with a closing pressure of 15 g (WPI, Sarasota, FL, USA) for 1 min, exerting pressure from side to side on the spinal cord. Sudden impact is produced by the rapid release of the vessel clip (Tator and Poon, 2008). After removal of the clip a hemorrhagic ring at the site of compression is present. This injury is analogous to the majority of lesions in humans, as the model constitutes both contusion and compression (Marques et al., 2014). The injury resulted in paralysis of the hindpaws below the level of injury and was considered mild-moderate based on the ability of animals to exhibit slight ankle movement 24 h following injury (a score of 1 on the Basso Mouse Scale for locomotion), their ability to recover additional ankle movement by day 4, but the inability of mice to recover complete bladder control by day 4 (Basso et al., 2006). The wound is closed with coated vicryl absorbable sutures (Ethicon, Somerville, NJ, USA). Mice were allowed to recover in warm cages for 24 h. All animals were administered antibiotics once immediately following surgery (5 mg/kg gentamicin), as well as subcutaneous saline for 4 days following surgery, without analgesics. Sham control mice received the same treatment excluding the vessel clip. Manual bladder expression on SCI mice was performed twice daily until mice were sacrificed. Mortality was less than 10% and typically occurred during the laminectomy due to excessive blood loss, or when postinjury weight loss required the animal to be sacrificed. Spinal cords collected at 4 days post-injury were most notably characterized by minimal cavitation and scar tissue that progressively diminished with increasing distance distal to the site of injury (~L1–2 vertebral level), analogous to previous studies that have characterized this injury model (Joshi and Fehlings, 2002; Marques et al., 2014).

Female C57BL/6 mice (6–10 weeks old) were fed powdered milled chow mixed to contain a final concentration of 0.2% bis (cyclohexanone) oxaldihydrazone (cuprizone; Sigma-Aldrich, St. Louis, MO, USA), with food and water available ad libitum. Each mouse received approximately 5 g of chow per day, fresh cuprizone containing chow was prepared every 7 days. Cuprizone feeding was maintained for 35 days, and tissue was collected for protein analyses on day 35. Digitized, non-overlapping electron micrographs of the corpus callosum were analyzed for unmyelinated axon frequency and g-ratios to assess effectiveness of cuprizone treatment (Wasko et al., 2019; N = 5).

Spinal cord segments at the level of laminectomy (T8–T11) were collected from naïve, sham, SCI, or cuprizone-treated mice immediately following perfusion with ice-cold 0.9% NaCl. Spinal cord segments were homogenized in ice-cold buffer containing 20 mM TES, pH 7.4, 10 mM mannitol, 0.3 mg/ml phenylmethylsulfonyl fluoride, 2 μg/ml leupeptin, 2 μg/ml pepstatin, 2 μg/ml benzamidine, 16 μg/ml benzamidine, and 50 μg/ml lima bean trypsin inhibitor at a concentration of 0.1 g tissue per 1 ml buffer (Mains et al., 2018). Homogenates were freeze-thawed three times, centrifuged (20 min, 17,400 g), and supernatants were collected. Approximately 60 μg of protein per sample was used for each ELISA. The ELISA assays were performed according to the manufacturer’s instructions (R&D systems mouse duo-sets IL-10, IL-6, IL-1β, TNF-α, completed with Ancillary Reagent Kit 2, Minneapolis, MN, USA). The sample absorbance was read with an ELISA plate reader at 450 nm; readings were also taken at 570 nm to subtract optical background. The concentration was determined based on a standard curve. All results were normalized to amount of protein added per sample and graphed as pg/mg. Naïve, 1 day Sham, 4 days Sham, and 1 day SCI N = 4 mice; 5 days and 7 days Sham conditions N = 3 mice.

To backlabel DRG L2-L6 projecting to the hairy hindpaw skin, mice were anesthetized with isoflurane. 0.3% wheat germ agglutinin conjugated to an AF-488 dye (WGA-488, Thermo Fisher Scientific, Waltham, MA, USA) in sterile PBS was injected into the sural, common peroneal, and saphenous nerve skin territories for retrograde labeling of DRG neurons (da Silva Serra et al., 2016; Berta et al., 2017). A total of 6 μl of WGA-488 was injected 2 days prior to surgery by three 2 μl injections in the lateral zones of each hindpaw (2 μl per nerve territory) using a 10 μl Hamilton Syringe and 30G needle. This was performed on both hindpaws of each mouse. This technique does not cause significant injury to the sensory afferents being studied.

Mice were anesthetized 4 days post-surgery with an intraperitoneal injection of ketamine (100 mg/kg) plus xylazine (10 mg/kg) and perfused with ice-cold 0.9% NaCl. A laminectomy was performed and L2-L6 DRG from both sides of the spinal column were collected into cold HBSS (KCl 5.4 mM, NaCl 137 mM, Glucose 5.6 mM, Hepes 20 mM, pH 7.35 NaOH), after which the mice were sacrificed by decapitation. Sensory neuron dissociation was performed as described (Malin et al., 2007). Briefly, following collection, tissue was treated with 60U papain (Worthington), 1 mg of cysteine, and 6 μl of NaHCO₃ in 1.5 ml HBSS at 37°C for 10 min. Tissue was then treated with 12 mg collagenase II (Worthington, Lakewood, NJ, USA) and 14 mg dispase (Roche, Basel, Switzerland) in 3 ml HBSS at 37°C for 20 min, washed, and triturated with fire-polished glass Pasteur pipettes in 1 ml of DMEM (Gibco Thermo Fisher Scientific, Waltham, MA, USA) supplemented with FBS (Hyclone, Logan, UT, USA) and pen/strep (Gibco). The cell suspension was pelleted (1 min, 80 g), DMEM was removed, and cells were re-suspended in a modified solution (Citri et al., 2011) containing 140 mM NaCl, 5 mM KCl, 10 mM Hepes, 10 mM glucose, 0.1% Bovine Serum Albumin, pH 7.4. After re-suspension, cells were strained through a 70 μm cell strainer and placed on ice in the modified solution until fluorescence-activated cell sorting (FACS).

Single cell suspensions of cells isolated from in situ WGA-488 labeled DRG were live-stained using Hoechst 33342 (10 μg/ml, Thermo Fisher Scientific, Waltham, MA, USA) and propidium iodide (PI; 1 μg/ml), and analyzed on an Amnis ImageStreamX Mark II imaging flow cytometer (Luminex Company, Austin, TX, USA). Fluorescent cell images were captured using a 60× objective lens with excitation from a 405 nm laser at 20 mW power and a 488 nm laser at 200 mW power. Images of in-focus nucleated WGA AF488-positive cells were identified and electronically gated using IDEAS software (Amnis, v6.2.183, Seattle, WA, USA).

Neurons labeled with WGA-488 dye in situ in the DRG were purified by FACS 4 days post-surgery. Following primary dissociation of DRGs L2-L6, single cell suspensions were analyzed and sorted using a BD FACS Aria II cell sorter (Becton Dickinson) set up with a 130 μm nozzle at 12 PSI in order to gently isolate cells between 10 and 30 μm. Single live neurons were defined by electronic gating in FACS DIVA software (BD, ver. 8.01) using forward and side-angle light scatter, omission of PI (1 μg/ml), and AF488 fluorescence. All fluorescence gates were confirmed using fluorescence minus one controls (e.g., a sample of cells from unlabeled DRG was used to gate for AF488 positive cells and a sample of cells not treated with PI was used to set the live cell gate). WGA-488 positive cells were sorted directly into lysis buffer (NucleoSpin RNA XS Kit, Machery-Nagel, Bethlehem, PA, USA) and immediately placed on dry ice until RNA extraction. Cell sorting took 5–10 min per sample and was performed at 4°C before placement of samples on dry ice.

RNA from FACS sorted cells was isolated using NucleoSpin RNA XS Kit, including a DNA digestion step but without carrier RNA step. Before library preparation, RNA quality and integrity was tested for each sample using the Agilent High Sensitivity RNA Screen Tape on the Agilent Tapestation 2200 (Agilent Technologies, Santa Clara, CA, USA). RNA with RIN values ≥6.7 (minimum 6.7, maximum 9.9, average 7.3) was further processed for RNA sequencing. It is important to note that these RIN values are based on relatively low RNA sample concentrations (RNA extracted from approximately 3,000 single cells per sample) and that the Agilent Tapestation software may have issues accurately calling RIN scores from sample concentrations below 50 ng/μl (Schroeder et al., 2006; Mueller et al., 2016). Because of this, RIN scores as low as 6.7 were considered suitable for RNA sequencing based on additional comparisons from the gel electrophoresis run on the Agilent Tapestation software that showed little to no RNA degradation, and previous work that has shown RNA samples with RIN scores below 7 are sufficient to identify differentially expressed genes (Gallego Romero et al., 2014). Library preparation was performed using the Illumina sequencing kit for high output 75-cycles for 25–30 M total single-end reads per sample. DESeq2 analyses¹ of differential expression were performed, and outliers beyond 30%–50% of the mean for each group of animals were eliminated (Love et al., 2014; Conesa et al., 2016; Labaj and Kreil, 2016; Wu and Wu, 2016).

Data were analyzed by Ingenuity Pathway Analysis (IPA; Qiagen, Germantown, MD, USA). An overlap of significance for DEseq2 comparisons plus an RPKM cut-off >10 were required for transcripts to be included for IPA analysis. We analyzed 125 transcripts for comparisons between SCI and naïve groups, and 560 transcripts for comparisons between SCI and sham groups.

Differences between groups were compared using student’s t-test or ANOVA, followed by Tukey’s post-test, Bonferroni’s multiple comparisons test, or by unpaired student’s t-test. P-values < 0.05 were considered statistically significant. Statistics on PCR data were conducted using delta CT values. Heat maps were generated by Microsoft software (Excel), hierarchical clustering was generated by Gene Cluster 3.0 and visualized using Java TreeView (de Hoon et al., 2004; Baek et al., 2017). All other data were plotted using Prism 6 (GraphPad Software, San Diego, CA, USA). R studio was utilized for differential expression analysis; Prism software was used for all other statistical tests.

To determine an optimal time point to observe transcriptional changes contributing to the transition from acute to chronic pain, we tested behavioral differences between naïve and sham-operated mice. The objective of this was to ensure that changes within the DRG were due to injury to the spinal cord itself, not to the laminectomy performed in both injured and sham mice. It is likely that early effects of sham surgery involve changes in gene expression that overlap with a subset of genes that also change following SCI. It is important to define how long the functional effects of sham surgery persist to better discern the effects of SCI. Specifically, we tested naïve and sham mice for open field locomotor differences 1, 3, 5, and 7 days post-surgery (Figure 1A1). Naïve and sham mice did not differ significantly at any of the time points tested, including as early as 1 day post-surgery. To ensure that SCI mice (T10 compression-clip injury) did exhibit behavioral differences and locomotor deficits following injury, we compared SCI mice to naïve and sham mice for open field behavior 1 day post-injury (Figure 1A2). SCI mice exhibited substantially decreased total ambulation (one-way ANOVA, p = 0.0059, Tukey’s multiple comparisons test, naïve **p < 0.005, sham *p < 0.05) following injury, as expected, although time spent in the periphery did not differ (Supplementary Figure S1A). Additional tests for mechanical (von Frey) and thermal (hot water tail-flick) hypersensitivity showed that naïve and sham conditions did not differ significantly at any time point (Figures 1B1,B2). Previous literature has demonstrated analogous results when testing hindpaw thermal sensitivity by radiant heat stimulus on sham animals at multiple time points (Christensen and Hulsebosch, 1997; Bedi et al., 2010; Gaudet et al., 2017). We did not examine SCI mice for thermal or mechanical sensitivity, as paralysis below the level of injury prohibited below-level sensitivity testing at the time points studied. Previous work using the clip-compression model has shown varying results regarding the occurrence of above-level sensitivity following SCI (Bruce et al., 2002; Densmore et al., 2010). While above-level pain may occur following injury, below-level pain is one of the most prevalent, and most difficult, pain types to treat, and is the focus of these studies (Yezierski, 2005). The development of chronic pain at later time points following clip-compression models of SCI (in addition to several other models) has already been well characterized. Therefore, we focused on early time points after injury, to identify contributions to the onset of pain, rather than changes after chronic pain has already developed (Bruce et al., 2002; Nakae et al., 2011; Gaudet et al., 2017).

Previous work has highlighted the importance of inflammatory cytokines within the CNS to facilitate the transduction of noxious stimuli in neuropathic pain (Cook et al., 2018). Thus, we examined post-surgical changes in inflammation, using ELISAs to analyze common cytokine markers to ensure that sham mice did not differ significantly from naïve mice following the removal of bone and muscle. We did not observe significant differences in inflammatory cytokine levels (TNF-α, IL-6, IL-1β, IL-10) in extracts of spinal cord segments (T8-T11) among naïve or sham mice at 1, 4, 5, and 7 days post-surgery (Figure 1C, Supplementary Figure S1B), similar to previous studies (Stammers et al., 2012). Spinal cords from SCI mice had significantly increased levels of IL-6 and IL-10 in comparison to naïve controls (one-way ANOVA, Bonferroni’s multiple comparisons test; *p < 0.05, **p < 0.005, respectively). When compared to sham conditions, spinal cords from SCI mice had significantly increased levels of IL-10 in comparison to 4 days, 5 days, and 7 days sham conditions (one-way ANOVA, Bonferroni’s multiple comparisons test; ***p < 0.001, **p < 0.005, **p < 0.005, respectively.

As a positive control, we tested cuprizone-treated mice (a model of multiple sclerosis), which are known to secrete proinflammatory cytokines at the end of 5 weeks of treatment (Mukhamedshina et al., 2017). Cuprizone-treated mice showed expected increases in levels of TNF-α, IL-1β, and IL-10 in comparison to naïve controls (one-way ANOVA, Bonferroni’s multiple comparisons test; *p < 0.05, **p < 0.005, **p < 0.005, respectively; Schmitz and Chew, 2008). Similarly, in comparison to sham mice, spinal cord extracts from cuprizone-treated mice had significant increases in IL-1β at the 1 day sham timepoint, and IL-10 at the 4 days, 5 days, and 7 days sham timepoints (one-way ANOVA, Bonferroni’s multiple comparisons test; *p < 0.05, ****p < 0.0001, **p < 0.005, ***p < 0.001, respectively). Behavioral testing and cytokine analyses did not reveal differences between naïve and sham mice at any time point.

Because the laminectomy did not significantly affect behavior or inflammatory responses in sham mice, we determined an optimal time point to study the transition of acute to chronic pain based upon established characteristics of nociceptors following SCI. The inflammatory cytokine assays demonstrated significant increases in inflammation 1 day post-SCI compared to 4 days sham controls (Figure 1C, Supplementary Figure S1B). We did not observe significant differences in inflammatory cytokine levels (TNF-α, IL-6, IL-1β, IL-10) in extracts of spinal cord segments (T8-T11) among naïve or sham mice at 1, 4, 5, and 7 days post-surgery (Figure 1C, Supplementary Figure S1B); this is similar to previous studies (Stammers et al., 2012). In addition, previous studies have documented onset of spontaneous activity in nociceptors distal to the site of SCI as early as 3 days post-injury in rats; increased activity persisted for at least 8 months (Bedi et al., 2010). Thus, it is possible that the transition to chronic pain begins around 3 days after SCI in rats and mice, and this transition is not due to laminectomy, as nociceptors isolated from sham rats show either no significant changes in spontaneous activity (Bedi et al., 2010; Yang et al., 2014) or only a minor increase compared to the effect of SCI (Odem et al., 2018). Subsequently, we chose 4 days post-injury to assess transcriptional changes in nociceptors below the level of injury in mice.

To perform transcriptional profiling on nociceptors that project to the cutaneous skin after injury, we injected wheat germ agglutinin conjugated to an AF-488 dye (WGA-488) into the sural, common peroneal, and saphenous nerve skin territories for retrograde labeling of DRG neurons. Next, we performed compression-clip SCI or sham surgeries at the T10 vertebral level 2 days post-WGA injection and collected L2-L6 DRG for dissociation 4 days post-injury (Figure 2A). Based on the vertebral level at which we performed SCI, the DRG collected were located below the level of injury, thus our analysis was comprised of non-injured nociceptors that do not project directly to the lesion site. Flow cytometry confirmed that our cell population of interest (cutaneous nociceptors) was positively labeled with WGA-488 (Figure 2B). We also observed non-labeled (WGA-488) cells and dead cells (propidium iodide, PI+), to be excluded from cell sorting and analysis (Figure 2B). As expected, a significant number of viable small nociceptor cells were excluded in this analytical approach, to avoid including RNAs from non-nociceptor cells or attached fragments of dead cells (examples in Figure 2C; Thakur et al., 2014; Lopes et al., 2017; Megat et al., 2019).

We performed FACS purification of nociceptor populations from naïve, sham, and SCI adult (8–12 weeks old) female mice (n = 5 per condition). We pooled DRGs from lumbar regions L2-L6 on either side of the spinal column to ensure all DRGs isolated would have projections to the hairy hindpaw skin. DRG cells were enzymatically dissociated and subjected to flow cytometry, to gently isolate positively labeled cells between 10 and 30 μm; PI staining was used to identify dead cells. All conditions were gated on DRG from naïve mice that did not receive WGA-488 injections, enabling purification of positively labeled cells (Figure 3A). Analysis of our flow cytometry data shows that, we were successful in retrogradely labeling DRG neurons projecting to hairy hindpaw skin. Many positively labeled neurons were part of cell aggregates, limiting retrieval of the single-cell population of interest to ~2% of all dissociated cells per animal (Figure 3B). This percentage amounted to approximately 3,000 cells per mouse, a suitable representation of this cell population in agreement with previous studies (Shields et al., 2012; Goswami et al., 2014; Thakur et al., 2014; Usoskin et al., 2015; Berta et al., 2017). DRG populations were sorted directly into lysis buffer and placed on dry ice to preserve transcriptional profiles at the time of isolation. RNA quality was tested using the Agilent TapeStation; a representative image of an RNA sample following FACS purification is shown (Figure 3C).

We used the RNA sequencing data to evaluate the neuronal population that had been isolated. Scn10a, which encodes Na_v1.8, is present in 80–90% of nociceptors and Trpv1 serves as a marker for the peptidergic population of nociceptors (Basbaum et al., 2009; Harriott and Gold, 2009; Wu et al., 2013). RPKM values of 1,000 for Scn10a and 400 for Trpv1 confirm that FACS purified cells express high levels of these nociceptor markers Figure 4A). The low RPKM values observed for Parvalbumin (Pvalb), a glial transcript also found in large diameter proprioceptors and Aβ neurons, and glial fibrillary acidic protein (Gfap), another marker for glial cells, suggests that the non-nociceptive sensory neurons responsible for touch and proprioception (Aβ and Aδ neurons) and satellite glial cells are largely absent from our purified nociceptor population (Figure 4B; Huang et al., 2013; Le Pichon and Chesler, 2014); this conclusion is further supported by the somatosensory profiling of transcripts in Figure 4C (below). Previous studies have confirmed that intact DRG as well as unsorted dissociated DRG yield much higher levels of non-neuronal markers (Thakur et al., 2014). This highlights the importance of excluding cell aggregates during cell sorting (Figure 3) to avoid the analysis of transcripts from cells outside of the target population of nociceptors.

We next analyzed gene expression patterns for known functional mediators of somatosensation (Chiu et al., 2014; Le Pichon and Chesler, 2014; Usoskin et al., 2015). The purified cutaneous nociceptors displayed high expression levels of genes involved in thermosensation and nociception, such as specific Trp channels (notably Trpv1), sodium channels (Scn9a, 10a, 11a) and peripherin (Prph; Figure 4C). The skin is innervated by both peptidergic nociceptors and non-peptidergic neurons (Zylka et al., 2005; Lallemend and Ernfors, 2012; Kjell and Olson, 2016). Markers for non-peptidergic nociceptors were abundant, such as Mas-Related G-Protein Coupled Receptor Member D (Mrgprd), Runt related transcription factor 1 (Runx1), Ret proto-oncogene (Ret), and transient receptor potential cation channel C3 (Trpc3). As expected, transcripts enriched in peptidergic nociceptors were present, such as Calca and Calcitonin Related Polypeptides (Calcb) and Tac1 (the tachykinin precursor for peptides such as Substance P), the peptide processing enzyme peptidylglycine α-amidating monooxygenase (PAM), as well as Npy1r, one of the most abundant Npy receptors (Julius and Basbaum, 2001; Basbaum et al., 2009). However, genes encoding proteins involved in itch, such as brain natriuretic peptide (Nppb) and histamine receptor H1 (Hrh1) were only expressed at low levels (Usoskin et al., 2015). Similarly, genes responsible for proteins involved in tactile function, including Transient Receptor Potential cation channel subfamily C Member 1 (Trpc1), and those responsible for proprioception, such as Runt related transcription factor 3 (Runx3), exhibited low expression levels (Le Pichon and Chesler, 2014). In contrast, nerve growth factor receptors (Neurotrophic Receptor Tyrosine Kinases 1–3, Ntrk), which are characteristic of many subtypes of nociceptors, were expressed at high levels for each condition and did not differ significantly when comparisons were made between all conditions. It is of significance that Ntrk2 is typically used as a marker to identify myelinated sensory neurons, and relatively high levels of Ntrk2 (and Ntrk3) may indicate the presence of some proprioceptive or low threshold mechanoreceptors within the isolated population of neurons (Lallemend and Ernfors, 2012; Usoskin et al., 2015). However, this seems unlikely based on the absence of other common proprioceptor and mechanoreceptor transcripts within our dataset (Figure 4C). It is possible that higher levels of Ntrk2 and Ntrk3 are unique to our isolated population of neurons that project to the hairy hindpaw skin. High expression levels were also observed for some of the transcripts of the two major subpopulations of nociceptors: peptidergic and non-peptidergic (Figure 4D).

To further assess expression profiles of the purified nociceptor population and differences among naïve, sham, and SCI conditions within this nociceptor-enriched population, we focused on expression patterns of gene families that mediate general neuronal functions (Chiu et al., 2014; Berta et al., 2017). We used differential expression analysis (DESeq2) to analyze significant changes between SCI and naïve or SCI and sham populations (Love et al., 2014). Pairwise comparisons of significant genes generated by DESeq2 analysis yielded many differentially expressed genes in each subset (Supplementary Figures S2A,B). We also assessed differences between sham and naïve populations. Any genes that differed significantly between the sham and naïve populations were considered significant transcript changes due to laminectomy, and thus were excluded from analysis in the SCI vs. Naïve and SCI vs. Sham comparisons (Supplementary Figure S2C).

We focused on expression patterns of gene families which mediate neuronal functions and contribute to pain phenotypes and found both high expression levels and significant differences within the chloride channel family, Trp channels, glutamate receptors, GABA receptors, potassium channels, sodium channels, and piezo channels (Figures 5A–G, p-values in Table 1). We also examined ASICs (acid-sensing ion channels), calcium channels, glycine receptors, and P2rx and P2ry families (purinergic receptors), because these are widely studied gene families known to be involved in the development or maintenance of chronic pain (Supplementary Figures S3A–E). Many channels and receptors were highly expressed within this cell population, but no significant changes among the transcripts were demonstrated for any condition (Supplementary Figure S3).

Several transcription factors were highly expressed, including signal transducer and activator of transcription 3 (Stat3), FJB osteosarcoma oncogene (Fos), and Jun Proto-Oncogene (Jun; Supplementary Figure S4). Evidence from various models of neuropathic pain implicates these transcription factors in the development or maintenance of chronic pain (Naranjo et al., 1991; Harris, 1998; Dominguez et al., 2008; Tsuda et al., 2011; Xue et al., 2014). Stat3 inhibitors are used to treat peripheral nerve injury-induced hyperexcitability within dorsal horn neurons, pain behaviors, chronic constriction injury, and signaling of IL-6 cytokines (Dominguez et al., 2008; Tsuda et al., 2011; Xue et al., 2014). Previous studies also show that Fos links extracellular events to long-term intracellular changes (such as noxious stimuli) and have established Fos expression as a valid tool to study nociceptive changes (Naranjo et al., 1991; Harris, 1998). Jun also contributes to persistent pain phenotypes following injury (Naranjo et al., 1991). DESeq2 analysis determined additional transcription factors to be significantly altered by SCI (Supplementary Figure S5). Interestingly, activating transcription factor-3 (Atf3), a cellular stress-inducible factor that also increases in non-injured neurons after axotomy, was present at similar levels as Atf4 in all conditions but did not differ significantly. This suggests that Atf3 may be upregulated by axonal injury during the DRG dissociation process, similar to other findings, but that the increases in Atf3 are not due to SCI itself (Tsujino et al., 2000; Huang et al., 2006; Seijffers et al., 2007).

To gain further insight into differentially regulated genes in our isolated population of neurons, we compared our dataset with similar studies on isolated DRGs from publicly available datasets which focused on nociceptors (expressing Na_v1.8) or DRG collected from the same lumbar spinal level (L2–L6; Figure 6; Thakur et al., 2014; Usoskin et al., 2015; Hu et al., 2016; Megat et al., 2019). Unsupervised hierarchical clustering of the top 260 genes revealed that a large number of genes display distinct patterns of expression dependent upon the technique used: isolated neurons from all DRG by translating ribosome affinity purification (TRAP) using the Na_v1.8^Cre mouse (Megat et al., 2019); single-cell isolation from L3 to L5 DRG (Hu et al., 2016); single-cell isolation from L4-L6 DRG (Usoskin et al., 2015); and magnetic cell sorting (MACS) using the Na_v1.8 TdTomato mouse (Thakur et al., 2014; Figure 6A). We found that, while our cutaneous nociceptor-enriched population clustered most closely with the datasets of both TRAP sorted and unsorted DRG from Megat et al. (2019), our isolated population was characterized by its own unique dataset. Notably, datasets from studies analyzing single-cell transcriptomes cluster together, while studies utilizing either TRAP or MACS cluster with their own respective unsorted DRG controls, suggesting that while techniques for isolating specific cell populations can be useful, some variations in gene expression may be attributed to individual differences in cell isolation and RNA extraction methods.

To better characterize enriched genes within our population, we graphed the expression of the 25 most significantly enriched genes in our dataset in comparison to the same datasets used for hierarchical clustering (Figure 6B). While many of the genes have similarly high expression levels, several genes were unique to our isolated population, including diacylglycerol kinase (Dgkh), ankyrin 2 (Ank2), PHD finger protein 24 (Phf24), serine/arginine repetitive matrix 2 (Srrm2), fatty acid synthase (Fasn), phosphoinositide-interacting regulator of transient receptor potential channels (Pirt), and pleckstrin homology domain containing, family A member 6 (Plekha6). We also compared the expression pattern of our naïve nociceptor-enriched population to the various datasets by again examining known neuronal markers of somatosensation (Figure 6C). As predicted, all of the sorted populations exhibit relatively high expression levels of gene transcripts associated with thermosensation, nociception, or neurotrophic receptors, and comparatively low levels of genes associated with itch, tactile function, or proprioception. It is evident that our cutaneous, nociceptor-enriched population can be defined by its distinct gene expression patterns, in particular, high expression levels of the Trp family of genes listed in Figure 6C, as well as Ntrk2, Ntrk3, and Gfra3 neurotrophic receptors.

Based on DESeq2 analysis, levels of several hundred transcripts in the nociceptor-enriched population were altered by SCI. We restricted IPA input lists to genes that had RPKM values greater than 10 and were statistically different (SCI vs. Naïve or SCI vs. Sham) by DESeq2 analysis (Figures 7A,B). Transcripts that exhibited large-fold changes included Mrgprb5, Hal, Chrnb4, Cap2, Sez6l, Calb1, Prokr2, Rxfp1, Nxpe2, and Arap3 (Figures 7A,B). Despite the noticeable changes in these individual genes, they did not appear in any common significant canonical pathways. Instead, we chose to focus on smaller-fold changes that may result in the coordinated functions of many genes. IPA identified several pathways that are considered important for inflammatory processes, pain transduction, or the maintenance of chronic pain (Figure 7C). This includes calcium signaling, Cxcr4 signaling, neuropathic pain signaling in dorsal horn neurons, opioid signaling, purinergic receptor signaling and synaptic long term potentiation (Figure 7C; Julius and Basbaum, 2001; Walters, 2012, 2018). The current study focused on significant changes due to SCI pain and not due to post-surgical pain (i.e., changes between sham vs. naïve groups).

We used qPCR to confirm changes in transcripts of interest from our RNAseq dataset. To validate qPCR and RNAseq comparisons, we compared qPCR or RPKM Log₂ transcript levels of SCI and sham genes of interest (Figure 7D). SCI and sham qPCR fold change results were analogous to the RNAseq fold change dataset. We focused on the synaptogenesis pathway in particular, as it includes several genes present in overlapping canonical pathways, including changes in receptors involved in organization of excitatory signaling (Ephb) and synapses which may be involved in development of chronic pain (TrkB and BDNF). We validated receptors for significant genes in the synaptogenesis pathway, and genes considered possible contributors to pain that also showed significant differences between conditions Gabrg3, Il6st, Kcng3, Piezo2, Scn5a, Trpc3; Waxman et al., 1999; Wickenden, 2002; De Jongh et al., 2003; Devor, 2006; Eijkelkamp et al., 2013; Guptarak et al., 2013; Xia et al., 2015; Deng et al., 2018; Szczot et al., 2018). Additional targets were chosen in order to validate isolation of the correct cell population (Scn10a). Samples for qPCR were collected by backlabeling cutaneous afferents and cell sorting, consistent with samples generated for RNAseq analysis. Following cDNA synthesis, samples were subjected to gene target-specific preamplification using the same primers used for qPCR (Supplementary Figure S6). Hprt was chosen as the housekeeping gene for qPCR analysis because it was the most constant normalizer transcript from the RPKM data across all 15 mice in the present study (Klenke et al., 2016; Lima et al., 2016).

We used IPA upstream network analysis to further interpret the function of the several hundred transcripts significantly altered determined by DESeq2. This method predicted several transcriptional regulators associated with altered expression levels of downstream target genes following SCI (Figure 7E). The upstream regulator molecules graphed do not show a significant change in RNA expression in response to injury. However, these targets are activated by posttranslational modifications that can alter many of the downstream molecules within its network.

We also considered post-surgical inflammation, using cytokine ELISAs to assess changes in the spinal cord at the level of laminectomy (T8–T11). We wanted to assess changes in cytokine levels for two main reasons; first, to confirm that sham mice did not exhibit differences from naïve mice at key time points in comparison to SCI mice, and second to determine whether sham mice exhibited a prolonged inflammatory response, which could potentially be correlated to the development of chronic pain (Guptarak et al., 2013; Krames, 2014). Both pro-inflammatory cytokines TNF-α and IL-1β have been studied in neuroprotection models of SCI. IL-6 has been implicated in neurodegeneration after CNS injury, and the anti-inflammatory cytokine IL-10 exhibits neuroprotective effects (Donnelly and Popovich, 2008; Schmitz and Chew, 2008; Zhang et al., 2019). However, we did not find any significant cytokine changes between naïve and sham mice within 7 days of injury, indicating that the laminectomy did not produce a significant inflammatory response at the time points tested. As a positive control for the cytokine ELISAs (Suzuki and Kikkawa, 1969), we used the cuprizone model for multiple sclerosis. Key pathological features of the treatment include secretion of proinflammatory cytokines such as TNF-α and IL-1β (Schmitz and Chew, 2008). Consistent with previous findings, mice treated with cuprizone exhibited significant increases in TNF-α, IL-1β, and IL-10 relative to naïve mice (Figure 1C; Schmitz and Chew, 2008).

Tissue injury can also lead to prolonged functional changes and hyperalgesia that are accompanied by behavioral changes due to increased spontaneous activity of nociceptors (Carlton et al., 2009; Bedi et al., 2010; Walters, 2012). Spontaneous activity in nociceptors following SCI begins at 3 days after injury and persists for at least 8 months (Bedi et al., 2010). This increase in nociceptor activity elicits changes within the spinal dorsal horn, which receives input from these nociceptors, ultimately contributing to spontaneous pain (Dubner and Ruda, 1992; Wu et al., 2001). However, the source of hyperexcitability of nociceptors after injury is still unknown. Because spontaneous activity in nociceptors begins by 3 days post-injury, and has been correlated with the generation of persistent pain, we chose to observe transcriptomic changes immediately after this time point at 4 days post-injury (Xie et al., 2005). We intentionally excluded large sensory afferents from our experimental model, as the response to SCI by these afferents is transient and large sensory afferents have not been directly correlated with pain transduction in other models of chronic pain (Huang et al., 2006; Hu et al., 2016). We focused on an anatomically defined population of nociceptors (projecting from below the level of SCI to hairy hindpaw skin) by backlabeling from peripheral afferent terminals and sorting based on both fluorescence and size. These neurons represent ~10% of dissociated DRG tissue and have a distinct transcriptome (Thakur et al., 2014). Our approach to identifying anatomically defined small nociceptors is distinct from single-cell transcriptome isolation based on the expression of the sodium channel Scn10a (Na_v1.8) or expression of advillin (Avil, a marker for all neural crest neurons; Lopes et al., 2017; Megat et al., 2019).

Different types of sensory neurons are distinct in their responses to injury. It is likely that, even within an identified subpopulation, cells will nonetheless exhibit heterogeneity (Hu et al., 2016). Because injury does not impact all afferents in the same way, we analyzed gene expression changes within the population of nociceptors projecting to the skin below the level of injury. By focusing on small nociceptors innervating dermatomes below the level of injury, we begin to address what unique set of genes within a specific population may be contributing to the burning, stabbing, and shooting pain reported in SCI patients suffering from below-level neuropathic pain (Siddall et al., 2003). Our goal was to better understand how SCI affects molecular changes within a specific population of neurons, and how this may contribute to hypersensitivity following SCI. We focused our analysis on sensory neurons from lumbar DRGs (below the level of the SCI) projecting to the hairy hindpaw skin (Figure 2A). After confirming, we had isolated the cell population of interest (Figures 2B, 3), we used RNAseq to identify changes in gene expression. The use of RNA-Seq has clear advantages over microarrays, since RNA-Seq is not limited to a set of pre-determined transcripts, has a larger dynamic range of transcript expression, and is highly reproducible (Usoskin et al., 2015). By utilizing this technology, we were able to identify transcript changes undetectable with traditional RT-PCR or microarrays (Wang and Zylka, 2009). Our RNAseq data from naïve, sham, and injured animals display distinct patterns of somatosensory genes present in this nociceptor-enriched population. In particular, RPKM values show high levels of Scn10a (a marker for nociceptors), purinergic receptor P2rx3, Mrgprd (markers for the non-peptidergic population of nociceptors), and Calca and Calcb (neuropeptide precursors), indicating that we isolated the desired nociceptor specific cell population, and also indicating important genes within the population of sensory neurons projecting to the hairy hindpaw skin (Figures 4A–D). Multiple gene transcripts important for itch, tactile function, and proprioception all had relatively low RPKM values, indicating again that we isolated the desired target cell population, and that injury did not induce modifications in the type of stimuli nociceptors transduce (Figure 4C). Our population-level analysis revealed significant changes after SCI in a number of ion channels and receptors that are already known to play a role in pain or hypersensitivity, such as Piezo2, and transcripts involved in excitatory signaling, such as Grik1 (Figure 5C, Table 1). However, there were also many genes whose expression and functional roles in persistent pain have yet to be characterized, including Trpc4 and Ttyh1 (Figures 5A,B, Table 1).

Many studies have utilized RNA-seq technology to gain insight into nociceptor transcriptomes within the DRG and how it changes relative to different pain models (Thakur et al., 2014; Usoskin et al., 2015; Hu et al., 2016; Megat et al., 2019). Although previous screens have yielded various nociceptor-specific genes, we have identified a unique pattern of gene expression within a population of nociceptors projecting to the periphery (Figure 6). By cross-comparing our dataset with similar studies, we were able to identify our cells of interest and confirm that, we obtained cell-type specificity. Notably, many of the current technologies require a transgenic mouse with a cell-type specific reporter. However, by taking advantage of the anatomical organization of the mouse, we were able to use backlabeling and FACS sorting in non-transgenic animals to isolate a nociceptor-enriched population with a low degree of contamination with other cell types. We further implemented this methodology to identify candidate genes in a specific population of neurons below the level of injury in a model of SCI to better understand the heterogeneous injury response among the many subtypes of DRG neurons.

Among these comparisons, we observed large fold changes in several genes (Figures 7A,B, highlighted in red and blue). However, instead of focusing on larger changes in a small subset of genes with individual functions, we concentrated our analysis on the interaction of many transcripts that were significantly altered after SCI and how these influenced intracellular signaling pathways (Figure 7C). Ingenuity Pathway Analysis implicated numerous pathways associated with the progression to persistent pain. We took particular interest in the synaptogenesis signaling pathway as a key player at this 4-day time point, suggesting a role for synaptic plasticity in the transition from acute to chronic pain after SCI. In addition to its relevance within our model, many of the transcripts involved in synaptic plasticity overlapped with several other pathways (Figure 7C) and had RPKM values that could be validated by qPCR (Figure 7D). Synaptogenesis is typically associated with developmental processes, including axon guidance and synapse formation (Klein, 2004; Price and Inyang, 2015). However, activation of various signaling pathways involved in synaptogenesis may also contribute to pain; for example persistent pain is supported via changes in synaptic signaling, neuronal plasticity, and long term potentiation, and may form memory-like networks for painful signals that allow persistent pain to occur long after the initial injury (Kobayashi et al., 2007; Khangura et al., 2019).

Included in the many of the genes of interest within the synaptogenesis pathway, Ephb2 was significantly down-regulated post-injury (Table 1). The gene transcript is part of the Ephrin tyrosine kinase receptor protein family that is expressed in laminae I-III of the spinal dorsal horn on small and medium-sized DRG neurons (presumably nociceptors; Bundesen et al., 2003). Ephb receptors regulate synaptic activity in the spinal cord and contribute to persistent pain associated with NMDA activity (Khangura et al., 2019). Numerous receptor tyrosine kinases localize to synapses and contribute to synaptogenesis in addition to EphB receptors, including Trk receptors (Biederer and Stagi, 2008). Ntrk2, the receptor for BDNF, was highly expressed in this cell population; commensurately, BDNF transcript levels were significantly upregulated in the injured population (Table 1). Camk2g transcript levels were significantly increased in the SCI population of cutaneous nociceptors as well, and recent work has shown phosphorylation of Camk2g induces Bdnf mRNA transcription (Melemedjian et al., 2013; Yan et al., 2016). This parallels increasing evidence that neuronal activity (such as increased activity or hyperexcitability) activates alternative neuronal circuits through activity-regulated genes, such as BDNF (Lu et al., 2009; Melemedjian et al., 2013). Many of the genes significantly altered in the synaptogenesis pathway may function together to generate neuropathic pain (Supplementary Figure S7).

Regulatory interrelationships predicted by the IPA program were also examined (Figure 7E). Many of these networks are related to inflammatory signaling mechanisms, suggesting a link between pro-inflammatory signaling and synaptic transmission (Alexander and Popovich, 2009; Walters, 2014; Medelin et al., 2018). Previous work has associated inflammatory mechanisms to diseases of the CNS, including multiple sclerosis, Alzheimer’s disease, and Parkinson’s disease (Medelin et al., 2018). The mechanisms through which inflammation prompts changes in synaptic transmission are not fully understood, but several of the classic pro-inflammatory cytokines (including TNF-α, IL-6, IL-1β) have been shown to contribute to a decrease in hippocampal neurogenesis, and could be playing a similar role within the spinal cord following injury (Kohman and Rhodes, 2013). Overall, pro-inflammatory conditions have been associated with increases in post-synaptic NMDA and AMPA receptors, and inhibition of GABAergic receptors (Medelin et al., 2018).