Han Zhou†
Xiaosheng Yang
Chenlong Liao
Fukai Ma
WenChuan Zhang- Department of Neurosurgery, Ninth People Hospital Affiliated to Shanghai Jiao Tong University School of Medicine, Shanghai, China
Mechanical allodynia (MA) is the main reason that patients with diabetic peripheral neuropathy (DPN) seek medical advice. It severely debilitates the quality of life. Investigating hyperglycemia-induced changes in neural transcription could provide fundamental insights into the complex pathogenesis of painful DPN (PDPN). Gene expression profiles of physiological dorsal root ganglia (DRG) have been studied. However, the transcriptomic changes in DRG neurons in PDPN remain largely unexplored. In this study, by single-cell RNA sequencing on dissociated rat DRG, we identified five physiological neuron types and a novel neuron type MAAC (Fxyd7+/Atp1b1+) in PDPN. The novel neuron type originated from peptidergic neuron cluster and was characterized by highly expressing genes related to neurofilament and cytoskeleton. Based on the inferred gene regulatory networks, we found that activated transcription factors Hobx7 and Larp1 in MAAC could enhance Atp1b1 expression. Moreover, we constructed the cellular communication network of MAAC and revealed its receptor-ligand pairs for transmitting signals with other cells. Our molecular investigation at single-cell resolution advances the understanding of the dynamic peripheral neuron changes and underlying molecular mechanisms during the development of PDPN.
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 GEXSCOPETM 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 GEXSCOPETM 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 GEXSCOPETM 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. 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. 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. 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. 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. 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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
Ingber, D. E. (2006). Cellular mechanotransduction: putting all the pieces together again. Faseb j. 20, 811–827. doi: 10.1096/fj.05-5424rev
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
Snyder, M. J., Gibbs, L. M., and Lindsay, T. J. (2016). Treating painful diabetic peripheral neuropathy: an update. Am Fam Physician. 94, 227–234.
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
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
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
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
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
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
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
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
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
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
Keywords: mechanical allodynia, neuropathic pain, painful diabetic peripheral neuropathy, somatosensory neurons, single-cell RNA sequencing
Citation: Zhou H, Yang X, Liao C, Chen H, Wu Y, Xie B, Ma F and Zhang W (2022) The Development of Mechanical Allodynia in Diabetic Rats Revealed by Single-Cell RNA-Seq. Front. Mol. Neurosci. 15:856299. doi: 10.3389/fnmol.2022.856299
Received: 17 January 2022; Accepted: 14 April 2022;
Published: 20 May 2022.
Edited by:
Fengxian Li, Southern Medical University, ChinaReviewed by:
Temugin Berta, University of Cincinnati, United StatesJunguk Hur, University of North Dakota, United States
Copyright © 2022 Zhou, Yang, Liao, Chen, Wu, Xie, Ma and Zhang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: WenChuan Zhang, zhangwench88@sjtu.edu.cn
†These authors share first authorship