Single-cell transcriptomic analysis of somatosensory neurons uncovers temporal development of neuropathic pain

Peripheral nerve injury could lead to chronic neuropathic pain. Understanding transcriptional changes induced by nerve injury could provide fundamental insights into the complex pathogenesis of neuropathic pain. Gene expression profiles of dorsal root ganglia (DRG) in neuropathic pain condition have been studied. However, little is known about transcriptomic changes in individual DRG neurons after peripheral nerve injury. Here we performed single-cell RNA sequencing on dissociated mouse DRG cells after spared nerve injury (SNI). In addition to DRG neuron types that are found under physiological conditions, we identified three SNI-induced neuronal clusters (SNIICs) characterized by the expression of Atf3/Gfra3/Gal (SNIIC1), Atf3/Mrgprd (SNIIC2) and Atf3/S100b/Gal (SNIIC3). These SNIICs originated from Cldn9+/Gal+, Mrgprd+ and Trappc3l+ DRG neurons, respectively. Interestingly, SNIIC2 switched to SNIIC1 by increasing Gal and reducing Mrgprd expression 2 days after nerve injury. Inferring the gene regulatory networks after nerve injury, we revealed that activated transcription factors Atf3 and Egr1 in SNIICs could enhance Gal expression while activated Cpeb1 in SNIIC2 might suppress Mrgprd expression within 2 days after SNI. Furthermore, we mined the transcriptomic changes in the development of neuropathic pain to identify potential analgesic targets. We revealed that cardiotrophin-like cytokine factor 1, which activates astrocytes in the dorsal horn of spinal cord, was upregulated in SNIIC1 neurons and contributed to SNI-induced mechanical allodynia. Therefore, our results provide a new landscape to understand the dynamic course of neuron type changes and their underlying molecular mechanisms during the development of neuropathic pain.


INTRODUCTION
Nerve injury could cause neuropathic pain, which affects 7%-10% of the general population. 1 Neuropathic pain severely debilitates patients' life quality. Peripheral nerve injury could induce a series of cellular and molecular responses including the gene regulation, axonal degeneration, activation of satellite glial cells and Schwann cells, infiltration of inflammatory cells, and the increase in vascular permeability in the injured dorsal root ganglion (DRG). 2,3 The dynamic regulation of the gene expression profiles in the injured DRG has been extensively studied, and could be an essential trigger in neuropathic pain. [4][5][6][7] For example, the expression level of voltage-dependent calcium channel subunit alpha-2/delta-1 (α2δ-1) encoded by Cacna2d1 was increased in DRG neurons after peripheral nerve injury. 4-7 α2δ-1 played a role in neuropathic pain by enhancing the insertion of NMDAR into the cell membrane. 8 Besides the neurons, non-neuronal cells also play a vital role in regulating neuropathic pain. 9 Satellite glial cells (SGCs) are peripheral glial cells embracing the neuron soma in DRGs. Peripheral nerve injury could induce a variety of changes in SGCs, including cell proliferation, the gap junctions formation, the expression of Kir4.1, the ATP signaling and the activation of cytokines. 10 Upregulated gap junctions in SGCs contributes to the pain transduction by enhancing activation of DRG neurons after peripheral nerve injury. 11 Taken together, the mechanisms underlying neuropathic pain are intricate. As a result, neuropathic pain is still hard to cure. Clinically, the drugs for treating neuropathic pain are usually addictive or less potent. 12 In addition, some may evoke considerable side effects. 12 Therefore, it is necessary to further explore the underlying mechanism for the occurrence and development of neuropathic pain at molecular and cellular levels.
Single-cell RNA sequencing (scRNA-seq) is a powerful technology to examine the transcriptome of a single cell. 13 Recently, scRNA-seq was applied to determine the heterogeneity of somatosensory neurons and identify the DRG neuron types. 14,15 Using Smart-seq2, we have identified at least 10 types and 14 subtypes of somatosensory neurons in the mouse DRG, and provide a sketch for sensory functions for each neuron type. 14 Besides the cell populations, gene regulatory networks can be inferred through cis-regulatory analysis at the single-cell resolution. 16 The algorithm for analyzing scRNA-seq data reveals the kinetics of single-cell gene expression in a complex biological progression. [17][18][19][20] Thus, we speculate that scRNA-seq could provide us with a deep insight into the dynamic progress on neuropathic pain. Moreover, it would be interesting to know whether and how the dynamic DRG neuron type change happens in response to peripheral nerve injury. Quantitative analysis on molecular and cellular alterations in individual cells and among different cell types may help search for potential targets for developing pain therapies.
In the present study, we analyzed the single-cell transcriptomes of mouse DRG by 10× Genomics scRNA-seq at different time points after the spared nerve injury (SNI), a widely used neuropathic pain model. 21 Distinct non-neuron cell types and DRG neuron types, including three new SNI-induced neuron clusters (SNIIC1-3), were identified with distinct transcriptome profiles. Three SNIICs were all closely related to nerve injury by the expression of Atf3, a well-known marker for neuron injury. We identified that the origins of three SNIICs were Cldn9 positive (Cldn9 + ), Mrgprd + and Trappc3l + DRG neurons, and found that the SNI-upregulated genes were mainly expressed in these SNIICs. We revealed the transient appearance of Atf3-high-expression (Atf3 high )/Mrgprd + neurons 24 h after SNI, and then they switched to Atf3 high /Gfra3 + /Gal + neurons within 2 days after SNI. The gene regulatory network closely related to nerve injury was also identified. The regulators of Gal and Mrgprd were evaluated. Finally, the pseudo-time analysis was performed to identify the differentially expressed genes (DEGs) in neuron type switches. The upregulated genes functioned in the synapse organization, ion transport, and inflammatory response. Cardiotrophin-like cytokine factor 1 (CLC) mRNA, one of the upregulated genes in SNIIC1, was required for the development of neuropathic pain. Our data provide a comprehensive landscape to reveal the dynamic changes of individual DRG neurons during the development of neuropathic pain and discover a new algogenic cytokine.

Cell heterogeneity in DRG
In the present work, we intended to determine the single-cell transcriptomes of DRG neurons in adult mice under neuropathic pain conditions. Rodents were subjected to the unilateral SNI, which could induce a conspicuous and long-lasting mechanical allodynia. 21 We performed 10× Genomics scRNA-seq on the cells dissociated from lumbar (L)4 and L5 DRGs ipsilateral to the SNI at different time points: 6 h, 24 h, 2 days (2d), 7d, and 14d after SNI (Fig. 1a). As a result, 36,810 cells, including 6935 neurons, were retained after removing the empty droplets, outlier, cell debris, and inferential doublets (Supplementary information, Table S1). The details were described in the "Material and Methods". Through principal component analysis, we performed an original Louvain algorithm to classify DRG cells and t-distributed stochastic neighbor embedding (t-SNE) to visualize the data by Seurat package. 17,20,22 DRG cells were classified into nine major types with distinct molecular markers (Fig. 1b, c; Supplementary information, Fig. S1b, c). Identified cells included primary sensory neurons, SGCs, Schwann cells, vascular endothelial cells (VEC), vascular smooth muscle cells (VSMC), vascular endothelial cells, capillary (VECC), fibroblasts, immune cells and red blood cells (RBC) (Fig. 1b). Averagely, 6407 genes and 53,904 unique molecular identifiers (UMIs) can be detected in a single DRG neuron and 1605 genes and 3723 UMIs in a single non-neuron cell using saturated sequencing (Supplementary information, Table S1 and Fig. S1a, b). SGCs surrounded the soma of DRG neurons and expressed Apoe and Fabp7 as shown by in situ hybridization (ISH) ( Fig. 1c; Supplementary information, Figs. S1c, S2a) and in previous report. 14 Schwann cells formed the myelin sheath in the peripheral nerves and were characterized by the expression of Mpz and Pllp ( Fig. 1c; Supplementary information, Figs. S1c, S2a). We also observed a cluster of fibroblast cells with the expression of Dcn and Apod ( Fig. 1c; Supplementary information, Figs. S1c, S2a). Here, we observed three clusters of vascular cells in the DRGs, including the VEC marked by Cldn5 and Ly6c1, VSMC by Tpm2 and Acta2 and VECC by S100a8 ( Fig.1c; Supplementary information, Fig. S1c). Besides, we identified immune cells expressing Cd74 and RBC with Hba-a1 ( Fig. 1c; Supplementary information, Fig. S1c). The classification of non-neuron cells after SNI was consistent with that under normal conditions (Supplementary information, Fig. S1d), and biological repetition yielded the same results (Supplementary information, Fig. S2b-d). Thus, we provided a complete gene expression profile for each DRG cell type.

DRG neuron types after nerve injury
To explore the heterogeneity of DRG neurons under normal and neuropathic pain conditions, we reclassified DRG neurons based on their transcriptional characteristics. Neurons in two control samples were evenly distributed on the t-SNE plot (Supplementary information, Fig. S1e). Nineteen neuron clusters were unbiasedly identified with distinct molecular markers (Fig. 1d, e). The representative genes were used to annotate the clusters ( Fig. 1e; Supplementary information, Table S2). There were 6 clusters of peptidergic neurons (clusters 1-6), 5 clusters of non-peptidergic neurons (clusters 7-11), 5 clusters of myelinated neurons types (clusters 12-16) and 3 SNIICs with high expression levels of Atf3 (Figs. 1d and 2b). To provide an adequate resource, 10× Genomics experiments for DRG neurons after SNI model at all time points were repeated, following integration analysis using Seurat. Consistent results with similar marker genes were identified in the repeated trials (Supplementary information, Fig. S2e). Due to Fig. 2 Integrated analysis of DRG neuron sequenced by 10× Genomics and Smart-seq2. a t-SNE plots of 3227 (479 for Smart-seq2 and 2748 for 10× Genomics) L4 and L5 DRG neurons sequenced by 10× Genomics and Smart-seq2 after alignment. Dots, individual cells. b Summary of neuron classification based on 10× Genomics and the previously reported classification. c ISH shows the Prokr2-and Smr2-expressing neurons in mouse L5 DRG (n = 3). Scale bar, 50 µm. d ISH shows Dcn expression in both DRG neurons (arrow) and fibroblast (arrowhead). Scale bar, 50 µm. e RNAscope result shows the Rxfp1-expressing DRG neuron (arrow). Scale bar, 50 µm. f Statistical results show the percentages of Dcn + and Rxfp1 + DRG neurons (n = 3). g RNAscope assay shows that the Zcchc12 (arrowhead) and Cldn9 (arrow) are not co-expressed in DRG neurons. Scale bar, 20 µm. h Dual fluorescent ISH shows that Dcn + neurons are a subset of tdTomato + neurons in Zcchc12-CreER::Ai9 mouse L5 DRG (arrow). Scale bar, 20 µm. i RNAscope assay shows that Rxfp1 + neurons are a subset of tdTomato + neurons in Zcchc12-CreER::Ai9 mouse L5 DRG (arrows). Scale bar, 20 µm. The data are shown as means ± SEM.
Article the increase of cell number, a small neuron subcluster of SNIIC1 marked by Fam19a4 and Atf3 was identified at SNI 24 h and SNI 2d (Supplementary information, Fig. S2e-g).
Previously, we separated DRG neurons into 10 types and 14 subtypes based on Smart-seq2 technology. 14 We also collected single-cell transcriptomes of DRG neurons on 14d after SNI by Smart-seq2 technology as reported. 23 We used Seurat to integrate the data from control and SNI 14d by 10× Genomics with the data by Smart-seq2 technology (Fig. 2a). As a result, we identified 18 clusters ( Fig. 2a; Supplementary information, Fig. S3a). The number of detected genes by Smart-seq2 was more than by 10× Genomics in all neuron types (Supplementary information, Fig. S3b). Then we matched 18 clusters with our previous report 14 (Fig. 2b). Most clusters were consistent with previous results (Fig. 2b). However, the previous C6 was characterized by the expression of Mrgprd in myelinated DRG neurons. C10 was not identified because 10× Genomics technology could not distinguish neuronal size. 14 Previously, C8 was subdivided into two clusters represented by Ntrk3 high and Ntrk1 high , respectively. In the present work, C8 could be subdivided into three clusters. Previous Ntrk1 high C8 neurons were classified into Prokr2 + (C8-2) and Smr2 + (C8-3) clusters (Fig. 2b). Prokr2 was expressed in~1.7% DRG neurons, while Smr2 in~1.6% DRG neurons (Fig. 2c).
To describe a complete atlas of DRG neuron types under both physiological and pathological conditions, we used the semisupervised analysis to modify the classification of our Smart-seq2 database. C6 neurons were separated from large neurons due to the shared marker genes of C5 neurons (Mrgprd) and large neurons (S100b and Nefh), while C2-2 neurons were sub-grouped within C2 due to the fact that they also expressed the marker genes of large neurons (Supplementary information, Fig. S4a). C10 neurons were separated from C1 according to the neuron size (Supplementary information, Fig. S4a, b). Finally, 16 DRG neuron types or subtypes were revealed by 10× Genomics and 19 types by Smart-seq2 in normal condition. And three SNI-induced clusters were revealed by 10× Genomics.
Emerging neuron types induced by nerve injury To investigate the development and progression of emerging neuron types induced by SNI, we split the t-SNE plot according to the time points. Neurons in clusters 1-16 were evenly distributed at all time points (Fig. 3a). Transcriptomic changes (q-value ≤ 0.01 and average logFC ≥ 1) in these 16 clusters were inconspicuous (Supplementary information, Table S3). Three emerging clusters appeared after SNI (clusters 17,18,19), namely SNIIC1-3 (Fig. 3a). Atf3 was extensively expressed in every cluster after SNI, but highly expressed in SNIIC1-3 (Supplementary information, Fig.  S6a). The abundance of SNIIC1 labeled by Atf3/Gfra3/Gal and SNIIC3 labeled by Atf3/S100b/Gal peaked at SNI 2d and maintained at SNI 14d (Fig. 3a), while SNIIC2 marked by Atf3/Mrgprd peaked at SNI 24 h and then disappeared (Fig. 3a). To study the duration of SNIIC1 and SNIIC3 appearance after SNI, we performed scRNA-seq of L4, L5 DRG cells at SNI 28d. Both SNIIC1 and SNIIC3 were maintained at SNI 28d (Supplementary information, Fig. S5a). However, the percentage of Atf3 + neurons were decreased at 28d after nerve injury. 27,28 Therefore, SNIICs might disappear in a long period after SNI. Hundreds of DEGs were identified in these three SNI-induced neuron clusters (myAUC > 0.7) (Fig. 3b). SNIICs was represented by Atf3, Sox11, Cacna2d1, Gadd45a, Cckbr, and Sprr1a (Fig. 3b), the well-known upregulated genes in DRGs after peripheral nerve injury. 4,7 Gene Ontology (GO) enrichment analysis showed that those representative genes in SNIIC1 and SNIIC3 were involved in inflammatory response, receptor signaling, and apoptotic process (Supplementary information, Fig. S5b). Then, we evaluated the expression of Atf3 in DRG neurons after SNI. ISH showed the presence of Atf3 in~50% neurons at SNI 2d, and~30% at SNI 14d (Supplementary information, Fig. S6b). Dual fluorescent ISH showed that Atf3 was co-expressed with Gfra3 iñ 5% DRG neurons at SNI 24 h and~12% at SNI 2d ( Fig. 3c, d). Atf3 was co-expressed with S100b in~15% DRG neurons at SNI 24 h and~17% at SNI 2d (Fig. 3c, d). The co-expression of ATF3 with Gfra3 or S100b at SNI 7d and 14d were also confirmed by fluorescent ISH combined with immunostaining (Supplementary information, Fig. S6g). Moreover, SNIIC1 and SNIIC3 were also marked by Gal (Fig. 1e), a neuropeptide that contributes to axonal regeneration and modulation of neuropathic pain. 29 Nppb + / Il31ra + neurons also expressed Atf3 at a moderate expression level, but not Gal (Supplementary information, Fig. S6c-f). Atf3 was co-expressed with Mrgprd in~6% DRG neurons at SNI 24 h but less than 1% at SNI 2d (Fig. 3c, d). Thus, our data showed that SNIIC2 neurons transiently existed. Noticeably, SNI could result in a decrease in the number of Mrgprd + DRG neurons and an increase in the number of Gfra3 + DRG neurons (Fig. 3d).
Lastly, we examined the origin of SNIIC1. SNIIC1 expressed Gfra3, Gal and Atf3. Given that Gfra3 was highly expressed in C1 and C2 neurons, we performed the lineage tracing using Zcchc12-CreER and Sst-Cre transgenic mice. The results revealed that only a few of Sst + neurons expressed Gal after SNI (Supplementary information, Fig.  S8d), and only a few of Zcchc12 + neurons expressed Atf3 after SNI (Supplementary information, Fig. S8e). Thus, neither Zcchc12 + neurons (C1-2) nor Sst + neurons (C2) could be the major origin of SNIIC1. The abovementioned data showed a decrease of Cldn9 + and Mrgprd + neurons induced by SNI. The results indicated that Cldn9 + and Mrgprd + neurons (C1-1 and C5) could be the major origin of SNIIC1. To evaluate the fate of transiently appeared SNIIC2 after SNI, we utilized a recombinase line, Mrgprd-CreER, for the lineage tracing of Mrgprd + neurons. We crossed Mrgprd-CreER mice with reporter line Ai9 to label the Mrgprd + neurons before SNI. We performed ISH and immunohistochemistry to evaluate the expression of Gal, Atf3 and Mrgprd in the DRGs of Mrgprd-CreER::Ai9 mice after SNI (Fig. 4e). Gal was not expressed in tdTomato-labeling neurons (Fig. 4e), and Mrgprd was expressed almost in all tdTomatolabeling neurons under normal condition (Fig. 4e). However,~18% of tdTomato-labeled neurons expressed Gal but not Mrgprd at SNI 2d (Fig. 4e, f), while~26% of tdTomato-labeled neurons expressed Atf3 but not Mrgprd at SNI 2d (Supplementary information, Fig. S8f). Thus, SNIIC2 might switch to SNIIC1 since SNI 2d. Three emerging neuron clusters appeared at different time points after SNI. a t-SNE plot shows somatosensory neuron clusters in control and at different time points after SNI. Dots, individual cells; colors, neuron clusters. b A heatmap shows the DEGs in the three emerging neuron clusters. The color represents the scaled transcript counts. c, d Dual fluorescent ISH confirms the co-expression of Atf3 with Mrgprd, Gfra3, or S100b at different time points after SNI (n = 3). Scale bar, 20 µm. The data are shown as means ± SEM. *P < 0.05, ***P < 0.001, n.s. no significance vs the control group.
Critical regulons involved in the dynamic process of neuron type switches To identify the critical regulators and gene regulatory network during the neuron type switches, we adopted SCENIC pipeline to infer the regulon activity and gene regulatory network. 16 The regulon activity matrices were displayed (Supplementary information, Fig. S9a and Table S4). Atf3, Jun and Sox11, wellknown transcription factors (TFs) expressed in nerve-injured neurons, 5,28-30 were also identified. Since the regulons of Mrgprd were missed using the default parameters, Normalized Enrichment Score (NES) threshold was reduced from 3 to 2. The t-SNE plot was generated based on the binary regulon activity matrix (Fig. 5a). The distinct segmentation of each neuron type indicated the different patterns of transcriptional activation. SCENIC also identified the predicted TFs regulating Gal and Mrgprd ( Fig. 5b; Supplementary information, Table S4), and the enriched DNAbinding motif (Fig. 5c). Those TFs were expressed and activated in the neurons of three SNIICs ( Fig. 5c; Supplementary information Fig. S9a, b). However, Atf3, Jun and Egr1 had broad expression patterns. Egr1 was also positive in a part of S100b + , Cldn9 + and Zcchc12 + neurons, while Atf3 and Jun were also activated in Cldn9 + or Zcchc12 + neurons (Fig. 5c). The activation of Crem was restricted in SNIIC1-3 neurons, while Bach1 only in SNIIC1 (Fig. 5c). Cpeb1 was activated in but not limited to Mrgprd + neurons, particularly in SNIIC2 (Fig. 5c).
To validate the gene regulatory network predicted by SCENIC, FACS was performed to collect tdTomato-labeled DRG neurons from Gal-Cre::Ai9 and Mrgprd-CreER::Ai9 mice after SNI. The tdTomato-labeled neurons were sorted and cultured for experiments. The siRNAs of TFs were transfected into the cultured DRG neurons. The expression of target genes was detected by RT-PCR (Fig. 5d). The results showed that knockdown of Atf3 or Egr1 downregulated Gal expression (Fig. 5e), while knockdown of Cpeb1 upregulated the Mrgprd expression (Fig. 5f). Thus, these identified regulons could play a role in regulating the expression of Gal and Mrgprd in DRG neurons after SNI, and be involved in the formation of SNIICs.
Nerve injury-induced dynamic gene regulation based on singlecell trajectories To identify the dynamic gene regulation in three SNI-induced neuron clusters, we performed the pseudo-time analysis to reconstruct the single-cell trajectories. 18,19,30 The reconstructed tree reflected the neuron fate and the change of gene expression during SNI (Fig. 6a, b; Supplementary information, Fig. S10a-d).
We identified a total of 6469 DEGs, including 1040 shared genes among the three processes of neuron-cluster switches induced by  Table S5, P-value < 10 −6 ). More DEGs (4887) were identified among the switch from Mrgprd + to Mrgprd + /Atf3 high , then to Gfra3 + /Atf3 high neurons, suggesting that injured Mrgprd + neurons have higher reprogramming capabilities (Supplementary information, Fig. S10b, e and Table S5). The DEGs in the switching process of Cldn9 + neurons were classified into 5 gene modules according to the expressing patterns across pseudo-time (Fig. 6b). To investigate the function of the regulated genes, we applied GO enrichment analysis to the DEGs. The DEGs in three switch processes have similar functions, such as synapse organization, ion transport, inflammatory response, and axon development ( Fig. 6c; Supplementary  information, Fig. S10f). Notably, the genes related to sensory perception of pain were decreased (Fig. 6b, c). These enriched genes in sensory perception of pain mostly were related to pain transduction, including Trpv1, Scn9a, Scn10a, and Scn11a. 31,32 Comparatively, the expression of genes associated with inflammatory responses was increased (Fig. 6b, c). The genes enriched in inflammatory response related to pathological pain included Csf1 and Fam19a4. 33,34 The regulated genes included the TFs, cytokines and related molecules (including receptors, ligands, suppressors and induced or regulated proteins), Gprotein-coupled receptors (GPCRs), ion channels, and growth factors and their receptors ( Fig. 6c; Supplementary information, Table S5). Real-time PCR proved that the expression of Crem and Crlf1 was increased after SNI, while Htr3a and Trpc4 were downregulated (Supplementary information, Fig. S11a). Atf3, Sox11, Sprr1a, and Gap43 were functionally associated with the axon regeneration. 5 CSF1 was induced in SNIIC1-3 neurons (Fig. 6c; Supplementary information, Fig. S11b) and played a role in neuropathic pain by activating the spinal microglia. 34 Blockade of cholecystokinin-B receptor encoded by Cckbr partially relieved neuropathic pain. 35 P2rx3 was upregulated in rat DRG neurons after peripheral nerve injury. 36 Downregulation of P2rx3 in DRG neurons could reduce neuropathic pain behavior. 37 Thus, the regulated genes in SNI participate in mediating nerve regeneration or neuropathic pain.

Inflammatory responses induced by peripheral nerve injury
Peripheral nerve injury induced the expression change of cytokines and related molecules, including receptors, ligands, suppressors, and induced or regulated proteins ( Fig. 6c; Supplementary information, Table S5). Previous reports showed that CSF1 and PAPs contributed to neuropathic pain by activating the spinal microglia. 34,38 Inflammatory responses contributed to neuropathic pain. 39 Then, we evaluated the results by RNAscope Fig. 6 Pseudo-time analysis reveals the gene expression changes and functional adaptation during Cldn9 + neuron type-switch. a Singlecell trajectory shows the fates of Cldn9 + neurons after SNI. Two branches were identified, one led to Atf3 + /Gfra3 + /Gal + neurons from Cldn9 + neurons (Fate1, F1), the other remained as it was (F2). The color indicates the time points after SNI (top) or represents neuron clusters (bottom). b Gene signatures during Cldn9 + neuron switch. The DEGs (rows) are shown using heatmap (left) with neurons (columns) in pseudotime from Root to F1. Gene-expression trends in each group (middle). GO terms associated with DEGs in the four kinetic clusters (right). c Heatmaps show the DEGs identified by BEAM (rows), and cells (columns) are shown in pseudo-time from Root to F1.
assay. Our data confirmed that Clcf1, Il17ra, Cd109, and Il34 were upregulated in DRG neurons after SNI ( Fig. 7a; Supplementary  information, Fig. S11c). Moreover, those genes were expressed mainly in one or more SNIICs (Supplementary information, Fig.  S11b). Clcf1 encodes cardiotrophin-like cytokine factor 1 (CLC) protein while Crlf1 encodes cytokine receptor-like factor 1 (CLF) protein. CLC belongs to interleukin (IL)-6 family cytokines, whose receptor complex contains gp130, a signaling receptor subunit. 40,41 CLC forms a stable complex with CLF, an orphan cytokine receptor, to mediate signal transduction. 42,43 The upregulation of Clcf1 induced by SNI has not been shown by previous reports. Thus, the function of CLC/CLF complex in pain modulation has not been studied. Real-time PCR showed that Clcf1 was upregulated in the DRG but not in the spinal cord after SNI (Supplementary information, Fig. S12a). RNAscope assay showed that the percentage of Clcf1 + DRG neurons was increased from~3.8% under normal condition to~14.8% at SNI 2d,~14.5% at SNI 7d and~5.6% at SNI 14d (Fig. 7a). Both Clcf1 + and Crlf1 + were present in SNIIC1 neurons (Supplementary information, Fig.  S11b), and~80% Clcf1 + DRG neurons contained Crlf1 at SNI 7d (Fig. 7b). ELISA analysis revealed the secretion of CLC protein in cultured DRG neurons. We found that the level of CLC secreted from the injured DRG neurons was higher than that from the control DRG neurons (Fig. 7c). To evaluate the function of CLC in the pain regulation, CLC/CLF complex was intrathecally injected in the wild-type mice. Then, the von Frey test was performed to examine the mechanical threshold in response to the pronociceptive mechanical stimulus. Intrathecal CLC/CLF complex evoked Fig. 7 Clcf1 is involved in the temporal development of neuropathic pain. a RNAscope results show that the percentage of Clcf1 + neurons was increased after SNI (n = 3). Scale bar, 20 μm. b RNAscope results show the co-expression of Clcf1 and Crlf1 after SNI (n = 3). Scale bar, 20 μm. c ELISA results showing CLC levels in cultured DRG neurons from control and SNI mice (n = 11 for each group). d Intrathecal injection (i.t.) of CLC/CLF dose-dependently decreased the mechanical thresholds in C57 mice (n = 14 for PBS group; n = 12, 11, 12 for CLC/CLF group of 6, 30, or 150 ng injection, respectively). e Immunostaining results show that the astrocytes are activated 4 h after CLC/CLF complex injection (30 ng per mouse). Scale bar, 50 μm. The enlarged ROI and outline are shown on the right. Scale bar, 20 μm. f, g The signal of activated astrocytes and their distribution in spinal cord area were both increased after CLC/CLF injection. The data are shown as means ± SEM (n = 3). h A flowchart shows that a basal threshold was tested before the injection, and the mice were injected with the siRNA for two consecutive days (4 μg each time). The von Frey test was performed on SNI 4d, and the knockdown efficiency of siClcf1 was detected from the mouse L4, L5 DRGs. i qPCR results show that the expression of Clcf1 was reduced in the DRGs after siClcf1 injection. j Knockdown of Clcf1 partially alleviated the mechanical allodynia induced by SNI (n = 18 for control and n = 18 for siClcf1). The data are shown as means ± SEM. *P < 0.05, **P < 0.01, ***P < 0.001, n.s. no significance vs the control group.
the mechanical allodynia in a dose-dependent manner (Fig. 7d). Previous reports demonstrated that the receptor of CLC/CLF complex, ciliary neurotrophic factor receptor (CNTFR), was expressed in astrocytes in CNS. 44 The spinal astrocytes also expressed Cntfr 45 (Supplementary information, Fig. S12b). Then, we examined the effect of CLC/CLF complex on the activation of spinal astrocytes. Intrathecal injection with CLC/CLF complex could activate the astrocytes but not microglia in the spinal dorsal horn (Fig. 7e-g, Supplementary information, Fig. S12c). Next, we evaluated the requirement of CLC/CLF complex in the regulation of neuropathic pain (Fig. 7h-j). Intrathecal injection of siClcf1 could decrease the expression of Clcf1 in the DRG but not the spinal cord (Supplementary information, Fig. S12d). Inhibition of Clcf1 expression in DRG neurons by siRNA reduced the mechanical allodynia at SNI 4d (Fig. 7j). Thus, CLC contributed to pain hypersensitivity induced by peripheral nerve injury.

DISCUSSION
Somatosensory neurons play essential roles in sensory transduction in both physiological and pathological conditions. Previously, DRG neuron types and subtypes were identified under physiological condition by low-coverage and high-coverage scRNA-seq based on Smart-seq2 technology, respectively. 14,15 Recently, 10× Genomics scRNA-seq has been applied to define the transcriptional regulations of neuronal specification. 46 However, distinct changes and regulatory mechanisms in different clusters of DRG neurons after peripheral nerve injury are not well understood. Here, we performed 10× Genomics scRNA-seq on DRG cells after SNI. Firstly, we identified the emerging clusters closely related to SNI. Then, we demonstrated the neuron type switches and identified the different origins of SNIICs. Furthermore, we reconstructed the gene regulatory networks related to neuron type switches. Finally, we examined the expression changes of functional genes and revealed that cytokine CLC is a novel potential target for the development of analgesic therapies.
Classification of DRG neurons with scRNA-seq The heterogeneity of neuron types include molecular, morphological, connectional, and functional properties. 47 For DRG neurons, the cell sizes are crucial because small DRG neurons give rise to unmyelinated afferent fibers which project to the lamina I-II of spinal dorsal horn, while the myelinated afferent fibers of large DRG neurons project mainly to the lamina III-V of spinal cord. 48,49 The scRNA-seq is a powerful method for examining the gene profiles within individual cells based on cell state and activity. 13 The number of average detectable genes by Smart-seq2 technology was 10,950 genes per neuron in our previous study, while by using 10× Genomics technology 6407 genes per neuron was detected in the present study. The major types of mouse DRG neuron and their marker genes identified by these two scRNA-seq approaches were largely consistent, suggesting that they are both reliable in cell typing. We also noticed that both Smart-seq2 and 10× Genomics displayed their technical advantages and disadvantages that may affect the data interpretation.
Although 10× Genomics scRNA-Seq could isolate and process a large number of cells, its technical limitation is the lack of the information about cell sizes. Using Smart-seq2 which processed the cells of all different diameters manually collected from the dissociated DRGs, we have identified at least 10 types and 14 subtypes of mouse DRG neurons, including the C10 large neurons marked by Gal which also serves as a marker gene for one type of small neurons (C1). 14 However, such neuron type can not be identified by 10× Genomics scRNA-seq which is not capable of providing the information about the sizes of processed neurons. According to the instructions of 10× Genomics, the microfluidic channels are approximately 100 μm in diameter, and beads are 50 μm in diameter. Thus, large DRG neurons with 40-50 μm in diameter may not be easily analyzed by 10× Genomics scRNA-seq, and cells larger than 50 μm in diameter can not pass the microfluidic device for automatic collection of cells for the analysis.
Therefore, for the classification of neurons with different diameters and numbers in the tissue, both technical advantages and disadvantages of Smart-seq2 and 10× Genomics scRNA-seq should be carefully considered and evaluated for the data interpretation. Based on the integrated analysis of the results from Smart-seq2 and 10× Genomics scRNA-seq, we propose that the mouse DRG neurons under normal circumstance could be classified into 19 types/subtypes. Our current data provide the insight into new DRG neuron subtypes. However, more DRG neuron types or subtypes may be identified from the study of connectional and functional properties, as well as from a deeper sampling and sequencing.
DRG neuron types after nerve injury Cumulative evidence shows that peripheral nerve injury could induce the changes in the gene expression profiles of DRG neurons. Such changes contribute to the pathophysiological and progression of neuropathic pain. 4,5,50 Here, we found that three new neuron types appeared after SNI, namely SNIIC1 (Atf3/Gfra3/ Gal), SNIIC2 (Atf3/Mrgprd), and SNIIC3 (Atf3/S100b/Gal). SNIIC2 transiently appeared 24 h after SNI, then switched to SNIIC1 within 2 days after nerve injury. Interestingly, Atf3, Gal and other nerve injury-regulated genes were mainly expressed in these emerging neuron types. 27,51 Analogous injured-related clusters were also identified in trigeminal neurons after partial infraorbital trigeminal nerve transection 52 and lumbar DRG neurons after spinal nerve transection, sciatic nerve transection + ligation (ScNT), or sciatic nerve crush. 28 By single-nucleus RNA sequencing (snRNA-seq), Renthal et al. have recently revealed the transcriptional program that responds to axonal regeneration in mouse models of peripheral axotomy. 28 However, snRNA-seq only captured RNA in nucleus which lead to lower sequencing depth (1284 unique genes per nucleus). In the present study, we described a complete atlas of DRG neuron types in both physiological and pathological conditions by Smart-seq2 (10950 genes per neuron) and 10× Genomics with high sequencing depth (6407 genes per neuron). Besides, we investigated the change of functional genes after SNI, explored the regulatory mechanisms underlying neuropathic pain and identify a novel target (Clcf1) for the developing potential analgesic therapies.
Gene expression profiles in various cells are regulated by TFs. 53,54 Therefore, the nerve injury-induced changes in gene expression seemed to mainly occur in distinct subgroups of DRG neurons, suggesting that such changes could be due to the injury-induced alterations of TF expression in the subgroups of DRG neurons.
Altered gene expression profiles in the individual neurons lead to the switch from one neuronal type or subtype to another one. By tracing the trajectories of neurons, we found that the C5 neurons under physiological condition could switch to SNIIC2 shortly after SNI. And then together with C1-1 they became SNIIC1, while C8 switched to SNIIC3. Moreover, the upregulated TFs, such as Atf3 and Egr1, could increase the expression of Gal, while the high level of Cpeb1 could reduce the expression of Mrgprd. Therefore, the fates of SNIICs could be further determined by the changes in the expression of TFs.
The induced DRG neurons in SNI may manifest the functions different from that of original neurons. SNIIC2 neurons were originated from Mrgprd + C5 itch-sensitive MHNs and were switched to SNIIC1 by increasing Gal and reducing Mrgprd expression 2 days after SNI. Thus, a certain subpopulation of Mrgprd + itch-sensitive mechanoheat nociceptor (MHN) switched to SNIIC1 neurons and no longer expressed MRGPRD to mediate the β-alanine-induced itch sensation. SNIIC1 neurons were originated from Cldn9 + /Gal + C1-1 subtype of C1 MHN, and C1-1 subtype showed the potential to respond to mechanical nociceptive stimulations. 14 After nerve injury, SNIIC1 neurons originated from C1-1 neurons and SNIIC2 may obtain new functions such as producing inflammatory cytokines including CSF1 and CLC. SNIIC3 neurons which originated from Trappc3l + C8 mechanoreceptor 14 also express inflammatory CSF1. The present study shows that SNIICs could be induced after nerve injury and play important roles in the development of neuropathic pain.
Gene regulations in the development of neuropathic pain Our previous study showed that the gene expression profile of DRG was altered 2 days after peripheral nerve injury, 4 suggesting that such changes could occur shortly after nerve injury. In the present study, we found that the gene expression profiles at 6 h after nerve injury could be changed markedly in specific subpopulations of DRG neurons, which led to the formation of emerging new neuron types. Three emerging neuron types distinctly appeared, and SNIIC1 and SNIIC3 were stabilized 2 days after nerve injury. Therefore, the molecular basis for developing persistent pain was initiated at single neuron level as early as 6 h after nerve injury.
Hyperexcitability of DRG neurons is a critical cause of neuropathic pain. Decades of work have been dedicated to the structure and function of voltage-gated sodium channels in pain transmission. 12,55 Notably, gain-of-function mutations in Na v 1.7, Na v 1.8, and Na v 1.9 (coded by Scn9a, Scn10a and Scn11a) can increase spontaneous firing and contribute to painful peripheral neuropathy. [56][57][58] Moreover, voltage-gated potassium channels also play a crucial role in pain processing. 59 For example, dysfunction of K v 2.1 and K v 2.2 (coded by Kcnb1 and Kcnb2) enhance DRG neuron excitability. 60 Here, we found the downregulation of Scn9a, Scn10a, Scn11a, Kcnb1, and Kcnb2 as in previous reports. 4,50,60 We also newly found the downregulation of Kcnv1 and Kcng2 and the upregulation of Kcnk12, Kcnk13 and Kcnk16. Upregulated CCK-B receptor and P2X3 were involved in the formation of neuropathic pain. 35,61 Downregulation of opioid receptors could reduce the neuronal sensitivity to opioid analgesics. 62 Thus, regulated genes could contribute to neuropathic pain. Our data provide more potential drug targets to develop analgesic strategies.
It has been well-accepted that the mechanisms of neuropathic pain rely on the changes in both neurons and glial cells. The spinal microglia mediate the immune responses during the development of neuropathic pain. 63 After peripheral nerve injury, the spinal microglia are quickly activated by DRG neuron-derived proinflammatory factors, such as CSF1 and chemokine (C-C motif) ligand 2 (Ccl2). 34,64,65 Activated microglia could aggregate around the central terminals of injured DRG neurons, secreting cytokines and chemokines to induce the central sensitization. 63,66 Chemokine receptors, such as C-C chemokine receptor type 2 (CCR2) and CX3C chemokine receptor 1, participate in this process. 67 At the later stage, nerve injury-induced PAP-I acted as a central proinflammatory factor and maintained SNI-induced tactile allodynia via activating microglial CCR2, mediating the neuron-microglial interaction in the spinal cord. 38 In this study, several inflammation-related genes, such as Clcf1, Il17ra, Cd109, and Il34, were upregulated in DRG neurons after peripheral nerve injury. The blockade of IL-17R could suppress the hyperexcitability of DRG neurons and neuropathic pain behavior induced by paclitaxel treatment. 68 IL34 functioned as the other ligand for CSF1 receptor. 69 We found that CLC expressed in DRG neurons after nerve injury could induce mechanical allodynia, a typical neuropathic pain behavior. CLC could activate the spinal astrocytes in the dorsal spinal cord. CTNFR was express in the astrocytes in brain. 44 Therefore, peripheral nerve injury could increase the expression of inflammatory factors in certain types of DRG neurons, and secretion of those factors in the dorsal spinal cord may activate either microglia or astrocytes, mediating the neuron-glial interactions to participate in the development of nerve injury-induced tactile allodynia.
The present study reveals that the single-cell transcriptomic alterations of somatosensory neurons after peripheral nerve injury regulate the gene expression patterns and generate emerging neuron types. The Clcf1 and genes expressed in these emerging neuron types may be involved in the temporal development of neuropathic pain. Thus, these transcriptome data of changes in DRG neuron types after peripheral nerve injury can be a resource for studying the mechanisms of neuropathic pain and developing pain therapy.

MATERIAL AND METHODS
Animals and genotyping C57BL/6J mice were used in the experiments, according to the guidelines of the Committee of Use of Laboratory Animals and Common Facilities, Institute of Neuroscience, Chinese Academy of Sciences, Shanghai, China. The male mice were used at 7-8 weeks of age.

Surgery and tamoxifen induction
The mice were anesthetized with isoflurane during the surgery. SNI model was made according to the previous protocol. 21 Briefly, the tibial and common peroneal nerve branches were cut off, leaving the remaining sural nerve intact. After SNI surgery, the C57BL/6J mice were allowed to survive for 6 h, 24 h, 2d, 7d, and 14d. For the initiation of tdTomato expression, tamoxifen (Sigma T5648) was dissolved into corn oil with 20 mg/mL. The mice were injected intraperitoneally with 100 μg per 1 g of body weight for 5 consecutive days Detecting empty droplets and quality control Empty droplets often contain RNA from the ambient solution, resulting in non-zero counts after debarcoding. The package DropletUtils 74,75 (v.1.0.3) implemented in R was applied to distinguish the empty drops and cells. Drops which had more than or equal to 1000 UMI counts were assumed to correspond to empty droplets. Quality control metrics were calculated using scater 76 (v.1.8.4) and the median absolute deviation (MAD) to defined outliers. We removed cells whose log-library sizes, the logtransformed number of expressed genes, were more than 3 MAD below or above the median value of each cluster, or the proportions of mitochondrial genes (mitochondrial genes were used as a proxy for cell damage) were more than 4 MAD above the median value of each cluster. Finally, we removed the clusters with ambiguous markers (e.g., cell debris).

Identification of cell clusters
The clean data were further processed with Seurat 17,22 (v.2.4.0 or v.3.1.5). In the first clustering, we found nine major cell types (primary sensory neurons, SGC, Schwann cells, fibroblasts, VEC, VSMC, VECC, immune cell, and RBC). The cell types were identified by the cell-type-specific/enriched markers genes that have been reported previously. 45 As the number of gene expression had a large difference among different cell types, we further filtered and classified the major cell types one by one, and deleted the inferential doublets (the cells were labeled by two major cell type markers). And the marker genes of each cell cluster were identified using ROC analysis. Finally, data of different batches were integrated using Seurat (v3.1.5). 20 Smart-seq2 library preparation and sequencing We performed the same Smart-seq2 protocol described in the previous publication. 14 Sequencing alignment and quality control for Smart-seq2 data The raw sequencing reads were preprocessed using fastp, 77 and low-quality reads (the length of reads < 35 or the quality of reads < 20) were removed. Then, the clean sequencing reads were mapped to the mouse reference genome (mm10) using HISAT2 (v.2.0.4) protocol. 78,79 Moreover, gene body coverage was performed using RSeQC, 80 and the samples with uniform reads coverage were retained.
Integrating analysis of Smart-seq2 and 10× Genomics We used Seurat 20 (3.0.0) to integrate the Smart-seq2 data and 10× Genomics data with default parameters. Then, the classification of Smart-seq2 data were modified according to neuron size and the expression of marker genes.
Pseudo-time analysis We used monocle 18,19,30 (2.10.1) packages to construct single-cell trajectories to discover the cell types' transitions following nerve injury. We used the DEGs identified by Seurat to sort cells in pseudo-time order. We performed "DDRTree" to reduce dimensions and "plot_complex_cell_trajectory" to plot the minimum spanning tree. "BEAM" was applied to analyze the DEGs depending on each branch, and protein-coding genes were selected for further analysis.
Pseudo-bulk RNA-seq data We split the clear single-cell UMI counts matrix of DRG neurons into individual data. And the pseudo-bulk RNA-seq data were created according to the mean UMI counts of each gene in predicted neuron types.
GO enrichment analysis GO enrichment analysis was performed using Metascape. 81 The protein-coding genes were selected as background.
Inference of regulons and their activity SCENIC was applied to infer the regulons and their target genes in each cell cluster as previously described. 16,82 First, the genes that were not present in the motif database Mouse (mm10, v9 that scored the motifs up to 500 bp upstream of each gene and 10 kb around the transcription start site (TSS)) were filtered out. Next, the retained genes were used for building the co-expression module by GENIE3. Then, TF regulons were identified through the cis-regulatory motif analysis. Finally, TF activity in individual cells was evaluated using AUCell. Cytoscape 83 was applied to visualize gene regulatory networks.
In vivo siRNA injection Before the experiment, each mouse was habited in the behavior room for 2d and the basal mechanical threshold was tested. Then the animals have received 2 consecutive days of intrathecal injection of 4 μg siClcf1 or scramble siRNA a time. The siRNAs were modified with 2′ OME and 5′-Chol. On the third day after the last injection, SNI model was made on the mice. von Frey test was performed in a double-blind way to detect the mechanical threshold at SNI 4d. Then the L4 and L5 DRGs were harvested and the RNA was extracted. RT-PCR was performed to check the efficiency of knockdown.

FACS
The DRGs were dissected from 8 mice of Gal-Cre::Ai9 or Mrgprd-CreER::Ai9 after SNI 2d. The neurons were dissociated into the single-cell suspension and resuspended in the sterilized PBS. They were kept on ice before sorting. The tdTomato-labeled neurons were sorted into the ice-cold DMEM by cell sorter (SONY MA900) and centrifuged at 200× g for 3 min. Then the neurons were planted on the Poly-D-lysine-coated glass dishes with F12 complete medium. siRNA was transfected after 4-5 h waiting for the cell adherent. The neurons were collected 48 h after transfection and the tdTomato + neurons were confirmed under the RFP channel and bright-field channel by fluorescence microscope.
Real-time PCR The total RNA was harvested from the transfected cells and it was reverse transcribed into cDNA with SuperScript II reverse transcriptase (Invitrogen, Cat No. 18064014). The sequence of primers for RT-PCR are listed in Supplementary information, Table S6.
Behavior tests Behavior tests were performed double-blindly. The mechanical threshold was determined by the von Frey test. We used an incremental mechanical stimulation to detect the hind paw responses of the mice. Every stimulation was performed five times around. If the mice withdrew and licked its hind paw with three or more times in a rest, we counted the stimulus is positive. Then the mechanical threshold can be defined.

ISH
We performed ISH as the previous report. 14 Briefly, the probes for the target genes were transcribed and labeled by DIG or FITC in vitro. The templates were amplified from the mouse cDNA by PCR. The primers for the probes were listed in Supplementary information, Table S7. L5 DRG were sectioned into 10 μm slides. Then the slides were fixed with 4% PFA in DEPC-PBS for 20 min, acetylated and prehybridized in hybridization buffer for 3 h at 67°C. The sections were incubated with the hybridization buffer containing the probes (1 μg/mL per probe). For the bright field ISH, the slides were washed, blocked, and incubated with the anti-DIG AP antibodies (1:1000) overnight at 4°C. Then the sections were developed with NBT/BCIP solution in alkaline phosphatase buffer.
As for double fluorescent ISH, 71 the probes were labeled by DIG and FITC, respectively. The pretreatment of slides in FISH on day 1 was the same as the bright field ISH. After washing and blocking, the slides were incubated with anti-FITC-HRP (1:4000) overnight at 4°C. Then the sections were washed with TNT buffer and treated with TSA-plus 2,4-DNP regents (1:100) for 10 min. The sections were incubated with anti-DIG-AP (1:1000), anti-DNP 488 (1:500) overnight at 4°C. Finally, the sections were developed with HNPP/ FR (1:100).
For the FISH combined with immunostaining, the primary antibody was applied at the same time with anti-DIG AP. The sections were incubated with the second antibody for 40 min at 37°C in the last day of this experiment. Then it was washed with PBS and treated with HNPP/FR for 40 min at RT.

Immunostaining
The antibody against galanin (GL Biochem) was obtained from the rabbit, which was immunized with the synthetic antigen, C-PGIPLATSSEDLEKS. To test the antibody specification, the antibodies (1:2000) and galanin (10 −5 M) were mixed and rotated for 1 h at RT and overnight at 4°C. The adult mice used in the experiments were perfused with 4% PFA and picric acid (1:1000). Sections of DRGs were processed for immunostaining. The images were collected with Leica SP8 confocal microscope. Immunolabeling of galanin in DRG neurons was found in the sections incubated with the galanin antibodies, but not with the pre-absorbed antibodies (Supplementary information, Fig. S8c).

RNAscope assay
The probes and detection kit of RNAscope were purchased from ACD. For the bright field assay, the RNAscope 2. ELISA assay The L4/L5 DRGs were collected from the mice of SNI 2d and control group. The ELISA kit was purchased from Novus Biologicals (NBP2-75321). For measurement of CLC secretion from primary sensory neurons, DRG neurons were treated with ECS solution (40 mM KCl) for 2 h after adherence. The supernatant was harvested and centrifugated at 1500 rpm for 20 min, 4°C. Then the proteins in the supernatant were processed as protocol described. The remaining neurons were collected using as lysate control.

Statistical analysis
The data were presented as means ± SEM. Statistical analysis for two groups was performed by Student's t-test. Comparisons between multiple groups were performed by two-way ANOVA with Bonferroni's post hoc test. The difference was considered statistically significant at P < 0.05.

DATA AVAILABILITY
Most datasets generated in the current study are available in the GSE155622. The other datasets analyzed in the current study are included in Li et al. 14 Processed data also available in zhanglab.shbrainproject.com: 100. R scripts are available from the corresponding author upon reasonable request.