The Development of Mechanical Allodynia in Diabetic Rats Revealed by Single-Cell RNA-Seq

Introduction

Diabetic peripheral neuropathy (DPN) is a common complication of diabetes with variable clinical presentations. Patients can suffer from a painless syndrome with loss of sense of touch and temperature, or neuropathic pain manifested by mechanical allodynia (MA) (Berti-Mattera et al., 2008). The latter is the main reason for patients to seek medical advice as it severely debilitates patients' life quality (Quattrini and Tesfaye, 2003). The pathogenesis of painful DPN (PDPN) is complex. Diabetes can lead to multiple damage to the peripheral nervous system, including abnormal activation of multiple metabolic pathways and imbalance of mitochondrial redox state (Feldman et al., 2017). In addition, nerve swelling caused by metabolic disturbances can result in nerve compression injury (Best et al., 2019). In the central nervous system, neuroplastic changes that involve the spinal cord and thalamus were observed in patients or experimental models (Selvarajah et al., 2014; Marshall et al., 2017). The dorsal root ganglia (DRG) contains most of cell bodies of primary sensory neurons, which transmits sensory neural signals through the peripheral nerves to the central nervous system (Maatuf et al., 2019). Due to exposure outside the blood-brain barrier, DRG is particularly vulnerable in diabetes and could be an important trigger for neuropathic pain (Sloan et al., 2018). Overall, the diversity of mechanisms determines that PDPN is difficult to cure. Current first-line treatment for PDPN is nonspecific, and serious adverse effects limit their clinical use (Todorovic, 2015; Snyder et al., 2016). Hence, deeper knowledge in molecular and cellular mechanisms of PDPN is urgently needed to provide a basis for pharmaceutical development.

Single-cell RNA sequencing (scRNA-seq) has progressed at a rapid pace for its advantages in determining tissue heterogeneity and identifying novel cell types (Rodriguez-Meira et al., 2019). Previous transcriptomic studies characterized changes that could contribute to PDPN or DPN in DRG or sural nerves (Hur et al., 2011, 2016; Yamazaki et al., 2013; Guo et al., 2020). These sequenced RNA samples were from whole tissue that inevitably contained mixed cell types. Meanwhile, the “average” gene expression profiles could possibly lead to the confounding interpretation of the results. The difference is that scRNA-seq can sequence thousands of cells in an unbiased manner to uncover both known and novel cell types. We can simultaneously classify and analyze numerous DRG neuron types to know which one is most associated with PDPN. Furthermore, thousands of cell sequencing data would help to investigate whether and how the dynamic DRG neuron type changes in response to hyperglycemia. In the current study, we performed scRNA-seq on rat DRG from PDPN models induced by streptozotocin (STZ) injection (Morrow, 2004). Our work established an unbiased classification of rat DRG neuronal types and identified a novel neuron type MAAC. MAAC was marked by Na, K -ATPase (NKA)-related genes Fxyd7, and Atp1b1. Compared with other neurons, MAAC was characterized by highly expressing genes related to neurofilaments and cytoskeletons. We inferred that it was derived from peptidergic neurons based on correlation analysis and used pseudo-time analysis to identify the gene expression kinetics in neuron type transitions. Furthermore, we studied the transcriptional regulatory networks and communication networks of MAAC. Our results provide comprehensive landscape to uncover the alteration of genes, cell types, and intercellular communication. These results advance our understanding of neuropathic pain in diabetes and provide a basis for the development of pain therapy.

Materials and Methods

Rat

All animal studies were approved by the Shanghai Jiao Tong University Animal Care and Use Committee and conducted in accordance with the animal policies of the Shanghai Jiao Tong University and the guidelines established by the National Health and Family Planning Commission of China. Sprague Dawley rats were obtained from Shanghai Lab Animal Research Center. Rats were housed in specific pathogen-free conditions with free access to water and rat chow.

Diabetes

To induce diabetes, male Sprague-Dawley rats (200–220 g) were injected intraperitoneally with 60 mg/kg STZ (Solarbio, China) dissolved in 1% citrate buffer (pH = 4.5). The control group received equal volumes of the vehicle. Three days after the injections, the glucose concentration was measured in tail vein blood samples using a blood glucose meter (Bayer HealthCare, USA). Only rats with the glucose concentration higher than 11.1 mmol/l were considered diabetic (Dominguez et al., 2015). All the animals were weekly weighed and daily observed during the study.

Mechanical Allodynia

Mechanical nociception was weekly evaluated using the von Frey filaments (North Coast, USA). Rats were separately placed in transparent plexiglass chambers on an elevated metal grid floor with 1 cm holes and probed from below. After a 15-min adaptation period, the hind paws were stimulated with varying forces (in this case 1.4, 2, 4, 6, 8, 10, 15, and 26 grams-force). Every stimulation delivered a constant pre-determined force for 5 s, and the interval of every stimulation was 15 s. The depression in the center of the rats' paws was accurately stimulated. Sensation in this area is within the innervation range of the tibial nerve, of which most fibers originate from L5 DRG in rats (Rigaud et al., 2008; Laedermann et al., 2014). The 50% paw withdrawal threshold (PWT) was determined by the up–down method of Chaplan et al. (1994). The experiment begins by testing the response to the 6 g filament. If a positive response occurred, the next filament with a lower force was applied. Conversely, if a negative response occurred, the next filament with a higher force was applied. This continues until at least four readings are obtained after different reactions appear for the first time. If the force to be used is higher than 15 g or lower than 1.4 g, the threshold on this side is directly recorded as 15 g or 1.4 g, respectively (Deuis et al., 2017). Finally, diabetic rats were divided into groups on day 28 after STZ injection. Diabetic rats with the PWT ≤ 8 g were considered to develop MA. Rats with the PWT ≥ 15 g were selected into the group without MA, and other diabetic rats were excluded from subsequent experiments.

Tissue Dissociation and Cell Purification

In total, 24 rats (eight in each group) were sacrificed 28 days after STZ injection. Among these animals, two control group rats, four diabetic rats with MA, and two diabetic rats without MA were used for scRNA-seq. DRGs, dissected from rats of the abovementioned groups, were washed with Hanks Balanced Salt Solution (HBSS; Sigma-Aldrich, USA) three times. Then, the tissue pieces were transferred into a 15 ml centrifuge tube and digested with 2 ml GEXSCOPE^TM Tissue Dissociation Solution (Singleron, China) at 37°C for 15 min under sustained agitation. After digestion, the samples were filtered using 40 μm sterile strainers and centrifuged at 1,000 rpm for 5 min. The supernatant was discarded and the sediment was resuspended in 1 ml phosphate-buffered saline (PBS; HyClone, USA). Two milliliters of GEXSCOPE^TM red blood cell lysis buffer (Singleron, China) was added to remove the red blood cells for 10 min at 25°C. The solution was then centrifuged at 500 × g for 5 min and suspended in PBS. The sample was stained with trypan blue (Sigma-Aldrich, USA) and microscopically evaluated.

Library Construction and Sequencing

The scRNA-seq libraries were constructed following the protocol of GEXSCOPE^TM Single-Cell RNA Library Kit (Singleron, China) (Dura et al., 2019). Briefly, the single-cell suspensions were loaded onto microfluidic devices. Then, millions of beads with unique cell barcodes and unique molecular identifiers (UMIs) were loaded similar to cells. Only one bead was loaded in each microwell. Following lysis, mRNA was captured onto the beads. More than 95% of the beads were recovered, and the captured mRNA was reverse-transcribed. After cDNA amplification and enrichment, the resulting scRNA-seq libraries were sequenced on Illumina HiSeq X10 instrument with 150 bp paired end reads. Sequencing data of cell populations was conserved in the National Center for Biotechnology Information (NCBI) database with the Gene Expression Omnibus (GEO) accession number GSE176017.

Pre-Processing of scRNA-Seq Data

The raw reads were processed to generate gene expression profiles by an internal pipeline, namely, FastQC for quality evaluation, fastp for trimming, STAR aligner (2.5.3a) for alignment, and featureCounts (1.6.2) for transcript counting (Dobin et al., 2013; Liao et al., 2014; Chen et al., 2018; Wingett and Andrews, 2018). Briefly, after filtering read 1 without polyT tails, cell-barcode and UMI were extracted. Adapters and polyA tails were trimmed before read 2 was mapped to RGSC rn6 reference genome with ensemble version 92 gene annotation (http://www.ensembl.org). Reads with the same cell barcode, UMI, and gene were grouped to calculate the number of UMIs per gene per cell. UMI count tables of each cellular barcode were used for further analysis.

Data Integration and the Dimensionality Reduction

The data were processed by the Seurat (version 4. 0. 4) (Hao et al., 2021) in R software (version 4. 1. 1). According to quality control metrics (Ilicic et al., 2016), cells with genes >2,500 or <200 and the cells that have >10% mitochondrial genes were filtered out. A total of 6,693 cells that were obtained (2,543 in the control group, 1,192 in diabetic rats without MA, and 2,958 in diabetic rats with MA) were used for further bioinformatics analysis. Gene expression of each cell was normalized by total expression, multiplied by a scale factor of 10,000, and log-transformed (NormalizeData function). The 2,000 highly variable genes (HVGs) were selected (FindVariableGenes). Then, we integrated different samples (IntegrateData), and the technical or batch effect was eliminated by canonical correlation analysis (CCA). The expression of each gene that was scaled to shift its mean/variance across cells is 0 and 1 (ScaleData). These results were used as input for dimensionality reduction via principal component analysis (PCA).

Unsupervised Clustering and Cell Type Identification

The top 30 principal components were chosen for cell clustering. The clustering analysis was performed based on FindClusters function after building the nearest neighbor graph using FindNeighbors function. The main cell clusters were identified and visualized with t-distributed stochastic neighbor embedding (t-SNE) plots or uniform manifold approximation and projection (UMAP) plots. Due to differences in algorithms, the distance between two points in the two-dimensional plane of UMAP plots can represent the difference in gene expression information between two cells of the high-dimensional space better. To annotate the cell clusters, cluster biomarkers with high discrimination abilities were identified (FindMarkers function). The cell groups were annotated based on SingleR (Aran et al., 2019) and conventional markers described in previous studies (Fallon, 1985; Arroyo et al., 1998; Le et al., 2005; Dhaka et al., 2007; Schroeter and Steiner, 2009; Usoskin et al., 2015; Li et al., 2016; Oikari et al., 2016; Urban-Ciecko and Barth, 2016; Wu et al., 2017; Donovan et al., 2018; Hockley et al., 2019; Ronning et al., 2019; Zhang et al., 2019; Avraham et al., 2020; Gerber et al., 2021). Similar procedures were applied in the subclusters identification of neurons.

Enrichment Analysis

The gene ontology (GO) enrichment analysis was explored in Database for Annotation, Visualization, and Integrated Discovery (DAVID) (Ashburner et al., 2000).

Pseudo-Time Trajectory Analysis

To discover the gene expression kinetics in neuron type transition in PDPN development, the pseudo-time trajectories were generated with the Monocle package (version 2.18.0) (Trapnell et al., 2014). It can arrange cells along the pseudo-time trajectory, where each cell corresponds to a distinct time point, showing their developmental trajectories such as cell differentiation and other biological processes. The gene expression matrix derived from the Seurat processed data were used as the inputs. The differentially expressed genes (DEGs) were identified (q-value < 0.1) to sort cells in a pseudo-time order. Dimension reduction was performed using the DDRTree method. DEGs were clustered into different categories according to the gene expression patterns along the pseudo-time (Plot_pseudotime_heatmap function). Biological processes were revealed using GO enrichment analysis in DAVID.

Transcription Factor Inference

The analysis of the single-cell gene regulatory network was performed using the single-cell regulatory network inference and clustering (SCENIC) package (Aibar et al., 2017). Rat gene symbols were converted in the corresponding mouse homologous genes using the homologene R package (github.com/oganm/homologene; www.ncbi.nlm.nih.gov/homologene). After initializing settings, the expression matrix was loaded onto GENIE3 for building the initial co-expression gene regulatory networks (GRN). The regulon data were analyzed using the RcisTarget package and the network activity was evaluated. The transcriptional network of TF and predicted target genes were visualized by Cytoscape (Shannon et al., 2003).

Cell Communications Analysis

Cell communications were analyzed based on CellChat (Jin et al., 2021). The ligand-receptor interaction database can be found at http://www.cellchat.org/. All the rat gene names were converted to the mouse ortholog using the homologene R package. The communication network was generated to discover the number of ligand and receptor (L-R) pairs and the signaling strength between each cell cluster. The strength was defined by the communication probability, which was calculated by a specific algorithm using the geometric mean of the normalized expression of receptors and ligands as input data. Outgoing and incoming interaction strength was calculated to compare the communication ability of each neuronal cluster as senders and receivers. The communication probabilities mediated by L-R pairs were compared and visualized in the bubble plot.

Statistical Analysis

All statistical analyses and the graph generation were performed in R (version 4. 1. 1).

Code Availability

The analysis codes are available on GitHub: https://github.com/SJTU-ZhouHan/SJTU-ZhouHan.

Result

Cellular Constitution of Rat DRG

In the present work, we intended to investigate the cell diversity of DRG neurons in the adult rat under PDPN conditions. A conspicuous MA typically becomes steady 4 weeks after diabetes induction in rats. We performed scRNA-seq on the cells dissociated from bilateral lumbar (L) 5 DRGs of diabetic rats with and without MA and normal control groups (Figure 1A; Supplementary Figure 1). After quality control, 6,693 cells, including 3,979 neurons, were obtained. Visualization of single-cell transcriptomes in t-SNE space was able to separate cells into clusters which we mapped to eight major cell types with distinct markers (Figure 1B). These cells included neurons, satellite glial cells (SGC), proliferating SGC (PSGC), Schwann cells, fibroblasts, vascular smooth muscle cells (VSMC), vascular endothelial cells (VEC), and microglia. SGC were identified due to highly expressed Fabp7 (Figure 1C) (Goncalves et al., 2017). PSGC specifically expressed Fabp7 and proliferation markers Top2a (Figure 1C) (Gerber et al., 2021). The cell population that highly expressed Mpz was annotated as Schwann cells (Figure 1C) (Le et al., 2005). Microglia highly expressed Lyz2 (Figure 1C) (Donovan et al., 2018; Ronning et al., 2019). In addition, we identified VEC marked by Cldn5 and VSMC by Tpm2 (Figure 1C) (Jang et al., 2011; Kalluri et al., 2019).

FIGURE 1

Figure 1. The heterogeneity of dorsal root ganglia (DRG) cells in the painful diabetic peripheral neuropathy (PDPN) model. (A) Workflow of the pain threshold evaluation, sample preparation, sequencing, and bioinformatic analysis. The red region in the center of the rats' paws was stimulated by von Frey filaments. Sensation in this area is within the innervation range of the tibial nerve, originating from L5 DRG in rat. (B) t-distributed scholastic neighbor embedding (t-SNE) plot of single cells profiled in the presenting work colored by cell types. Each colored dot represents a cell. (C) Feature heatmap shows the marker genes in each cell type. The color represents the gene expression level after batch effect correction and normalization. (D) t-SNE plot shows the unsupervised clustering DRG neurons. Dots, individual cells; colors, neuron clusters. (E) Dot plot shows the differentially expressed genes (DEGs) of each DRG neuron cluster. The size of the dot means the percentage of cells expressing the gene, and the color indicates the average expression level.

Subclusters-Specific Analysis of Neurons

To detect discrete neuron subclusters, we reclassified neurons based on the canonical DRG neuron markers (Usoskin et al., 2015; Liguz-Lecznar et al., 2016; Kupari et al., 2021). Eleven neuron clusters and their DEGs were unbiasedly identified by Seurat program (Figures 1D,E). These clusters were annotated as non-peptidergic nociceptors (NP), peptidergic nociceptors (PEP), somatostatin-positive neurons (SOM), C-fiber low-threshold mechanoreceptors (C-LTMR), and Trpm8-positive neurons (TRPM8) (Figure 2A). We identified MA-associated clusters (MAAC) for the obviously increased cell number in diabetic rats with MA (Figure 2A). The specific neuron number for each cluster is given in Supplementary Table 1. Neuron clusters were matched with a previous scRNA-seq study of macaques DRG (Figure 2B) (Kupari et al., 2021). NP was named based on Purinergic receptor P2X3 (P2rx3), and PEP was named based on tropomyosin receptor kinase A (TRKA, Ntrk1). SOM was identified using somatostatin (Sst) and interleukin 31 receptor A (Il31ra). Although SOM contains markers for both peptidergic and nonpeptidergic neurons (Tac1 and P2rx3), this cluster was named NP3 in previous literature (Usoskin et al., 2015). Exoc1l, P2ry1, and Gfra2 were used to assign C-LTMR. Trpm8 expression was used for naming the TRPM8 neuron cluster.

FIGURE 2

Figure 2. Novel neuron cluster mechanical allodynia-associated cluster (MAAC) appeared in the development of PDPN. (A) t-SNE plot shows somatosensory neuron clusters in control versus in diabetic rats without or with mechanical allodynia (MA). The dots in the control group are grayed out when they are compared with the dots in the other two groups. Dots, individual cells; colors, neuron clusters. (B) Summary of neuron type classification in this work and previous study. (C) Heatmap shows the DEGs of MAAC. The genes with the largest fold difference are marked. Color, expression level. (D) Bar plot shows gene ontology (GO) terms of biological processes [false discovery rate (FDR) <0.3] enriched for the DEGs of MAAC.

Emerging Novel Neuron Cluster Associated With MA

The t-SNE plot was split based on different groups to present the progression of the emerging novel neuron cluster MAAC (Figure 2A). Neurons in other clusters were evenly distributed in different groups (Figure 2A). Hundreds of DEGs of MAAC were identified (log2FC > 0.25 and adj. p-value < 0.05) (Figure 2C; Supplementary Table 2). Among these genes, Fxyd7 and Atp1b1 were highly expressed in MAAC. The paralog of Fxyd7, Fxyd1, and the paralog of Atp1b1 and Atp2b4 were also included (Supplementary Table 2). We next annotated the DEGs with biological processes [false discovery rate (FDR) < 0.3] from GO databases to obtain transcriptomic characteristics of MAAC. These cell-type specific genes recapitulated neurofilament, axon, synapse, and receptor-related biological processes (Figure 2D; Supplementary Table 3). “Regulation of cell shape,” “receptor internalization,” and “neurofilament bundle assembly” were the most significant enrichment terms. “Neurofilament bundle assembly” was the most enriched biological process with the highest fold enrichment.

The Origin and Transition of MAAC

To explore the origin of MAAC, we calculated the transcriptomic correlations of different neuron clusters. The heatmap showed that the gene expression of MAAC was close to PEP (Figure 3A), and the distance between these two clusters on the UMAP plot also provides corroborative evidence (Supplementary Figure 2). Then, we reconstructed the pseudo-time trajectories and identify genes whose expression changed as the neurons underwent transition (Figures 3B,C). Genes in the neuronal switch process from PEP to MAAC were classified into four modules based on their expressing patterns (Figure 3C). GO enrichment analysis were used to analyze the biological processes of each gene module. Upregulated genes in Module 1 were associated with “cell morphogenesis,” including Tgfb2, Tenm4, Hgf , and Atp2b2. Module 2 reflected the genes that were upregulated from a lower level across pseudo-time. These genes were enriched in “cell adhesion” and “regulation of neuron projection development.” Comparatively, the expression of genes related to “inflammatory response” eventually decreased. Notably, some genes involved in “sensory perception of pain” were downregulated, including Ndn, Trpv1, Asic1, Aqp1 (Module 3), Adcyap1, Trpa1, Npy1r, Oprk1, Tac1 (Module 4).

FIGURE 3

Figure 3. The origin and transition of MAAC. (A) Heatmap shows the Pearson correlation of each neuron cluster based on their genes expression files. (B) Pseudo-time trajectories show the fate of MAAC originated from peptidergic nociceptors (PEP). Dot, cells; colors, neuron clusters or pseudotime. (C) Heatmaps show the DEGs clustered based on their dynamic expression characteristic, which were shown in pseudo-time from PEP to MAAC.

Pivotal Regulons Involved in Neuron Type Transition

The SCENIC analysis was performed to identify the critical regulators and to reconstruct gene regulatory network. The binary regulon activity matrix was generated and two MAAC-specific regulons, Hobx7 and Lar1p, were identified (Figure 4A). The binary t-SNE plot showed that these TFs were expressed and activated in MAAC (Figure 4B). Other TFs like Fosb, Cebpd, Fosl1, Junb, and Ddit3 had broad expression patterns. They were mainly expressed and activated in PEP and MAAC. In addition, we constructed the gene regulatory network of Hobx7 and Lar1p (Figure 4C). We observed that Hobx7 regulated the most DEGs of MAAC, including genes with high expression level, like Atp1b1 and S100b. However, Larp1 only regulated the expression of Gap43 and Arid5b. Thus, Hobx7 could serve as a critical regulator in the formation of MAAC.

FIGURE 4

Figure 4. Specifically expressed transcription factors and regulatory network of MAAC. (A) The binary heatmap shows the activity of inferred transcription factors in different neuron clusters. (B) t-SNE plots show the activities cells distribution of two specific expressed regulons of MAAC. Dots, individual cells; colors, activated cells. (C) Gene regulatory networks of Hobx7 and Larp1 inferred by single-cell regulatory network inference and clustering (SCENIC). Colors indicated the log2FC of DEGs of MAAC.

Intercellular Communication Analysis of MAAC

The intercellular communication was inferred based on the information of L-R pairs in CellChat database (Jin et al., 2021). The numbers of interactions between cell clusters were calculated (Figures 5A,B). There is a close link between MAAC and SGC, but not between other neuron clusters (Figures 5A,B). Then, we compared the outgoing and incoming interaction strength of each cell type (Figure 5C). The interaction strength of glia appeared to be greater than other cells, and the strength of MAAC was obviously greater than their PEP cluster origin. To identify the ligands or receptors involved, we respectively compared the communication probabilities mediated by L-R pairs between MAAC and other cells (Figures 5D,E). Only the L-R pairs with the most pronounced changes (p-value < 0.05) were displayed. When MAAC served as the signal source, communication probabilities of PTN ligand were prominent. MAAC may communicate with SGC and PSGC via Ptn-Sdc4, Ptn-Ptprz1, Ptn-Ncl, and Bdnf-Ntrk2. The communication between MAAC and VEC was mediated by Calcb-Calcrl and Calca-Calcrl in low communication probabilities. Besides, MAAC could communicate with microglia and VSMC by Ptn-Ncl. When MAAC served as the signal target, the Ncl receptor was playing an important role. The biologic behavior of MAAC may be regulated by Ptn-Ncl and Mdk-Ncl.

FIGURE 5

Figure 5. The intercellular communication of MAAC. (A) Intercellular communication network reflects the number of interactions between clusters. Nodes, neuronal clusters; node size, cell counts; edge width, number of interactions. (B) Intercellular communication network shows the number of interactions between MAAC and other neuron clusters. Nodes, neuronal clusters; node size, cell counts; edge width, number of interactions. (C) Dot plot shows the incoming and outgoing strength of each cell type. The outgoing /incoming interaction strength is defined by the comprehensive communication probability between the signal sending /target cells and all cell types. (D) Dot plot shows the communication probabilities of ligand-receptor pairs when MAAC are sender cells. Only the ligand-receptor pairs with significant changes (p <0.05) are shown in the picture. Commun. Prob, communication probability. The communication probability here equals to the interaction strength and is not exactly a probability. Dot color means communication probabilities and dot size represents computed p-values. (E) Dot plot shows the communication probabilities of ligand-receptor pairs when MAAC serves as target cells.

Discussion

Dorsal root ganglia are vulnerable in diabetes due to lack of the protection of blood-nerve barrier. Injured somatosensory neurons could be an important trigger for pain. Here, using scRNA-seq, we generated a cell-type specific transcriptome atlas and revealed the complexity of the cellular landscape inside DRG tissues. Then, we identified a PDPN-related cluster MAAC and their origin. The transcriptomic characteristics of neuron type transition and pivotal regulons were demonstrated. Finally, we constructed the intercellular communication networks and revealed associated L-R pairs.

In the current study, we used STZ-induced diabetic rats to evoke MA. The extent of nerve damage exhibited individual differences among diabetic rodent models. In detail, in the same batch, the PWT of some rats dropped below 1.4 g, while some still stayed above 15 g on Day 28 after STZ injection. We compared L5 DRG cells from mechanically sensitive (diabetes with MA; PWT ≤ 8 g), insensitive rats (diabetes without MA; PWT ≥ 15 g), and normal rats.

Traditionally, DRG neuron classification was based on size and neuron-type genetic marker, but they were not always consistent with the functional heterogeneity of DRG neurons (Murray and Cheema, 2003). In contrast, comprehensive transcriptome analysis of scRNA-seq data can classify neurons in an unbiased manner. Usoskin et al. (2015) clustered mice L4-L6 DRG neurons into PEP, NP, and neurofilament-containing and tyrosine hydroxylase-containing clusters. Subsequently, these clusters were further classified into 11 subtypes including thermosensitive, itch sensitive, low-threshold mechanosensitive, and nociceptive neurons, etc. Compared with this study, we also identified cold thermoreceptors (TrpM8), C-LTMR, and SOM neurons which expressed the itch-related biomarker Il31ra. The SOM was named NP3 in the study of Usoskin et al. (2015). Because it expressed both peptidergic and non-peptidergic neuronal markers, and the distance in space is relatively independent in the t-SNE plot, we think it is more reasonable to label it as SOM. Kupari et al. (2021) classified macaques' primate DRG sensory neurons. The results of their neuron subtypes identification were broadly consistent Usoskin et al.'s (2015). A difference was that Nefh⁺ neurons were not clustered into one subcluster as NF cluster in the latter, which was closer to our classification result.

In previous single-cell or single-nucleus sequencing studies of different models of pain, injury-related subpopulations of neurons have been reported (Nguyen et al., 2019; Renthal et al., 2020; Wang et al., 2021). In the present study, we found that a new neuron type appeared in the development of PDPN, namely, MAAC (Fxyd7⁺/Atp1b1⁺), and identified the transcriptomic characteristics, origin, transition, and intercellular communication of them. Two DEGs with the highest fold change were considered as the biomarkers of MAAC. Interestingly, they are both closely associated with the NKA, and the dysfunction of NKA has been reported in patients with DPN (Vague et al., 1997; Krishnan et al., 2008). Fxyd7 proven to be a regulator of NKA isozymes, which could significantly affect the apparent affinity for extracellular K⁺ of NKA α1–β1 complexes (Béguin et al., 2002). The paralog of Fxyd7 and Fxyd2 could regulate neuronal activity by modulating NKA activity, which may be a fundamental mechanism underlying the persistent hypersensitivity to pain (Wang et al., 2015). Atp1b1 encodes the protein of NKA family. Its paralog, Atp1b3, has been reported to be involved in the formalin pain behavior (LaCroix-Fralish et al., 2009).

The enrichment analysis investigated multiple processes of MAAC. Some processes with high significance were involved in neurofilament and cytoskeleton, including “regulation of cell shape,” “neurofilament bundle assembly,” and “neurofilament cytoskeleton organization.” The mechanisms of how the cytoskeleton affects mechanical nociception are not fully understood, but it is known that through cytoskeleton force, it can be transmitted to cell membrane receptors, making them more responsive to pressure (Ingber, 2006). Bhattacherjee et al. (2017) reported that ganglion cytoskeletal genes were critical in determining mechanosensory properties in Rett syndrome. In addition, Dina et al. (2003) suggested that inflammatory mediator-induced MA was differentially dependent on the neuron cytoskeleton and that cytoskeletal disruptors could attenuate epinephrine-induced hyperalgesia in rat paws. Thus, activated cytoskeletal genes are another important feature of MAAC.

Our results suggested MAAC could originate from PEP by the alteration of gene expression profiles. The similar cellular transition was found in spared nerve injury (SNI) models. Wang et al. (2021) found three novel neuronal clusters induced by SNI. Most of them originate from the peptidergic neuron cluster. Then, through tracing the inferred trajectories, the dynamic characteristics were revealed. Notably, some enriched genes in sensory perception of pain were downregulated during MAAC formation. A similar phenomenon was also present after SNI injury (Wang et al., 2021). Further, based on the binary heatmap of regulon activity, we found the fate of MAAC could be determined by Hobx7 and Larp1. Hobx7 played a dominant role and was involved in the regulation of many DEGs of MAAC.

Next, the cellular communication activity of MAAC was inferred, and we found abundant communication relationships between MAAC and SGC and endothelial cells and microglia. Soma of neurons could not form synaptic contact with one another (Pannese, 1981), and each soma is tightly enwrapped by SGC. Therefore, the interactions of MAAC with other neurons were not found, and how hyperglycemia influences the signaling between SGC and neurons should be of particular concern. Both MAAC and SGC can send or receive signals to or from each other via Ptn-Ncl. SGC can also influence MAAC via Mdk-Ncl. Ncl encodes nucleolin, which is involved in the synthesis and maturation of ribosomes (Hirano et al., 2005). Previous studies have pointed out that nucleolin in neuronal cell bodies can restrict axon growth (Perry et al., 2016). Hence, we hypothesize that SGC may regulate abnormal axon growth through Ptn-Ncl and Mdk-Ncl in PDPN.

In summary, the present study reveals the single-cell transcriptomic alterations of somatosensory neurons in PDPN, identifies a novel neuron type MAAC, and investigates the novel neuron type's transcriptomic characteristics, origin, transition trajectory, regulons, and cellular communication. Thus, the study's findings advanced our current understanding of PDPN and neuropathic pain, which can be a resource for developing pain therapy. However, there were some limitations in the present study. For example, further basic experiments are needed to validate our results, and how the novel cluster MAAC induces pain symptom remains to be further clarified.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Ethics Statement

The animal study was reviewed and approved by the Ethics Committee, Ninth People Hospital Affiliated to Shanghai Jiao Tong University School of Medicine.

Author Contributions

HZ and XY built models, collected specimens, analyzed data, and wrote the manuscript. CL and HC analyzed data. YW, BX, and FM contributed to the discussion. WZ reviewed the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by the National Natural Science Foundation of China (Grant No. 81771320).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We sincerely thank Dr. Jie Yu (Ninth People Hospital, Shanghai Jiao Tong University School of Medicine) for helpful discussion and Dr. Peng Jin (Ruijin Hospital, Shanghai Jiao Tong University School of Medicine) for his bioinformatics advice.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnmol.2022.856299/full#supplementary-material

References

Aibar, S., González-Blas, C. B., Moerman, T., Huynh-Thu, V. A., Imrichova, H., Hulselmans, G., et al. (2017). SCENIC: single-cell regulatory network inference and clustering. Nat Methods. 14, 1083–1086. doi: 10.1038/nmeth.4463

PubMed Abstract | CrossRef Full Text | Google Scholar

Aran, D., Looney, A. P., Liu, L., Wu, E., Fong, V., Hsu, A., et al. (2019). Reference-based analysis of lung single-cell sequencing reveals a transitional profibrotic macrophage. Nat Immunol. 20, 163–172. doi: 10.1038/s41590-018-0276-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Arroyo, E. J., Bermingham, J. R. Jr., Rosenfeld, M. G., and Scherer, S. S. (1998). Promyelinating Schwann cells express Tst-1/SCIP/Oct-6. J Neurosci. 18, 7891–7902. doi: 10.1523/jneurosci.18-19-07891.1998

PubMed Abstract | CrossRef Full Text | Google Scholar

Ashburner, M. M., Ball, C. A. C., Blake, J. A. J., Botstein, D. D., Butler, H. H., Cherry, J. M. J., et al. (2000). Gene ontology: tool for the unification of biology Consortium TGO. Nat Genet. 25, 25–29. doi: 10.1038/75556

PubMed Abstract | CrossRef Full Text | Google Scholar

Avraham, O., Deng, P. Y., Jones, S., Kuruvilla, R., Semenkovich, C. F., Klyachko, V. A., et al. (2020). Satellite glial cells promote regenerative growth in sensory neurons. Nat Commun. 11, 4891. doi: 10.1038/s41467-020-18642-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Béguin, P., Crambert, G., Monnet-Tschudi, F., Uldry, M., Horisberger, J. D., Garty, H., et al. (2002). FXYD7 is a brain-specific regulator of Na,K-ATPase alpha 1-beta isozymes. Embo j. 21, 3264–3273. doi: 10.1093/emboj/cdf330

PubMed Abstract | CrossRef Full Text | Google Scholar

Berti-Mattera, L. N., Kern, T. S., Siegel, R. E., Nemet, I., and Mitchell, R. (2008). Sulfasalazine blocks the development of tactile allodynia in diabetic rats. Diabetes. 57, 2801–2808. doi: 10.2337/db07-1274

PubMed Abstract | CrossRef Full Text | Google Scholar

Best, T. J., Best, C. A., Best, A. A., and Fera, L. A. (2019). Surgical peripheral nerve decompression for the treatment of painful diabetic neuropathy of the foot—A level 1 pragmatic randomized controlled trial. Diabetes Res Clin Pract. 147, 149–156. doi: 10.1016/j.diabres.2018.08.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Bhattacherjee, A., Mu, Y., Winter, M. K., Knapp, J. R., Eggimann, L. S., Gunewardena, S. S., et al. (2017). Neuronal cytoskeletal gene dysregulation and mechanical hypersensitivity in a rat model of Rett syndrome. Proc Natl Acad Sci USA. 114, E6952–e6961. doi: 10.1073/pnas.1618210114

PubMed Abstract | CrossRef Full Text | Google Scholar

Chaplan, S. R., Bach, F. W., Pogrel, J. W., Chung, J. M., and Yaksh, T. L. (1994). Quantitative assessment of tactile allodynia in the rat paw. J Neurosci Methods. 53, 55–63. doi: 10.1016/0165-0270(94)90144-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, S., Zhou, Y., Chen, Y., and Gu, J. (2018). fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics. 34, i884–i890. doi: 10.1093/bioinformatics/bty560

PubMed Abstract | CrossRef Full Text | Google Scholar

Deuis, J. R., Dvorakova, L. S., and Vetter, I. (2017). Methods used to evaluate pain behaviors in rodents. Front Mol Neurosci. 10, 284. doi: 10.3389/fnmol.2017.00284

PubMed Abstract | CrossRef Full Text | Google Scholar

Dhaka, A., Murray, A. N., Mathur, J., Earley, T. J., Petrus, M. J., and Patapoutian, A. (2007). TRPM8 is required for cold sensation in mice. Neuron. 54, 371–378. doi: 10.1016/j.neuron.2007.02.024

PubMed Abstract | CrossRef Full Text | Google Scholar

Dina, O. A., McCarter, G. C., de Coupade, C., and Levine, J. D. (2003). Role of the sensory neuron cytoskeleton in second messenger signaling for inflammatory pain. Neuron. 39, 613–624. doi: 10.1016/s0896-6273(03)00473-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Dobin, A., Davis, C. A., Schlesinger, F., Drenkow, J., Zaleski, C., Jha, S., et al. (2013). STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 29, 15–21. doi: 10.1093/bioinformatics/bts635

PubMed Abstract | CrossRef Full Text | Google Scholar

Dominguez, J. M 2nd, Yorek, M. A., and Grant, M. B. (2015). Combination therapies prevent the neuropathic, proinflammatory characteristics of bone marrow in streptozotocin-induced diabetic rats. Diabetes. 64, 643–653. doi: 10.2337/db14-0433

PubMed Abstract | CrossRef Full Text | Google Scholar

Donovan, K. M., Leidinger, M. R., McQuillen, L. P., Goeken, J. A., Hogan, C. M., Harwani, S. C., et al. (2018). Allograft inflammatory factor 1 as an immunohistochemical marker for macrophages in multiple tissues and laboratory animal species. Comp Med. 68, 341–348. doi: 10.30802/aalas-cm-18-000017

PubMed Abstract | CrossRef Full Text | Google Scholar

Dura, B., Choi, J. Y., Zhang, K., Damsky, W., Thakral, D., Bosenberg, M., et al. (2019). scFTD-seq: freeze-thaw lysis based, portable approach toward highly distributed single-cell 3' mRNA profiling. Nucleic Acids Res. 47, e16. doi: 10.1093/nar/gky1173

PubMed Abstract | CrossRef Full Text | Google Scholar

Fallon, J. R. (1985). Neurite guidance by non-neuronal cells in culture: preferential outgrowth of peripheral neurites on glial as compared to nonglial cell surfaces. J Neurosci. 5, 3169–3177. doi: 10.1523/jneurosci.05-12-03169.1985

PubMed Abstract | CrossRef Full Text | Google Scholar

Feldman, E. L., Nave, K. A., Jensen, T. S., and Bennett, D. L. H. (2017). New horizons in diabetic neuropathy: mechanisms, bioenergetics, and pain. Neuron. 93, 1296–1313. doi: 10.1016/j.neuron.2017.02.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Gerber, D., Pereira, J. A., Gerber, J., Tan, G., Dimitrieva, S., Yángüez, E., et al. (2021). Transcriptional profiling of mouse peripheral nerves to the single-cell level to build a sciatic nerve ATlas (SNAT). Elife.10doi: 10.7554/eLife.58591

PubMed Abstract | CrossRef Full Text | Google Scholar

Goncalves, N. P., Vaegter, C. B., Andersen, H., Ostergaard, L., Calcutt, N. A., and Jensen, T. S. (2017). Schwann cell interactions with axons and microvessels in diabetic neuropathy. Nat Rev Neurol. 13, 135–147. doi: 10.1038/nrneurol.2016.201

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, K., Eid, S. A., Elzinga, S. E., Pacut, C., Feldman, E. L., and Hur, J. (2020). Genome-wide profiling of DNA methylation and gene expression identifies candidate genes for human diabetic neuropathy. Clin Epigenetics. 12, 123. doi: 10.1186/s13148-020-00913-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Hao, Y., Hao, S., Andersen-Nissen, E., Mauck, W. M. 3rd, Zheng, S., Butler, A., Lee, M. J., et al. (2021). Integrated analysis of multimodal single-cell data. Cell. 184, 3573-3587.e3529. doi: 10.1016/j.cell.2021.04.048

PubMed Abstract | CrossRef Full Text | Google Scholar

Hirano, K., Miki, Y., Hirai, Y., Sato, R., Itoh, T., Hayashi, A., et al. (2005). A multifunctional shuttling protein nucleolin is a macrophage receptor for apoptotic cells. J Biol Chem. 280, 39284–39293. doi: 10.1074/jbc.M505275200

PubMed Abstract | CrossRef Full Text | Google Scholar

Hockley, J. R. F., Taylor, T. S., Callejo, G., Wilbrey, A. L., Gutteridge, A., Bach, K., et al. (2019). Single-cell RNAseq reveals seven classes of colonic sensory neuron. Gut. 68, 633–644. doi: 10.1136/gutjnl-2017-315631

PubMed Abstract | CrossRef Full Text | Google Scholar

Hur, J., O'Brien, P. D., Nair, V., Hinder, L. M., McGregor, B. A., Jagadish, H. V., et al. (2016). Transcriptional networks of murine diabetic peripheral neuropathy and nephropathy: common and distinct gene expression patterns. Diabetologia. 59, 1297–1306. doi: 10.1007/s00125-016-3913-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Hur, J., Sullivan, K. A., Pande, M., Hong, Y., Sima, A. A., Jagadish, H. V., et al. (2011). The identification of gene expression profiles associated with progression of human diabetic neuropathy. Brain. 134, 3222–3235. doi: 10.1093/brain/awr228

PubMed Abstract | CrossRef Full Text | Google Scholar

Ilicic, T., Kim, J. K., Kolodziejczyk, A. A., Bagger, F. O., McCarthy, D. J., Marioni, J. C., et al. (2016). Classification of low quality cells from single-cell RNA-seq data. Genome Biol. 17, 29. doi: 10.1186/s13059-016-0888-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Ingber, D. E. (2006). Cellular mechanotransduction: putting all the pieces together again. Faseb j. 20, 811–827. doi: 10.1096/fj.05-5424rev

PubMed Abstract | CrossRef Full Text | Google Scholar

Jang, A. S., Concel, V. J., Bein, K., Brant, K. A., Liu, S., Pope-Varsalona, H., et al. (2011). Endothelial dysfunction and claudin 5 regulation during acrolein-induced lung injury. Am J Respir Cell Mol Biol. 44, 483–490. doi: 10.1165/rcmb.2009-0391OC

PubMed Abstract | CrossRef Full Text | Google Scholar

Jin, S., Guerrero-Juarez, C. F., Zhang, L., Chang, I., Ramos, R., Kuan, C. H., et al. (2021). Inference and analysis of cell-cell communication using CellChat. Nat Commun. 12, 1088. doi: 10.1038/s41467-021-21246-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Kalluri, A. S., Vellarikkal, S. K., Edelman, E. R., Nguyen, L., Subramanian, A., Ellinor, P. T., et al. (2019). Single-cell analysis of the normal mouse aorta reveals functionally distinct endothelial cell populations. Circulation. 140, 147–163. doi: 10.1161/circulationaha.118.038362

PubMed Abstract | CrossRef Full Text | Google Scholar

Krishnan, A. V., Lin, C. S., and Kiernan, M. C. (2008). Activity-dependent excitability changes suggest Na+/K+ pump dysfunction in diabetic neuropathy. Brain. 131, 1209–1216. doi: 10.1093/brain/awn052

PubMed Abstract | CrossRef Full Text | Google Scholar

Kupari, J., Usoskin, D., Parisien, M., Lou, D., Hu, Y., Fatt, M., et al. (2021). Single cell transcriptomics of primate sensory neurons identifies cell types associated with chronic pain. Nat Commun. 12, 1510. doi: 10.1038/s41467-021-21725-z

PubMed Abstract | CrossRef Full Text | Google Scholar

LaCroix-Fralish, M. L., Mo, G., Smith, S. B., Sotocinal, S. G., Ritchie, J., Austin, J. S., et al. (2009). The beta3 subunit of the Na+,K+-ATPase mediates variable nociceptive sensitivity in the formalin test. Pain. 144, 294–302. doi: 10.1016/j.pain.2009.04.028

PubMed Abstract | CrossRef Full Text | Google Scholar

Laedermann, C. J., Pertin, M., Suter, M. R., and Decosterd, I. (2014). Voltage-gated sodium channel expression in mouse DRG after SNI leads to re-evaluation of projections of injured fibers. Mol Pain. 10, 19. doi: 10.1186/1744-8069-10-19

PubMed Abstract | CrossRef Full Text | Google Scholar

Le, N., Nagarajan, R., Wang, J. Y., Araki, T., Schmidt, R. E., and Milbrandt, J. (2005). Analysis of congenital hypomyelinating Egr2Lo/Lo nerves identifies Sox2 as an inhibitor of Schwann cell differentiation and myelination. Proc Natl Acad Sci USA. 102, 2596–2601. doi: 10.1073/pnas.0407836102

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, C. L., Li, K. C., Wu, D., Chen, Y., Luo, H., Zhao, J. R., et al. (2016). Somatosensory neuron types identified by high-coverage single-cell RNA-sequencing and functional heterogeneity. Cell Res. 26, 83–102. doi: 10.1038/cr.2015.149

PubMed Abstract | CrossRef Full Text | Google Scholar

Liao, Y., Smyth, G. K., and Shi, W. (2014). featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics. 30, 923–930. doi: 10.1093/bioinformatics/btt656

PubMed Abstract | CrossRef Full Text | Google Scholar

Liguz-Lecznar, M., Urban-Ciecko, J., and Kossut, M. (2016). Somatostatin and somatostatin-containing neurons in shaping neuronal activity and plasticity. Front Neural Circuits. 10, 48. doi: 10.3389/fncir.2016.00048

PubMed Abstract | CrossRef Full Text | Google Scholar

Maatuf, Y., Geron, M., and Priel, A. (2019). The role of toxins in the pursuit for novel analgesics. Toxins (Basel) 0.11. doi: 10.3390/toxins11020131

PubMed Abstract | CrossRef Full Text | Google Scholar

Marshall, A. G., Lee-Kubli, C., Azmi, S., Zhang, M., Ferdousi, M., Mixcoatl-Zecuatl, T., et al. (2017). Spinal disinhibition in experimental and clinical painful diabetic neuropathy. Diabetes. 66, 1380–1390. doi: 10.2337/db16-1181

PubMed Abstract | CrossRef Full Text | Google Scholar

Morrow, T. J. (2004). Animal models of painful diabetic neuropathy: the STZ rat model. Curr. Protoc. Neurosci. Chapter 9, Unit 9.18. doi: 10.1002/0471142301.ns0918s29

PubMed Abstract | CrossRef Full Text | Google Scholar

Murray, S. S., and Cheema, S. S. (2003). Constitutive expression of the low-affinity neurotrophin receptor and changes during axotomy-induced death of sensory neurones in the neonatal rat dorsal root ganglion. J Anat. 202, 227–238. doi: 10.1046/j.1469-7580.2003.00151.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Nguyen, M. Q., Le Pichon, C. E., and Ryba, N. (2019). Stereotyped transcriptomic transformation of somatosensory neurons in response to injury. Elife. 8:e49679. doi: 10.7554/eLife.49679

PubMed Abstract | CrossRef Full Text | Google Scholar

Oikari, L. E., Okolicsanyi, R. K., Qin, A., Yu, C., Griffiths, L. R., and Haupt, L. M. (2016). Cell surface heparan sulfate proteoglycans as novel markers of human neural stem cell fate determination. Stem Cell Res. 16, 92–104. doi: 10.1016/j.scr.2015.12.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Pannese E. (1981). The satellite cells of the sensory ganglia. Adv. Anat. Embryol. Cell. Biol. 65, 1–111. doi: 10.1007/978-3-642-67750-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Perry, R. B., Rishal, I., Doron-Mandel, E., Kalinski, A. L., Medzihradszky, K. F., Terenzio, M., et al. (2016). Nucleolin-mediated RNA localization regulates neuron growth and cycling cell size. Cell Rep. 16, 1664–1676. doi: 10.1016/j.celrep.2016.07.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Quattrini, C., and Tesfaye, S. (2003). Understanding the impact of painful diabetic neuropathy. Diabetes Metab Res Rev. 19, S2–8. doi: 10.1002/dmrr.360

PubMed Abstract | CrossRef Full Text | Google Scholar

Renthal, W., Tochitsky, I., Yang, L., Cheng, Y. C., Li, E., Kawaguchi, R., et al. (2020). Transcriptional reprogramming of distinct peripheral sensory neuron subtypes after axonal injury. Neuron. 108, 128-144.e129. doi: 10.1016/j.neuron.2020.07.026

PubMed Abstract | CrossRef Full Text | Google Scholar

Rigaud, M., Gemes, G., Barabas, M. E., Chernoff, D. I., Abram, S. E., Stucky, C. L., et al. (2008). Species and strain differences in rodent sciatic nerve anatomy: implications for studies of neuropathic pain. Pain. 136, 188–201. doi: 10.1016/j.pain.2008.01.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Rodriguez-Meira, A., Buck, G., Clark, S. A., Povinelli, B. J., Alcolea, V., Louka, E., et al. (2019). Unravelling intratumoral heterogeneity through high-sensitivity single-cell mutational analysis and parallel RNA sequencing. Mol Cell. 73, 1292–1305.e1298. doi: 10.1016/j.molcel.2019.01.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Ronning, K. E., Karlen, S. J., Miller, E. B., and Burns, M. E. (2019). Molecular profiling of resident and infiltrating mononuclear phagocytes during rapid adult retinal degeneration using single-cell RNA sequencing. Sci Rep. 9, 4858. doi: 10.1038/s41598-019-41141-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Schroeter, M. L., and Steiner, J. (2009). Elevated serum levels of the glial marker protein S100B are not specific for schizophrenia or mood disorders. Mol Psychiatry. 14, 235–237. doi: 10.1038/mp.2008.85

PubMed Abstract | CrossRef Full Text | Google Scholar

Selvarajah, D., Wilkinson, I. D., Maxwell, M., Davies, J., Sankar, A., Boland, E., et al. (2014). Magnetic resonance neuroimaging study of brain structural differences in diabetic peripheral neuropathy. Diabetes Care. 37, 1681–1688. doi: 10.2337/dc13-2610

PubMed Abstract | CrossRef Full Text | Google Scholar

Shannon, P., Markiel, A., Ozier, O., Baliga, N. S., Wang, J. T., Ramage, D., et al. (2003). Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 13, 2498–2504. doi: 10.1101/gr.1239303

PubMed Abstract | CrossRef Full Text | Google Scholar

Sloan, G., Shillo, P., Selvarajah, D., Wu, J., Wilkinson, I. D., Tracey, I., et al. (2018). A new look at painful diabetic neuropathy. Diabetes Res Clin Pract. 144, 177–191. doi: 10.1016/j.diabres.2018.08.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Snyder, M. J., Gibbs, L. M., and Lindsay, T. J. (2016). Treating painful diabetic peripheral neuropathy: an update. Am Fam Physician. 94, 227–234.

PubMed Abstract | Google Scholar

Todorovic, S. M. (2015). Is Diabetic Nerve Pain Caused by Dysregulated Ion Channels in Sensory Neurons? Diabetes. 64, 3987–3989. doi: 10.2337/dbi15-0006

PubMed Abstract | CrossRef Full Text | Google Scholar

Trapnell, C., Cacchiarelli, D., Grimsby, J., Pokharel, P., Li, S., Morse, M., et al. (2014). The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells. Nat Biotechnol. 32, 381–386. doi: 10.1038/nbt.2859

PubMed Abstract | CrossRef Full Text | Google Scholar

Urban-Ciecko, J., and Barth, A. L. (2016). Somatostatin-expressing neurons in cortical networks. Nat Rev Neurosci. 17, 401–409. doi: 10.1038/nrn.2016.53

PubMed Abstract | CrossRef Full Text | Google Scholar

Usoskin, D., Furlan, A., Islam, S., Abdo, H., Lönnerberg, P., Lou, D., et al. (2015). Unbiased classification of sensory neuron types by large-scale single-cell RNA sequencing. Nat Neurosci. 18, 145–153. doi: 10.1038/nn.3881

PubMed Abstract | CrossRef Full Text | Google Scholar

Vague, P., Dufayet, D., Coste, T., Moriscot, C., Jannot, M. F., and Raccah, D. (1997). Association of diabetic neuropathy with Na/K ATPase gene polymorphism. Diabetologia. 40, 506–511. doi: 10.1007/s001250050708

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, F., Cai, B., Li, K. C., Hu, X. Y., Lu, Y. J., Wang, Q., et al. (2015). FXYD2, a γ subunit of Na⁺, K⁺-ATPase, maintains persistent mechanical allodynia induced by inflammation. Cell Res. 25, 318–334. doi: 10.1038/cr.2015.12

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, K., Wang, S., Chen, Y., Wu, D., Hu, X., Lu, Y., et al. (2021). Single-cell transcriptomic analysis of somatosensory neurons uncovers temporal development of neuropathic pain. Cell Res. 31, 904–918. doi: 10.1038/s41422-021-00479-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Wingett, S. W., and Andrews, S. (2018). FastQ Screen: A tool for multi-genome mapping and quality control. F1000Res. 7, 1338. doi: 10.12688/f1000research.15931.2

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, Y. E., Pan, L., Zuo, Y., Li, X., and Hong, W. (2017). Detecting Activated Cell Populations Using Single-Cell RNA-Seq. Neuron. 96, 313-329.e316. doi: 10.1016/j.neuron.2017.09.026

PubMed Abstract | CrossRef Full Text | Google Scholar

Yamazaki, S., Yamaji, T., Murai, N., Yamamoto, H., Matsuda, T., Price, R. D., et al. (2013). FK1706, a novel non-immunosuppressive immunophilin ligand, modifies gene expression in the dorsal root ganglia during painful diabetic neuropathy. Neurol Res. 34, 469-477. doi: 10.1179/1743132812Y.0000000029

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, M., Wang, Y., Geng, J., Zhou, S., and Xiao, B. (2019). Mechanically Activated Piezo Channels Mediate Touch and Suppress Acute Mechanical Pain Response in Mice. Cell Rep. 26, 1419–1431.e1414. doi: 10.1016/j.celrep.2019.01.056

PubMed Abstract | CrossRef Full Text | Google Scholar