Abstract

Background

   Neuropathic pain (NP), which results from injury or lesion of the
   somatosensory nervous system, is intimately associated with glial
   cells. The roles of microglia and astrocytes in NP have been broadly
   described, while studies on oligodendrocytes have largely focused on
   axonal myelination. The mechanisms of oligodendrocytes and their
   interactions with other glial cells in NP development remain uncertain.

Methods

   To explore the function of the interaction of the three glial cells and
   their interactions on myelin development in NP, we evaluated changes in
   NP and myelin morphology after a chronic constriction injury (CCI)
   model in mice, and used single-cell sequencing to reveal the
   subpopulations characteristics of oligodendrocytes, microglia, and
   astrocytes in the spinal cord tissues, as well as their relationship
   with myelin lesions; the proliferation and differentiation trajectories
   of oligodendrocyte subpopulations were also revealed using pseudotime
   cell trajectory and RNA velocity analysis. In addition, we identified
   chemokine ligand-receptor pairs between glial cells by cellular
   communication and verified them using immunofluorescence.

Results

   Our study showed that NP peaked on day 7 after CCI in mice, a time at
   which myelin lesions were present in both the spinal cord and sciatic
   nerve. Oligodendrocytes, microglia, and astrocytes subpopulations in
   spinal cord tissue were heterogeneous after CCI and all were involved
   in suppressing the process of immune defense and myelin production. In
   addition, the differentiation trajectory of oligodendrocytes involved a
   unidirectional lattice process of OPC-1-Oligo-9, which was arrested at
   the Oligo-2 stage under the influence of microglia and astrocytes. And
   the CADM1-CADM1, NRP1-VEGFA interactions between glial cells are
   enhanced after CCI and they had a key role in myelin lesions and
   demyelination.

Conclusions

   Our study reveals the close relationship between the differentiation
   block of oligodendrocytes after CCI and their interaction with
   microglia and astrocytes-mediated myelin lesions and NP. CADM1/CADM1
   and NRP-1/VEGFA may serve as potential therapeutic targets for use in
   the treatment of NP.

Supplementary Information

   The online version contains supplementary material available at
   10.1186/s12974-024-03207-3.

   Keywords: Intercellular communication, Myelin lesion, Neuroglia,
   Neuropathic pain, Single-cell RNA sequencing

Background

   Neuropathic pain (NP) is a chronic, persistent and intractable pain
   resulting from injuries or lesions of the somatosensory nervous system
   [[42]1]. The incidence of NP due to spinal stenosis and lumbar disc
   herniation increases with aging and large-scale shifts in lifestyle
   changes, reaching 18% in developing countries [[43]2]. Accordingly, NP
   substantially reduces the quality of life and increases the public
   health burden [[44]3].

   Glial cells play essential roles in the development of central
   sensitization. Activated microglia induce a sustained experience of NP
   by releasing inflammatory factors and activating inflammatory pathways,
   as well as by producing an imbalance in the M1/M2 polarization ratio
   after nerve injury [[45]4, [46]5]. Astrocytes regulate NP through
   changes in glial signaling pathways, receptors and channel protein
   expression and glia-derived factors [[47]6]. In addition, interactions
   mediated by IL-18/IL-18R between microglia and astrocytes have also
   been shown to promote NP [[48]7].

   Oligodendrocytes participate in NP progression by secreting
   neurotrophic factors that influence the survival of neuronal cells and
   astrocytes, as well as through their ability to influence to their own
   survival [[49]8]. However, current studies of oligodendrocytes have
   primarily focused on the mechanisms of central nervous system (CNS)
   axonal myelination [[50]9]. Microstructural changes in the myelin
   sheath contribute to aberrant neuronal firing by affecting nerve
   conduction and functions of neuronal circuits [[51]10, [52]11]. NP, as
   resulting from abnormal spontaneous neuronal firing and peripheral
   nerve demyelination, represent key factors in the development of
   diabetic NP [[53]12]. Consequently, we hypothesized that myelin lesions
   can promote NP initiation and progression, effects which may be
   associated with oligodendrocyte activity. However, whether
   oligodendrocytes and the generation of NP are associated with myelin
   lesions remain unclear. Moreover, potential relationships of
   interactions among oligodendrocytes, microglia and astrocytes, as
   related to NP, also require further investigation.

   In this study, we utilized a chronic constriction injury (CCI) model of
   the sciatic nerve. This model, which has been employed worldwide,
   better simulates the NP resulting from clinical peripheral nerve injury
   and has several advantages, including reduced invasiveness and greater
   reliability. Electron microscopy was used to observe any structural
   changes in the myelin sheath, while single-cell RNA sequencing
   (scRNA-seq) provided a means to compare gene expression levels in glial
   cells at a single-cell resolution using high-throughput methods
   [[54]13]. The heterogeneity and functional diversity among cell
   subpopulations were assessed, and results as obtained with
   receptor–ligand interactions served to reveal some of the mechanisms of
   interaction among oligodendrocyte, astrocyte and microglial
   subpopulations. Most notably, we unraveled the specific molecular
   mechanisms through which oligodendrocytes, astrocytes and microglia
   affect myelin sheaths as related to lesion-mediated NP development. The
   findings of this study offers crucial insights into the identification
   of safe and effective targets and interventions that can be applied in
   the treatment of NP.

Methods

CCI model and experimental protocols

   Within the spinal cord, male and female rodents process pain through
   different mechanisms, with microglia being primarily involved in males
   and T cells in females [[55]14]. In this report, male mice were
   selected as our model to study the effects of glial cell interactions
   on pain. Male C57 mice (18–20 g) were purchased from the Shandong
   University Animal Center and maintained with three mice per cage with
   food and water provided ad libitum. For induction of the sciatic nerve
   injury, mice were anesthetized with pentobarbital sodium via an
   intraperitoneal injection at 35 mg/kg and the surgical area was
   disinfected. The skin of the left femur area was incised and the muscle
   layers were bluntly dissected until the sciatic nerve was exposed.
   Without blocking the superficial vascular circulation, 6 − 0 silk was
   used to ligate the exposed sciatic nerve until a slight twitching of
   the ipsilateral hind limb was observed [[56]15]. After rinsing the
   wound with saline, the skin incision was sutured layer-by-layer. None
   of the mice died or were excluded from the experiment. The Animal Care
   and Use Committee at Shandong University of the Qilu Hospital reviewed
   and approved all experimental protocols (Approval number:
   KYLL-2020(KS)-363). Our study was performed following the Guide for the
   Care and Use of Laboratory Animals of the National Institutes of
   Health, USA.

Behavioral testing

   Mice were randomly divided into two groups, control (CON) and treatment
   (CCI), with 7 mice per group. In the CON group, the sciatic nerve was
   exposed but not ligated. The paw-withdrawal mechanical threshold (PWMT)
   and paw-withdrawal thermal latency (PWTL) were performed before (0) day
   and at 3, 7, 14, 21 days after CCI treatment using an electronic Von
   Frey apparatus (KW-CT-1, NJKEWBIO, Nanjing, China) and a thermal
   plantar test kit (KW-600, NJKEWBIO), respectively. Prior to testing,
   the mice were placed in the Plexiglas enclosures for a 30 min period to
   accustom them to the test environment. The Von Frey and thermal plantar
   tests were conducted consecutively on the same day. There was a 6–8 h
   interval between these two tests to prevent any potential for sensory
   thresholds to be affected. Each test was performed five times with
   a > 5 min interval between each of these five tests. The average score
   obtained from these five tests was then used in the statistical
   analyses. All procedures were performed by the same investigator.
   Differences between the two groups were compared using an unpaired
   Student’s t-test as achieved with use of the GraphPad Prism statistical
   software program (version 9.5.0). Data are presented as
   means ± standard error of the means (SEM). A P < 0.05 was required for
   results to be considered as statistically significant.

Transmission electron microscopy (TEM)

   TEM was used to assess any changes in the myelin sheath within the
   spinal cord and sciatic nerve. Three mice in each group (CON and CCI)
   were transcardially perfuse-fixed with 4% paraformaldehyde and 1.25%
   glutaraldehyde in phosphate buffer. The injured/sham and contralateral
   sides at the L4/L5 spinal cord level and the sciatic nerves at
   approximately 5 mm proximal to the sciatic nerve trifurcation were
   isolated and placed in the fixative solution for > 24 h. Sections were
   then sliced at 1 mm thickness. Targeted tissue blocks were trimmed from
   tissue samples and placed into a mixed fixative for fixation. After
   rinsing, osmic acid fixation, dehydration and electron staining, the
   resin was embedded, polymerized and cut into ultrathin sections. The
   sections were photographed while under electron microscopy (HT7700,
   Hitachi, Tokyo, Japan). Six intact myelin sheaths were randomly
   selected from each mouse and Image-pro plus 6.0 (Media Cybernetics,
   Inc., Rockville, MD, USA) was used to measure the diameters of
   myelinated nerve fibers and axons. The G-ratio (axon/myelinated fiber
   diameter) was calculated to evaluate the degree of demyelination.
   Differences between the two groups were analyzed using an unpaired
   Student’s t-test as described above.

Immunofluorescence (IF)

   IF analyses of neuropilin-1 (NRP-1), vascular endothelial growth factor
   A (VEGFA) and cell adhesion molecule 1 (CADM1) were performed within
   paraffin-embedded spinal cord tissue sections. Five sections from the
   L4/L5 spinal cord level of each group were randomly obtained for
   statistical analysis as described above. IF staining was performed by
   dewaxing the paraffin sections, antigen repair, sealing, primary
   antibody incubation, secondary antibody incubation and 4′,
   6-diamidino-2-phenylindole (DAPI) staining. Imaging was performed using
   a panoramic slide scanner (3DHISTECH, Budapest, Hungary) and the data
   from these images were analyzed using CaseViewer 2.4 (3DHISTECH). The
   primary antibodies included anti-NRP-1 (1:100; cat. no. ET1609-69;
   HUABIO, Hangzhou, China), anti-VEGF (1:100; cat. no. ET1604-28; HUABIO)
   and anti-CADM1 (1:100; cat. no. DF6679; Affinity, Cincinnati, OH, USA).
   The secondary antibody was goat anti-rabbit antibody (1:1,000;
   Proteintech, Chicago, IL, USA).

Quality control (QC) and preprocessing of scRNA-seq data

   On day 7 post-injury/sham treatment, L4/L5 spinal cord section samples
   were obtained from three randomly selected mice of each group for
   scRNA-seq assay. Cells were loaded onto the Chromium single-cell
   controller (10× Genomics) to generate single-cell gel beads in the
   emulsion according to the manufacturer’s protocol. ScRNA-seq libraries
   were constructed using Single Cell 3′ Library and Gel Bead Kit version
   3.1 and sequenced at an average-read target depth of 50,000 reads per
   cell from total gene expression libraries using the NovaSeq 6000
   sequencer (Illumina). The number of genes per cell are listed in
   Supplementary Table [57]1.

   A raw gene-expression matrix was generated using Cell Ranger (version
   5.0.0; 10× Genomics, Pleasanton, CA, USA), with the mm10 (version
   2020-A) mouse genome as a reference and the other parameters set at the
   default value. The gene-cell unique molecular identifier (UMI) matrix
   was analyzed using R software (version 3.6.1; The R Foundation for
   Statistical Analysis, Vienna, Austria) and the Seurat package (version
   3.2.0). To perform the QC, genes expressed in fewer than three cells
   were filtered out and those cells with > 200 genes were used for
   further analysis. Mitochondrial genes were expressed at a proportion of
   < 10%. Doublets were defined by Scrublet (version 0.2.2), and cells
   with a doublet score > 0.2 were removed from the downstream analysis.
   After performing the QC, the dataset comprised 39,711 cells with 23,567
   genes. The UMI count matrix was normalized using the NormalizeData
   function. For this calculation, the number of UMIs in each gene was
   divided by the total number of UMIs in each cell, multiplied by 10,000
   and then converted to a log scale.

   We selected the 2,000 genes with the highest standardized intercellular
   variation (HVGs) to reduce noise using the FindVariableFeatures
   function with use of the vst selection method [[58]16]. The percent of
   mitochondrial and heat-shock protein UMI counts were regressed. The
   normalized data were scaled to z-scores and the percent of
   mitochondrial reads was regressed using the ScaleData function,
   followed by the principal component analysis (PCA) with the RunPCA
   function. An unsupervised graph-based clustering algorithm was used to
   determine cell heterogeneity. The ElbowPlot and DimHeatmap functions
   were used to identify the appropriate PCA dimensions.

   We applied the FindNeighbors and FindClusters functions with a
   resolution of 1 to perform a first-round clustering and annotated each
   cluster as based on canonical marker genes. Small clusters expressing
   the dual-lineage genes, Csf1r-Mog, were defined as doublets. Erythroids
   and doublets were removed from the downstream analysis. We identified
   16 major subsets, including four immune cell types, 10 nonimmune cell
   types and two cycling cell types. A second-round clustering for each
   cell type was performed to identify the functional cellular subsets
   within these cell types. The resolution selection was based on the data
   characteristics, with the top n principal components (PCs) used to
   identify clusters. For the clustering of oligodendrocyte precursor cell
   (OPC) and OligoDendrocytes (Oligos), the top 12 PCs were selected, for
   microglia the top six PCs were selected, for macrophages the top 15 PCs
   were selected and for astrocytes the top six PCs were selected. For the
   clustering of T cells, PCs with a substantial contribution from CD4 and
   CD8A were prioritized for dimensional reduction and subclustering.
   Clustering results were visualized by uniform manifold approximation
   and projection (UMAP) using the RunUMAP function with default settings.

Assessment of associations between cell clusters within the CON and CCI
groups

   A mixed-effects modeling of associations of single cells (MASC, version
   0.1.0-alpha) algorithm was applied to account for random effects when
   testing proportional differences within each cell cluster between the
   CON and CCI groups [[59]17]. The donor was specified as a random effect
   covariate. A detailed description of statistics for the MASC of the 16
   major subsets and 45 subpopulations (34 subgroups with secondary
   clustering and 11 major subsets without secondary clustering) is
   contained in Supplementary Table [60]2. Events with adjusted P-values
   of < 0.05, as based on the Benjamini–Hochberg (BH) procedure, were
   considered as statistically significant.

Identification of signature genes

   Differentially expressed genes (DEGs) were identified using the limma
   package (version 3.42.2). Briefly, DEGs within one cell cluster were
   determined by comparing the cells of its cluster with those of other
   clusters. The donor effect was considered and added to the end of the
   call to model the matrix. Only genes with a minimum log[2] fold change
   (FC) of 0.25 and maximum adjusted P-value (via the Bonferroni method)
   of 0.05 were considered as DEGs for the cell cluster (Supplementary
   Table [61]3). The Wilcoxon rank-sum test was applied to determine
   differences in DEGs between the CON and CCI groups. Only DEGs with a
   minimum percent of 0.1 in expressing cells, a minimum log[2] FC of 0.2
   and a maximum adjusted P-value (via the Bonferroni method) of 0.05 were
   considered as DEGs [[62]18].

Gene set enrichment analysis

   Pathway enrichment analysis was performed using the R package fgsea
   (version 1.17.1). DEGs for enrichment between the groups were analyzed
   using the Wilcoxon rank-sum test, with only genes demonstrating a
   minimum proportion of expressing cells of 0.1 being retained. For
   cell-type comparisons, the DEGs for enrichment were calculated using
   limma with an FDR of < 0.05. For the Gene Ontology (GO), gene sets of
   biological processes (BPs) from MsigDB (version 7.1) were used for GO
   enrichment. GO terms with BH-adjusted P-values of < 0.05 were
   illustrated in a heatmap using the ggplot2 R package (version 3.3.3).
   For the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway
   enrichment analysis, canonical pathway gene sets from MsigDB (version
   7.1) were used, with pathways showing BH-adjusted P-values < 0.05 being
   retained.

Developmental trajectory inference

   To construct the OPC-oligodendrocyte cell-state transition, Monocle3
   (version 1.0.0) was applied to the expression matrix, and UMAP was
   embedded to preserve the local relations of cells and pseudotime from
   beginning to end [[63]19]. Briefly, the principal graph was displayed
   in UMAP using the learn_graph function with a minimal_branch_len = 15,
   indicating differentiation trajectories. The principal graph node
   containing the highest fraction of OPC-1 was assigned as the root, with
   the pseudotime then calculated using the order cell function.

   RNA velocity analysis was performed on the OPC-Oligo data to confirm
   trajectory results from Monocle3. The spliced and unspliced UMIs of
   each gene within each cell were counted using the Python package
   Velocyto (version 0.17.17). Subsequent analyses were performed using
   ScanPy (version 1.7.2) and scVelo [[64]20]. The count matrices were
   first normalized by library size and then filtered to retain only those
   genes that could be commonly detected in more than 30 cells for both
   spliced and unspliced matrices and those genes exhibiting high
   variability. A k-nearest-neighbor graph (k = 30) was constructed using
   the top 30 PCs. For each cell, the moments (means and uncentered
   variances) of the normalized spliced/unspliced counts were computed
   using the 30 nearest neighbors with the function scv.pp.moments. The
   moments facilitated the RNA velocity estimation implemented in function
   scv.tl.velocity with the mode set to “dynamical” and the percent set to
   5–95%. Estimated velocities were used to construct a velocity graph
   representing the transition probabilities among cells using the
   scvelo.tl.velocity_graph function with default parameters. Finally, a
   velocity graph was used to embed the RNA velocities into the UMAP or
   diffusion map on a grid using the scv.pl.velocity_embedding_grid
   function with default parameters.

Definitions of macrophage and microglial cluster phenotypes

   Phenotypic scores were defined as the average expression levels of the
   signature genes. The signature genes used to define M1/M2 phenotype
   scores of macrophage and microglial clusters are listed in
   Supplementary Table [65]4 [[66]21].

Similarity analysis of clusters from different datasets

   To compare similarities between the OPC and Oligo clusters in the
   different studies, the data from our study were aligned with those from
   a previously published nervous system study using the integration
   method in Seurat. In total, 16,200 cells from our dataset and 31,062
   cells from the previously published study were included [[67]22]. We
   searched the anchors between the two datasets in the top 30 canonical
   correlate spaces using the FindIntergrationAnchors function and
   integrated the data using the IntegrateData function with default
   parameters. Correlations were employed to assess the relationship
   between cell subtypes for our OPC-Oligos with those from the published
   study using an integrated matrix of 2,000 highly variable genes. The
   correlation matrix of the corresponding cell subtypes was used for
   heatmap visualization using the pheatmap function in the pheatmap R
   package (version 1.0.12).

Identification of significant chemokine receptor–ligand pairs

   Chemokine receptor–ligand pairs were evaluated using CellPhoneDB
   (version 2.1.1) as based on a previously described curated chemokine
   receptor–ligand interaction map [[68]23]. The linking line displayed
   only the interaction pairs with permutation-test P-values < 0.05. The
   gene-pair-related pathways were annotated by CellchatDB (version
   1.1.3), and comparisons between the CON and CCI groups were performed
   using the rankNet function with default parameters.

Results

NP and ultrastructural changes within the myelin sheath of the spinal cord
and sciatic nerve after CCI

   To evaluate the impact on NP following CCI, behavioral tests were
   performed before (0) day and at 3, 7, 14, 21 days after CCI treatment.
   Before testing, the mice were placed in the Plexiglas enclosures for
   30 min to acclimatize them to the test environment. There was a minimum
   of a 5 min interval between tests. On day 7, mice were euthanized and
   spinal cord and sciatic nerve tissue samples were extracted for
   evaluation using scRNA-seq, IF and TEM. The specific processes are
   illustrated in Fig. [69]1A.

Fig. 1.

   [70]Fig. 1
   [71]Open in a new tab

   Effects of CCI on NP and myelin ultrastructure. (A) Protocol workflow
   illustrating the temporal sequences for the CCI model, behavioral
   tests, scRNA-seq, IF and TEM. Time-varying curves of (B) PWMT and (C)
   PWTL in mice. ####, P < 0.0001, CCI group vs. CON group, n = 7 per
   group. (D) TEM of spinal cord samples from CON and CCI groups, n = 3
   per group. (E) Bar graph showing quantitative analysis of g-ratios from
   spinal cord samples of CON and CCI groups. Data are presented as
   means ± SEMs. #, P < 0.05. (F) TEM of sciatic nerve samples from CON
   and CCI groups, n = 3 per group. (G) Bar graph showing quantitative
   analysis of g-ratios from sciatic nerves of CON and CCI groups. Data
   are presented as means ± SEMs. ##, P < 0.01

   There were no statistically significant differences in pre-operative
   pain thresholds between the two groups (Fig. [72]1B-C, P > 0.05) or
   over any of the time periods sampled in the CON group (Fig. [73]1B-C,
   P > 0.05). On postoperative day 3, the mechanical and thermal pain
   thresholds in the CCI group were significantly decreased as compared
   with those of the CON group (Fig. [74]1B-C, P < 0.0001). On
   postoperative day 7, mechanical and thermal pain thresholds decreased
   to their lowest levels in the CCI group as compared with those in the
   CON group (Fig. [75]1B-C, P < 0.0001). On postoperative day 14,
   mechanical and thermal pain thresholds remained at significantly lower
   levels in the CCI versus CON group (Fig. [76]1B-C, P < 0.0001),
   although the CCI group showed a slight rebound. On postoperative day
   21, the overall mechanical and thermal pain thresholds in the CCI group
   remained significantly decreased as compared with those in the CON
   group (Fig. [77]1B-C, P < 0.0001). As a full neuropathic phenotype
   (lowest pain thresholds and most pronounced hyperalgesia) was present
   on day 7 post-CCI and there is a report indicating that a full
   neuropathic phenotype develops at this interval [[78]24], this time
   point was selected for subsequent analyses.

   To determine the extent of myelin sheath damage after CCI, TEM was used
   to compare myelin morphology in the spinal cord and sciatic nerve
   within the two groups, with the unpaired Student’s t-test used for
   analysis of the g-ratio (Fig. [79]1D-G). Not surprisingly, the degree
   of nerve injury was less severe in the CON group. In this CON group,
   the myelin sheath structure was complete and dense, the boundary was
   clear and there was no demyelination. In contrast, myelin sheath
   thickness in the CCI group was significantly reduced and lamellar
   myelin sheaths were sparsely arranged and disordered. The periaxonal
   space with a large gap between the axolemma and the myelin sheath, and
   g-ratios were significantly increased (P < 0.05) and demyelination,
   along with incomplete remyelination, was present in the spinal nerves
   of the CCI group (Fig. [80]1D-E).

   The degree of sciatic nerve demyelination after CCI was also more
   severe. As compared with the CON group, the g-ratio of the CCI group
   was significantly increased (P < 0.01). Notably, the sciatic nerve
   myelin sheath within the CON group was sparse and disordered, with the
   degree of nerve injury being greater than that observed within the
   spinal cord in this group (Fig. [81]1F-G). This disruption may be
   related to the unintentional induction of small injuries resulting from
   procedures involved with exposing (but not ligating) the sciatic nerve
   in the CON group.

Heterogeneous expression of 16 major cell types in mouse spinal cord tissues
after CCI and participation of glial cells in immune defense and myelin
sheath lesions

   We next attempted to determine the basis for the demyelinated spinal
   cord lesions in response to CCI as assessed in six samples of spinal
   cord tissues from CON and CCI mice. A scRNA-seq assay was employed to
   analyze the heterogeneity of spinal cord cells as isolated from the CON
   versus CCI group. After completing the QC, we performed UMAP
   dimensionality reduction and clustering across 39,711 cells divided
   into 33 subpopulations (Fig. [82]2A). As grounded on results from a
   previous study [[83]25], the original subcluster was merged into 16
   main subsets based on the most-significant DEGs (Fig. [84]2B), which
   consisted of four immune cell types (T cells, B cells, macrophages, and
   neutrophils), 10 nonimmune cell types (OPCs, oligodendrocytes,
   microglia, astrocytes, fibroblasts, Schwann cells, pericytes,
   endothelial cells, ependymal cells and neurons) and two cycling cell
   types (proliferative OPCs and proliferative cells) (Fig. [85]2C). The
   DEGs and their distributions were obtained for the OPCs expressing
   Pdgfra and C1ql1; proliferative OPCs expressing Pdgfra, C1ql1 and
   Mki67; and oligodendrocytes expressing Mobp, Mog, Mag, Mbp and Plp1
   (Fig. [86]2D). Notably, cells identified as OPCs, oligodendrocytes and
   microglia belonged to multiple clusters, suggesting that these cell
   types may share a heterogeneity (Fig. [87]2C).

Fig. 2.

   [88]Fig. 2
   [89]Open in a new tab

   ScRNA-seq reveals changes in mouse spinal cord cell composition after
   CCI. (A) UMAP visualization of 39,711 cells within 33 subpopulations in
   mouse spinal cord tissue. Each point represents a single cell, colored
   by cell clusters. Subpopulations with low-quality or doublet scores
   were removed. (B) Bubble plots of DEGs in different cell clusters. Dot
   size represents the percentage of cells expressing DEGs. Color
   represents the average scaled RNA expression of that gene in the
   cluster. (C) UMAP revealing 16 main cell subsets in mouse spinal cord
   tissue. (D) UMAP plots showing expression of marker genes by
   oligodendrocytes, OPCs and proliferative OPCs. Red bars indicate
   expression levels of corresponding genes. (E) Bar plots showing percent
   of each cell type in six spinal cord samples. Each stacked bar
   represents a cell cluster, with the total number of cells being scaled
   to 1. (F) Cell counts for different cell types in six spinal cord
   samples. (G) Heat map showing differential proportions of cell types
   between the CON and CCI groups. P-values were determined using MASC. *,
   P < 0.05, CCI group vs. CON group. (H) Bar graph showing percent of
   information flow from different signaling pathways in the CON and CCI
   groups. Row names indicate pathway identities. Blue indicates
   statistical significance in the CON group and red in the CCI group, as
   determined using CellChat

   Heterogeneous cell expression was observed in the CON and CCI groups as
   based on the composition ratio and cell numbers of the various cell
   types from different sources (Fig. [90]2E-F). Among these, due to the
   influence of factors such as sampling site and tissue source, the
   proportions and numbers of spinal cord cells vary as reported in
   different studies [[91]25, [92]26]. In this study, myelination-related
   cells (Schwann cells, pericytes, endothelial cells and ependymal cells)
   were more frequently expressed in the CON versus CCI group
   (Fig. [93]2E-G). Glial cells (OPCs, proliferative OPCs,
   oligodendrocytes, microglia and astrocytes) and immune cells
   (macrophages and T cells) showed higher levels of expression in the CCI
   versus CON group (Fig. [94]2E). In addition, due to their large size
   and irregular shapes, neurons rarely show a preserve intact morphology
   in cytoplasmic sequencing, thus resulting in a lower number
   (Fig. [95]2F).

   Results obtained from our pathway-interaction analysis within the CON
   and CCI groups revealed a heterogeneity in pathway expressions. There
   was a statistically significant increase in the CCI group for levels of
   intercellular adhesion molecule (ICAM). ICAM-1, a ligand for ICAM,
   participates in inflammatory processes and host-defense mechanisms
   mediated by T cells. Midkine, insulin-like growth factor, tumor
   necrosis factor superfamily member 14, large T antigen, somatostatin,
   hypocretin neuropeptide precursor, pyroglutamylated RFamide peptide,
   pituitary adenylate cyclase-activating polypeptide, nerve growth
   factor, thyrotropin-releasing hormone and major histocompatibility
   complex-1 were also statistically significant (Fig. [96]2H).

   These results suggest that the myelin sheath lesions and inflammatory
   immunity that occurs in mouse spinal cord tissues in response to CCI
   may be related to glial cell activation.

Heterogeneity of 12 subpopulations of spinal cord oligodendrocytes in mice
after CCI and expressions of marker genes associated with demyelination

   The findings as presented above suggest that oligodendrocytes were the
   most abundant cell type in spinal cord tissue. Therefore, in our next
   series of experiments we assessed the heterogeneity and functional
   changes in oligodendrocytes following CCI as a means to determine their
   roles in demyelination (Fig. [97]2A).

   We identified 12 oligodendrocyte subpopulations as based on UMAP
   visualizations (Fig. [98]3A) and used multiple markers to reveal
   subpopulation characteristics based on heat mapping of the top 20 DEGs
   in each subpopulation (Fig. [99]3B-C). Results were expressed as
   follows: OPC-1 (Cspg5, Ptprz1, Pdgfra, C1ql1), OPC-2 (Fyn, Gpr17, Cd9,
   Tnr), proliferative OPCs (Hmgb2), Oligo-1 (Plekha1, Tspan2, Ptgds),
   Oligo-2 (Glul, Apoe, Car2), Oligo-3 (Il33, Opalin, Qdpr), Oligo-4
   (Serpinb1a, Anln, Trf), Oligo-5 (Plin3), Oligo-6 (Enpp6, Klk6), Oligo-7
   (Hopx, Mgst3), Oligo-8 (Cdkn1c, Anxa5) and Oligo-9 (Uchl1, S100b).

Fig. 3.

   [100]Fig. 3
   [101]Open in a new tab

   ScRNA-seq reveals heterogeneous expression of spinal cord
   oligodendrocyte subpopulations in mice after CCI. (A) UMAP plots
   showing oligodendrocyte subpopulations in second-run clustering
   analyses. (B) UMAP plots showing levels of specific marker genes in
   each subpopulation of oligodendrocytes. (C) Heat map showing mean
   levels of the top-20 DEGs in each subpopulation. Each column
   corresponds to a subpopulation and each row corresponds to a gene. (D)
   UMAP plots showing expression distributions of oligodendrocytes in six
   spinal cord samples. (E) Bar graph showing expression of
   oligodendrocyte subpopulations in CCI and CON groups presented as the
   percent of total number of cells within each subpopulation. Stacked
   bars indicate cell frequencies of subpopulations in the CCI and CON
   groups. (F) Bar graph showing the number of different oligodendrocyte
   subpopulations in the CON and CCI groups. (G) Heat map showing
   differential proportions of cell subsets in the CON and CCI groups. *,
   P < 0.05; ***, P < 0.001; CCI group vs. CON group, as determined using
   MASC. (H-J) Heat map showing DEGs among proliferative OPCs, OPCs and
   oligodendrocytes in the CON and CCI groups. (K-M) Bar graphs showing GO
   enrichment pathway annotations corresponding to DEGs between the CON
   and CCI groups as presented in H-J. (N, O) UMAP plots showing
   expression distributions of K-M differentially expressed marker genes
   and their expression levels in the CON and CCI groups

   We found that oligodendrocytes were widely expressed in different
   samples (Fig. [102]3D), but expression ratios and cell numbers within
   each subpopulation varied as a function of the stage examined
   (Fig. [103]3E-F). While all subpopulations were identified in both the
   CON and CCI groups, Oligo-5, -7, -8, and − 9 were mainly expressed in
   the CON group (Fig. [104]3E). When analyzing the differential
   expressions within each subpopulation, we found that significant
   differences were present in the expressions of OPC-1, proliferative
   OPCs, Oligo-7, Oligo-2 and Oligo-9 between the CON and CCI groups
   (Fig. [105]3G). Of these, the expressions of proliferative OPCs, OPC-1
   and Oligo-2 were significantly increased in CCI groups.

   Comparisons of DEGs in CCI- versus CON-derived cells were performed as
   an approach to dissect the roles of OPCs, proliferative OPCs and
   oligodendrocytes in spinal cord demyelination (Fig. [106]3H-J). Our GO
   enrichment analysis showed that the genes upregulated in CCI-derived
   OPCs, proliferative OPCs and oligodendrocytes were significantly
   enriched in genetic information signaling, cellular metabolism and
   immune defense (Fig. [107]3K-M). The differential marker genes Opalin,
   Ptprz1, Klk6 and Pmp22 for each pathway were analyzed as shown in
   Fig. [108]3N-O.

   Taken together, these findings reveal that CCI promoted the
   differentiation and heterogeneous expression of spinal oligodendrocyte
   subpopulations and that oligodendrocytes may play a role in the
   demyelination of spinal cord neurons to then mediate pain sensation.

Heterogeneity of nine subpopulations of spinal microglia after CCI and their
involvement in changes of immune defense and myelin structure

   Given the abundance of microglia in CCI samples (Fig. [109]2A, C),
   their roles in promoting oligodendrocyte differentiation after spinal
   cord damage and the necessity for remyelination following
   demyelination, we next directed our efforts toward investigating the
   intrinsic structure and function of the entire microglial population
   [[110]27]. Cluster analysis was used to separate the microglia into
   nine subgroups (Microglia-1 to Microglia-9), each demonstrating unique
   characteristics (Fig. [111]4A-B). The subgroups Microglia-1, -2 and − 7
   were mainly enriched in the CCI group, whereas Microglia-3, -4, -5 and
   − 8 were expressed at higher levels in the CON group (Fig. [112]4C-D),
   and their differential expressions were significantly different between
   subpopulations (Fig. [113]4E). Results from UMAP analysis showed that
   microglia within different stages formed distinct heterogeneous
   clusters (Fig. [114]4F-G). However, findings from the heat map of the
   top 10 DEGs also revealed considerable similarities among
   subpopulations (Fig. [115]4B).

Fig. 4.

   [116]Fig. 4
   [117]Open in a new tab

   ScRNA-seq reveals transcriptional differences in mouse spinal microglia
   after CCI. (A) UMAP plots showing microglial subpopulations as obtained
   using subclustering analyses. (B) Heat map showing levels of the top-10
   DEGs in each subpopulation, with each column and row representing a
   single subpopulation and gene, respectively. (C) Bar graph showing
   expressions of microglial subpopulations in the CON and CCI groups,
   presented as the percent of total number of cells in this
   subpopulation. Stacked bars indicate frequencies of cell subpopulations
   in the CON and CCI groups. (D) Bar graph showing cell counts for
   different microglial subpopulations in the CON and CCI groups. (E) Heat
   map showing differential proportions of microglial subsets in the CON
   and CCI groups. *, P < 0.05; **, P < 0.01; ***, P < 0.001; CCI group
   vs. CON group. (F) UMAP plots showing expression distributions of
   microglia in six spinal cord samples. (G) UMAP plots showing expression
   distributions of microglia in the CON and CCI groups. (H, I) Bars
   showing enrichment of top-20 GO terms and KEGG pathways in CCI
   microglia. (J) Heat map showing DEGs in microglia within the CON and
   CCI groups. (K) Bar graph showing GO annotations of DEGs in microglia
   within the CON and CCI groups

   To verify the presence of myelin sheath-associated marker genes among
   DEGs in each subpopulation of microglia, GO and KEGG pathway enrichment
   analyses were performed for the microglial subpopulations. According to
   results of the GO analysis, microglial genes were found to be enriched
   in cell metabolism, immune defense and the transmission of genetic
   information after CCI (Fig. [118]4H). Results of the KEGG pathway
   analysis showed enrichment in ribosome production, pro-inflammatory
   factor production induction, immune activation, cell metabolism and
   signaling, revealing that microglia may participate in inflammatory
   immunity after CCI, and correlate with changes in myelin(Fig. [119]4I).
   When further analyzing the genes in each pathway we found that Plp1 was
   frequently expressed in the top 20 GO and KEGG pathways.

   To assess potential differences in expression status between microglia
   originating from CON versus CCI mice, DEGs and GO-related enrichment
   pathways were analyzed in microglia. We detected 16 upregulated and 10
   downregulated genes which were significantly expressed (Fig. [120]4J).
   Pathway expressions of immune-promoting factors and astrocyte
   development were found to be significantly downregulated in CCI mice
   (Fig. [121]4K).

   Collectively, these results suggest that the expression of microglial
   subpopulations in spinal cord tissue was heterogeneous and that
   microglia were associated with a reduction in immune defenses and
   myelin structure changes after CCI.

Subpopulations of spinal cord astrocytes and inhibition of immune defense and
myelin sheath synthesis by astrocytes following CCI

   During the neuroinflammation resulting from CCI, astrocytes function as
   heterogeneous agents. Therefore, we performed a clustering analysis and
   visualization using UMAP plots to evaluate the changes in astrocytes
   that occur in response to CCI [[122]28]. From these analyses we found
   two distinct subpopulations of astrocytes, Astrocytes-1 and
   Astrocytes-2 (Fig. [123]5A). The top 20 specific markers for each
   subpopulation are shown in Fig. [124]5B. In addition, when
   investigating the expression of astrocytes in multiple samples and
   states, as based on specific subpopulation markers (Fig. [125]5C-F), we
   observed that cell numbers within the astrocyte subpopulations
   increased after CCI (Fig. [126]5D, F). While this increase was greater
   for Astrocytes-1 (Fig. [127]5E-F), differences between Astrocytes-1 and
   − 2 failed to achieve statistical significance with regard to astrocyte
   subpopulations (Fig. [128]5G, P = 0.259).

Fig. 5.

   [129]Fig. 5
   [130]Open in a new tab

   ScRNA-seq reveals that highly heterogeneous subpopulations of
   astrocytes exist in mouse spinal cords after CCI. (A) UMAP plots
   showing subpopulations of astrocytes as resulting from subclustering
   analyses. (B) Heat map showing levels of DEGs in each subpopulation.
   Each column and row represent individual subpopulations and gene
   expressions, respectively. (C) UMAP plots showing distributions of
   astrocyte expression in six spinal cord samples. (D) UMAP plots showing
   distributions of astrocyte expression in the CON and CCI groups. (E)
   Bar graph showing expressions of astrocyte subpopulations in the CON
   and CCI groups as a percent of total number of cells within this
   subpopulation. Stacked bars indicate cell frequencies of subpopulations
   in the CON and CCI groups. (F) Bar graph showing cell counts for
   different astrocyte subpopulations in the CON and CCI groups. (G) Heat
   map showing differential expressions of astrocyte subpopulations in the
   CON and CCI groups. (H) Heat map showing significantly up- and
   down-regulated genes in astrocytes in the CON and CCI groups. (I) Bar
   graph showing GO annotations for DEGs in astrocytes between the CON and
   CCI groups

   DEGs in astrocytes derived from the CON and CCI groups were analyzed as
   a means to assess the function of astrocytes in the CCI-induced
   demyelination of the spinal cord. We identified 14 downregulated and 16
   upregulated genes (Fig. [131]5H). Using these DEGs, we further analyzed
   astrocyte functional variations in various states and found that
   pathways linked to immune cell activation and axonal myelin production
   were downregulated after CCI (Fig. [132]5I).

   These findings demonstrate that in response to CCI a heterogeneity in
   function exists between different astrocyte subpopulations, though
   their expression levels remain similar. Such findings suggest that
   interactions between astrocytes and other cells may result in torsion
   of their functions after CCI.

Pseudotime analysis revealed that oligodendrocyte subpopulations are involved
in myelin sheath formation and the blocking of this process after CCI

   As CCI induced a heterogeneity in the expression of individual
   oligodendrocyte subpopulations within the spinal cord, we next
   investigated the transformation of each oligodendrocyte subpopulation
   by reconstructing the pseudotemporal trajectory of oligodendrocyte
   differentiation in the spinal cord using a Monocle3 pseudotime analysis
   based on the clustering of each subpopulation (Fig. [133]6A). Our
   findings indicated that oligodendrocyte differentiation started with
   OPC-1 cells and progressed through a unidirectional lattice process of
   pathways involving OPC-2 to Oligo1-Oligo9 (Fig. [134]6A), results which
   were corroborated from the proposed time series as shown in
   Fig. [135]6B. Therefore, we evaluated the DEGs showing significant
   variation (Fig. [136]6C) and found that Cspg4 [[137]29], Gria2
   [[138]30], Gnai1 [[139]31], Sh3kbp1 [[140]32], Tanc2 [[141]33], Gal3st1
   [[142]34], and Arntl2 [[143]35] all demonstrated increased levels of
   gene expression.

Fig. 6.

   [144]Fig. 6
   [145]Open in a new tab

   Proposed pseudotime analysis of oligodendrocyte differentiation and
   mitosis in mouse spinal cord tissues. UMAP plots of pseudotime
   trajectories of oligodendrocyte differentiation in mouse spinal cord
   tissues as performed using (A) Monocle3 and (B) scVelo. (C)
   Pseudotemporal heat map showing gene expression dynamics during
   oligodendrocyte differentiation. Rows represent genes and columns
   pseudotime ordering of oligodendrocyte differentiation. (D) UMAP plots
   showing expression patterns of oligodendrocytes in mouse spinal cord
   tissues at different stages of the cell cycle: G1 - G1 phase, G2/M -
   G2/M phase and S - S phase. (E) UMAP plot showing distribution of
   proliferative OPC expression in the S and G2/M phases. (F) UMAP plots
   showing average expression of proliferative OPCs in the S or G2/M
   phases. These cell cycle-related genes were from the Seurat R package.
   (G) Diffusion plot showing pseudotemporal trajectory of OPC-1 and
   proliferative OPCs along cell cycle after CCI. (H) Pseudotemporal heat
   map showing gene expression dynamics during mitotic process of OPC-1
   and proliferative OPCs. Rows represent genes and columns pseudotemporal
   ordering of mitotic process. (I) RNA velocity analysis of expression
   patterns of representative marker genes in the G1, G2/M, and S phases.
   (J, K) Diffusion plots showing pseudotemporal trajectories of
   oligodendrocyte differentiation in the (J) CON and (K) CCI groups. (L,
   M) Heat map showing expression of DEGs (rows) and pseudotemporal
   oligodendrocyte differentiation (columns) in the (L) CON and (M) CCI
   groups

   When examining the status of oligodendrocytes in the cell proliferation
   cycle, we found that there was a heterogeneous expression of
   oligodendrocytes as a function of the various cell-cycle stages, as
   demonstrated via clustering analysis and visualization using UMAP plots
   (Fig. [146]6D). Notably, proliferative OPCs tended to be more
   frequently expressed during the S and G2/M phases (Fig. [147]6E-F).
   When separating OPC-1 from the proliferative OPCs, primarily expressed
   during the cytokinesis cycle, for pseudotime analysis (Fig. [148]6G),
   we observed that oligodendrocytes exhibited a normal cytokinesis cycle
   and were capable of normal proliferation.

   We further analyzed the DEGs showing varying expression patterns in
   cytokinesis according to pseudotime (Fig. [149]6H) and conducted an RNA
   velocity analysis on marker genes to characterize the G1, S and G2/M
   phases (Fig. [150]6I). With this analysis Zbtb20 [[151]36], Ezh2
   [[152]37], Brca1 [[153]38], Hdgf [[154]39], Pex5l [[155]40], Stmn1
   [[156]41], and Vegfa [[157]42] were all found to be highly expressed in
   the cell cycle.

   Next, the differentiation characteristics of each subpopulation of
   oligodendrocytes in the CON and CCI groups were analyzed to identify
   any potential changes in their differentiation trajectory after CCI. As
   shown in the diffusion plots of Fig. [158]6J-K, a heterogeneous
   differentiation of oligodendrocytes was present between the CON and CCI
   groups. Through pseudotemporal analysis, we further analyzed changes in
   DEGs of the two groups and found that Serpine2 [[159]43], Ntrk2
   [[160]44], Fads2 [[161]45], Kcnmb4 [[162]46], Rgcc [[163]47], Slc1a1
   [[164]48], Dusp10 [[165]49] and Car2 [[166]50] were all significantly
   expressed in the CON group (Fig. [167]6L). In the CCI group, Cntn4,
   Adgrb3, Rgs7, Sox6, Slc1a1, Ttll7 and Car2 were significantly expressed
   (Fig. [168]6M). In contrast to that observed in the CON group, the
   expression of Car2 gradually decreased in the Oligo-2 stage in the CCI
   group, while the expression of Sox6 was increased.

   Overall, these findings suggest that oligodendrocyte proliferation and
   differentiation were both required for myelin sheath synthesis and
   axonal signaling in spinal cord tissues, and CCI inhibited this
   process.

ScRNA-seq reveals that developmental oligodendrocyte differentiation is
blocked in the Oligo-2 stage after CCI

   As the differentiation process of oligodendrocytes in the CON and CCI
   groups was shown to be heterogeneous (Fig. [169]6J-K), we next
   evaluated the differentiation process of the oligodendrocyte
   subpopulations in the CON and CCI groups using a clustering analysis
   (Fig. [170]7A-B). The visualization of results for the different
   samples and sources assessed is shown in Fig. [171]7C-D. Significant
   differences were observed in the expression of proliferative OPCs,
   OPC-1 and Oligo-2 between the CCI and CON groups, with Oligo-2 showing
   the largest difference (Fig. [172]7E-F). Importantly, during the
   differentiation of each subgroup, the expression of each subtype in the
   CON group, after the Oligo-2 phase, gradually increased whereas that of
   each subgroup in the CCI group from differentiation to the Oligo-2
   stage showed a downward trend. Such findings suggest that Oligo-2 acted
   as a node in blocking the differentiation of oligodendrocytes after CCI
   (Fig. [173]7G).

Fig. 7.

   [174]Fig. 7
   [175]Open in a new tab

   ScRNA-seq during differentiation of each subpopulation of
   oligodendrocytes in the CON and CCI groups. UMAP plots showing
   expression distributions of oligodendrocytes in the (A) CON and (B) CCI
   groups. Cells were grouped into 12 subgroups by clustering, with
   subgroups being marked with different colors in B. (C) UMAP plots
   showing expression patterns within each subpopulation of
   oligodendrocytes in the CON and CCI groups. (D) UMAP plots showing
   expression status of individual subpopulations of oligodendrocytes from
   different sample sources. (E) Heat map of differential proportions of
   oligodendrocyte subsets in the CON and CCI groups. *, P < 0.05; ***,
   P < 0.001; CCI group vs. CON group. (F) UMAP plots showing expression
   patterns of proliferative OPCs, OPC-1 and Oligo-2 stages. (G) Box plots
   of trends in cell proportions within each oligodendrocyte subpopulation
   in the CON and CCI groups. (H) Bubble plot showing GO enrichment
   pathway analysis for each oligodendrocyte subpopulation. (I) Bar graph
   of the top-90 GO enrichment pathway terms for Oligo-2 stage. (J)
   Kinetic trend graph showing pseudotime analysis of Oligo-2 marker genes

   We verified this inference at the gene level as achieved using a
   cluster analysis of GO functions across the subgroups. We found that
   myelin sheath-related GO functions were markedly expressed after the
   Oligo-1 stage (Fig. [176]7H). Upon further investigation of the
   biological properties of Oligo-2, we observed that myelin sheath
   generation, axon extension, cell metabolism and signal transduction
   were all enriched in Oligo-2 cells during differentiation after CCI, as
   based on results from our GO enrichment analysis (Fig. [177]7I). Based
   on these findings, we examined the marker gene changes in Oligo-2 stage
   along with pseudotemporal trajectories (Fig. [178]7J), with the result
   that the expressions of Apoe, Glul, Pex5l and Ptgds showed similar
   trends as those presented in Fig. [179]7G.

   When collating these findings, it seems that Oligo-2 acted as a node in
   blocking the oligodendrocyte differentiation and myelin sheath
   regeneration after CCI.

Integration analysis confirms oligodendrocyte differentiation arrest after
CCI

   We have provided robust evidence that oligodendrocytes were blocked at
   the Oligo-2 stage after CCI using two integrated and validated
   scRNA-seq datasets to assess the consistency of these data. When
   combining the integrated results for the spinal cord oligodendrocyte
   subpopulations in this study with those as reported by Russ et al.,
   [[180]25]it can be seen that the distribution of the spinal cord
   oligodendrocyte subpopulations was, in general, consistent in the two
   studies (Fig. [181]8A). We analyzed the expression of relevant marker
   genes in the oligodendrocyte subpopulations between this and previously
   published studies. With this comparison we found a good agreement
   between the subpopulations in the two studies, where the Oligo1-Oligo3
   stage matched OProg-2 in the published study (Fig. [182]8B), which was
   in accord with the results from the dendrogram and heat map analyses
   (Fig. [183]8C).

Fig. 8.

   [184]Fig. 8
   [185]Open in a new tab

   Consistency analysis of mouse spinal cord tissue oligodendrocytes with
   published validated cell types. (A) UMAP plots showing integration of
   oligodendrocyte data and each subpopulation in this and published
   study. (B) Dot plot showing expression status of oligodendrocyte marker
   genes in each oligodendrocyte subpopulation of this and published
   study. Marker genes of oligodendrocytes were from published Russ et al.
   (C) Dendrogram and heat map showing similarities between this and
   published study for each subpopulation of oligodendrocytes. (D) UMAP
   plots showing integration of data for all oligodendrocytes and each
   subpopulation in this and another published study. (E) Heat map showing
   correlation analysis of RNA-seq data for each subpopulation of
   oligodendrocytes in this and published study. (F) Dendrogram showing
   similarities between this and published study for each subpopulation of
   oligodendrocytes

   Similarly, when comparing the oligodendrocytes within spinal cord
   tissues between our data and that of Zeisel et al. [[186]22]. , the
   comparisons of these data revealed that a large extent of integration
   connections were present for all oligodendrocytes (Fig. [187]8D).
   Accordingly, we further evaluated whether any notable relationship
   exists between cell subpopulations in these studies (Fig. [188]8E). As
   compared with the results from a previously published study, Oligo-2/3
   and MOL1 were found to belong to the same stage, as indicated by a
   phylogenetic tree analysis (Fig. [189]8F). Notably, although the
   oligodendrocyte subpopulations in our current study did not perfectly
   match with those of the previously published results, a relatively good
   concurrence was present within these results.

   Therefore, our overall characterization of oligodendrocyte
   subpopulations was in good agreement with that of previously published
   studies. Such findings reveal that the significantly increased
   oligodendrocyte cell numbers in the Oligo-2 stage in our study were not
   due to the emergence of a new cell subtype, rather oligodendrocytes
   appear to be blocked at the Oligo-2 stage of differentiation.

Analysis of cell communication reveals that CCI enhanced glial cell
interactions in mouse spinal cord tissue

   As an approach to examine the roles of each cell type in blocking
   oligodendrocyte differentiation we first evaluated the interaction
   intensities between different cell types and their subpopulations in
   CON and CCI spinal cord tissue (Fig. [190]9A). Results of this assay
   demonstrated that there was an increase in spinal cord glial cell
   interactions after CCI.

Fig. 9.

   [191]Fig. 9
   [192]Open in a new tab

   Analysis of interactions in different cell types and glial
   subpopulations in CON and CCI groups. (A) Heat map showing the number
   of potential receptor–ligand pairs between different cell types in
   spinal cord tissue of the CON and CCI groups. (B) Communication network
   showing interactions of different cell types in spinal cord tissue of
   the CON and CCI groups. (C, D) Communication network showing
   differences in numbers and strength of receptor–ligand pair-based
   interactions among cells in spinal cord tissue of the CON and CCI
   groups. Blue represents CON group and red the CCI group. (E, F)
   Communication network showing interactions between subpopulations of
   glial cells in spinal cord tissue of the CON and CCI groups. The left
   side of figure represents all receptor–ligand pairs and the right side
   represents the 40% receptor–ligand pair cut-off

   We further analyzed the differences in receptor–ligand interactions
   among various cell clusters between the CON and CCI groups. We found
   that receptor–ligand interactions associated with spinal cord function
   were reduced after CCI (Fig. [193]9B). However, the number and
   intensity of interacting receptor–ligand pairs among glial cells in the
   CCI group, particularly astrocytes and oligodendrocytes, were superior
   to that observed in the CON group (Fig. [194]9C-D).

   The interactions among various glial cell subpopulations were then
   assessed, after removing 40% of the receptor–ligand pairs with weaker
   interaction intensities as based on all receptor–ligand pairs. In this
   way, it was possible to achieve a convenient visual analysis
   (Fig. [195]9E-F). The results showed that at the Oligo-2 stage, a
   receptor was subjected to stronger ligand actions after CCI, with both
   astrocytes and microglia being involved at this stage.

   These results indicate that microglia and astrocytes are involved in
   the Oligo-2 stage of differentiation arrest in oligodendrocytes after
   CCI.

Cellular interaction analysis and IF staining reveal key targets of glial
cell action in myelin sheath production, blockade and demyelination

   To validate the role of glial cell interactions in Oligo-2
   differentiation arrest, we evaluated receptor–ligand interaction pairs
   among different glial cells in the CON and CCI groups to determine
   target characteristics related to differentiation arrest. Initially,
   receptor–ligand interactions between astrocyte and oligodendrocyte
   subpopulations in the CON and CCI groups were analyzed. In response to
   CCI, the signaling interactions of CADM1/CADM1 and NRP-1/VEGFA were
   enhanced between the astrocytes and oligodendrocytes (Fig. [196]10A-B).
   GO enrichment analysis of astrocytes indicated that pathways associated
   with cholesterol synthesis were downregulated after CCI
   (Fig. [197]10C).

Fig. 10.

   [198]Fig. 10
   [199]Open in a new tab

   Analysis of receptor–ligand interactions among spinal cord glial cells
   in CON and CCI groups. Bubble diagram showing differential expression
   patterns of receptor–ligand interaction pairs between astrocyte and
   oligodendrocyte subpopulations in the CON and CCI groups: (A) Astrocyte
   ligands, (B) Oligodendrocyte ligands. (C) Bar plots showing GO pathway
   annotation in astrocytes in the CON and CCI groups. (D-F) Bubble
   diagram showing differential expression patterns of receptor–ligand
   interaction pairs between microglia and oligodendrocyte subpopulations
   in the CON and CCI groups. Ligands are Microglia-1, -2 and − 8.
   Receptors are oligodendrocytes. (G) UMAP plots showing distributions of
   differentially expressed marker genes in spinal cord tissue. (H) Bubble
   plots showing expression status and interactions of receptor–ligand
   pairs of glial cell subsets in the CON and CCI groups. (I) Heatmap
   showing interactions between microglia and astrocyte subpopulations in
   the CON and CCI groups. (J) Bubble diagram showing expression status
   and interactions of receptor–ligand pairs for microglia and astrocytes
   in the CON and CCI groups. (K) Representative images of IF staining for
   CADM1 in spinal cord tissue of the CON and CCI groups. From left to
   right: CADM1, CADM1 channel image; DAPI, DAPI channel image: Merge,
   merged image between CADM1 and DAPI. Scale bar = 100 μm. n = 3; 3
   sections per animal. (L) Representative VEGFA and NRP-1 staining images
   of spinal cord sections from the CON and CCI groups. From left to
   right: VEGFA, VEGFA channel image; NRP-1, NRP-1 channel image; DAPI,
   DAPI channel image; Merge, merged image among VEGFA, NRP-1 and DAPI.
   Scale bar = 100 μm. n = 5. (M) Bar graphs showing quantification of
   fluorescent images of CADM1, NRP1 and VEGFA as normalized within the
   CON and CCI groups. Data are presented as means ± SEMs. ###, P < 0.001;
   ####, P < 0.0001

   We subsequently examined the differences in receptor–ligand
   interactions between microglia and oligodendrocytes (Fig. [200]10D-F).
   Similar to that as observed in astrocytes, interaction intensities of
   CADM1/CADM1 and NRP-1/VEGFA between microglia and oligodendrocytes were
   enhanced.

   We next analyzed the expression of DEGs in myelin sheath tissue during
   cellular interactions (Fig. [201]10G) and re-evaluated receptor–ligand
   interactions among glial cells in the CON and CCI groups
   (Fig. [202]10H). CADM1/CADM1 interactions between astrocytes and
   oligodendrocytes and the NRP-1/VEGFA interactions between
   Microglia1-Microglia2 and oligodendrocytes were all significantly
   enhanced after CCI. In addition, a NRP-1/VEGFA communication
   relationship was established between Microglia-8 and oligodendrocytes
   after CCI.

   Interactions between astrocytes and microglial subpopulations
   (Fig. [203]10I) and the expression of NRP-1/VEGFA in the two types of
   glial cells were then analyzed in the CON and CCI groups
   (Fig. [204]10J). While the expression of NRP-1 ligands in microglia and
   VEGFA receptors in astrocytes decreased after CCI, NRP-1/VEGFA
   interactions between Microglia-8 and astrocytes were enhanced as
   compared with that observed in the CON group.

   To verify that CADM1/CADM1 and NRP-1/VEGFA interactions among the three
   glial cell types after CCI were consistent with that of scRNA-seq
   results, we used IF staining to assess the expression and localization
   of CADM1, NRP-1 and VEGFA in spinal cord tissue (Fig. [205]10K-M).
   CADM1 expression was significantly increased after CCI (P < 0.0001) and
   was concentrated in the cytoplasm of glial cells, with small
   intercellular spaces and tight connections (Fig. [206]10K, M). When
   using the same protocol to determine the expression and localization of
   NRP-1 and VEGFA in spinal cord tissue, we found that levels of NRP-1
   (P < 0.0001) and VEGFA (P < 0.001) in the CCI group were significantly
   greater than those in the CON group, with enhanced interaction signals
   being primarily localized to glial cells (Fig. [207]10L, M).

   In summary, the enhanced CADM1/CADM1 and NRP-1/VEGFA interactions
   between the spinal cord-tissue glial cells after CCI were key in
   blocking myelin sheath development and demyelination.

Discussion

   NP, which has a complex pathogenesis, is a somatosensory injury or
   disease primarily resulting from abnormal glial cell activation and
   ectopic neuronal discharge. In our mouse model, findings from our
   scRNA-seq, TEM and molecular assays revealed that neuroglial cell
   interactions were enhanced after CCI. Moreover, myelin sheath lesions,
   due to an abnormal differentiation of oligodendrocytes, can lead to
   development of the pain experienced with this condition. In response to
   CCI, there was a heterogeneous expression of several subpopulations of
   astrocytes, microglia and oligodendrocytes in spinal cord tissue and
   NRP-1/VEGFA and CADM1/CADM1 receptor–ligand binding enhanced the
   interactions among these three cell types. Additionally,
   oligodendrocytes affected myelin sheath development by blocking Oligo-2
   subtype differentiation, thus leading to neuronal exposure and ectopic
   firing, critical effects which can then promote NP development.

   In our study, 39,711 cells were clustered and dimensionally reduced by
   UMAP from the data resulting from our scRNA-seq assay and were then
   divided into 33 cell subgroups. Finally, 16 major subsets were
   identified, of which glial cells accounted for the largest proportion.
   We also found a heterogeneous expression of 12 subpopulations of
   oligodendrocytes, 9 subpopulations of microglia and 2 subpopulations of
   astrocytes in the spinal cord tissue after CCI. Such findings suggest
   that glial cells played important roles in NP development, however, it
   was unclear as to whether the effects of these cells involved cellular
   interactions.

   Activated microglia and astrocytes can induce neuroinflammation by
   excessive production of pro-inflammatory mediators, ultimately leading
   to cognitive dysfunction, brain edema and secondary brain damage
   [[208]51, [209]52], effects which can reflect interactions occurring
   between microglia and astrocytes. Gap junction-mediated communication
   between connexins expressed in oligodendrocytes and astrocytes is
   important for myelin sheath formation, as a loss of connexins promotes
   demyelination, which then affects myelin sheath function [[210]53]. In
   the thalamus, oligodendrocytes perform their functions mainly by
   assisting astrocytes in transferring metabolites to synapses [[211]54].
   These findings suggest that interactions between astrocytes and
   oligodendrocytes are critical for demyelination.

   Overall, glial cells play an essential role in neurological disorders
   through various modes of action, however, their precise role and
   mechanisms in NP remain unclear. From our scRNA-seq data, we found that
   enhanced levels of associations and mechanisms of action for the three
   glial cells in the CCI versus CON group. In particular, notable
   increases in the number and intensity of the receptor–ligand pairs
   interacting between astrocytes and oligodendrocytes were observed in
   response to CCI. And, interactions among astrocytes-1, astrocytes-2,
   microglia-2 and Oligo-2 were stronger when both astrocytes and
   microglia were involved.

   Oligodendrocytes are myelin sheath-forming cells of the CNS,
   originating from OPCs. These OPCs differentiate into oligodendrocytes
   by extending multiple protrusions that individually wrap around axons
   to generate concentric layering of modified cell membranes which then
   form myelin sheaths [[212]55]. The structural and functional integrity
   of the myelin sheath depends on the coordination of its components,
   which include lipoproteins, phospholipids, cholesterol, cephalin and
   basic proteins, all serving as the basis for its functional integrity
   [[213]56]. An absence or hypoplasia of the myelin sheath can lead to
   neuron exposure, which can then trigger ectopic firing and produce pain
   [[214]57]. Accordingly, we hypothesized that myelin sheath lesions, as
   resulting from an abnormal differentiation of oligodendrocytes, may
   trigger NP. In support of this hypothesis were our results
   demonstrating that normal oligodendrocyte proliferation and
   differentiation were inhibited by CCI. In addition, GO enrichment
   analysis showed that Oligo-2 cells were enriched in myelin sheath
   production and other aspects of oligodendrocyte differentiation after
   CCI, whereas Oligo-2 marker genes showed a decreasing trend in the
   pseudotime trajectory. Moreover, results from our integrated analysis
   verified that oligodendrocyte differentiation and myelin sheath
   regeneration were blocked at Oligo-2, supplying robust evidence in
   support of this hypothesis.

   Plp1, a major component of CNS myelin, is essential for the axonal
   support function of myelin, with different types of Plp1 mutations
   resulting in Pelizaeus–Merzbacher disease and allelic disease of type
   II spastic paraplegia [[215]58]. Our single-cell clustering analysis
   results demonstrated that Plp1 was differentially expressed in
   oligodendrocyte subtypes after CCI and was involved in NP development.
   The peripheral myelin protein, PMP22, which is mainly expressed in the
   dense myelin sheath of the peripheral nervous system has been shown to
   be involved in orthomyelin formation, myelin sheath stability and
   axonal maintenance during peripheral nerve development [[216]59].
   Ptprz, which encodes a receptor-like protein tyrosine phosphatase,
   exerts a negative role in oligodendrocyte differentiation in early CNS
   development and myelin sheath regeneration in CNS demyelinating
   diseases [[217]60]. Here, we found a differential expression of Ptprz1
   and Pmp22 in oligodendrocytes after CCI, effects which may be
   associated with the blockade of myelin development. These results
   suggest that a differential expression of Plp1, Ptprz1 and Pmp22 in
   oligodendrocytes after CCI may be related to the formation and
   development of the myelin sheath and production of myelin sheath
   lesions.

   Results from our subclass analysis, revealed that specific interactions
   exist between Astrocytes-1 and Oligo-2 subtypes in response to CCI.
   These interactions were significantly stronger than that observed in
   the CON group and were largely dependent on the interactions of
   CADM1/CADM1. The CADM family of cell-adhesion proteins regulates myelin
   sheath development by modulating myelination and myelinated axons
   [[218]61]. CADM1, a cell-adhesion molecule expressed in glial cells and
   axons, induces functional synapses by promoting the formation of
   presynaptic terminals and also plays an important role in facilitating
   astrocyte–astrocyte and astrocyte–neuron communications [[219]62,
   [220]63]. We also showed that the Astrocytes-1 pathway was enriched in
   phospholipid and cholesterol synthesis-related pathways in the CON, but
   not in the CCI group. These findings suggest that myelin sheath
   development in the CCI group was influenced by the intensity and
   function of Oligo-2/Astrocytes-1 interactions and that a blocking of
   the Oligo-2 subtype differentiation further promotes abnormal myelin
   sheath development. Thus, CADM1/CADM1-mediated interactions in the
   Astrocytes-1 and − 2 phases within the CCI group primarily mediated
   pain by affecting myelin sheath development. Based on these findings
   CADM1 may serve as a potential therapeutic target for use in the
   treatment of NP.

   Demyelination is a typical manifestation of myelin sheath lesions, and
   various factors contribute to the link between demyelination and degree
   of NP. Findings from our analyses revealed that enhanced interactions
   are present between astrocytes/oligodendrocytes,
   Microglia1-Microglia2/oligodendrocytes and Microglia-8/astrocytes for
   NRP-1/VEGFA. Such results suggest that the strong interactions present
   between these three glial cell pairs were maintained via the
   NRP-1/VEGFA receptor–ligand pair.

   VEGFA is a major factor that initiates and regulates angiogenesis and
   promotes neurogenesis, neuronal migration, neuronal survival and axonal
   guidance during neurodevelopment [[221]31]. Hiratsuka et al. [[222]64]
   found that inhibiting VEGF signaling downregulates the proliferation of
   OPCs through a lysophosphatidylcholine-induced demyelination, effects
   which demonstrate that VEGFA plays a role in demyelination and nerve
   conduction.

   NRP-1, a transmembrane glycoprotein receptor for VEGFA, exerts a
   specific role in guiding vascular bud endothelial tip cells [[223]65].
   Interactions between NRP-1 and VEGFA enhance signaling and promote
   angiogenesis [[224]66]. In systemic sclerosis (SSc), the reduced
   expression of NRP-1 disrupts the VEGFA/VEGFR2 system, leading to
   peripheral microangiopathy and defective angiogenesis in SSc [[225]67].
   These findings suggest that the interactions present between
   NRP-1/VEGFA play important roles in various nervous system diseases.
   With use of TEM, we demonstrated that demyelination was observed in the
   CCI group. And, results from our scRNA-seq and IF assays suggested that
   NRP-1/VEGFA were responsible for the strong interactions between the
   three glial cell pairs, effects which were then closely associated with
   the demyelination seen in the CCI group. The specific mechanisms of
   action involved in these effects will require further investigation.

   In summary, in this report we demonstrate that oligodendrocyte
   differentiation-mediated myelin sheath lesions and their interactions
   with microglia and astrocytes were closely related to NP development
   and that an intertwined “network” of glial cell DEGs regulate myelin
   sheath lesion formation. Current treatments for NP focus on
   pharmacotherapy, physiotherapy and surgery [[226]68], all of which
   require prolonged treatments, high rates of trauma and/or
   unsatisfactory results. Results from our single-cell clustering
   analysis identified some of the neurobiological changes underlying NP.
   Such findings provide important, new directions and predictive targets
   to achieve maximal analgesic effects in the treatment of NP. In
   specific, CADM1 inhibitors could reduce the intensity of interactions
   between Astrocytes-1 and Oligo-2 stages to improve myelin sheath
   development and slow NP development. Second, based on the associations
   observed between demyelination and NP, various drugs, including
   olesoxime [[227]69], clemastine [[228]70], GnbAC1 [[229]71] and
   etazolate [[230]72], which have been employed to alleviate
   demyelination and treat multiple sclerosis, could be promising
   treatment strategies. While the findings present in this report offer
   great promise for future applications, it is important to note that
   practical applications of the molecular targets identified in this
   study will require further evaluation in follow-up studies.
   Additionally, as the spatial information was lost during the
   single-cell dissociation, the expression distribution and function of
   glial cells and ligand–receptor pairs in spinal cord tissue after CCI
   could not be identified. Accordingly, a follow-up study will need to be
   performed to validate intercellular communication mechanisms by spatial
   transcriptomics and further elucidate the pathogenic mechanisms of NP
   as based on the spatial information of cells in spinal cord tissue.

Conclusion

   Our current results indicate that myelin lesions were present in the
   spinal cord of the CCI mouse model and that myelin sheath development
   was blocked or demyelinated. From our scRNA-seq assay we demonstrate
   that oligodendrocyte development in the spinal cord of our NP mouse
   model was blocked and could mediate myelinopathy through CADM1/CADM1
   and NRP-1/VEGFA receptor–ligand interactions. Such receptor–ligand
   interactions, in conjunction with astrocyte and microglia interactions,
   can then trigger the pain associated with NP. Based on these findings,
   it seems that drugs which can inhibit CADM1 or other demyelinating
   disorders may provide new clinical options for NP treatment.

Electronic supplementary material

   Below is the link to the electronic supplementary material.
   [231]Supplementary Material 1^ (9.8MB, txt)
   [232]Supplementary Material 2^ (18.6KB, xlsx)
   [233]Supplementary Material 3^ (2.6MB, xlsx)
   [234]Supplementary Material 4^ (7.9KB, xlsx)

Acknowledgements