Single-cell sequencing reveals glial cell involvement in development of neuropathic pain via myelin sheath lesion formation in the spinal cord

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.

Background

Neuropathic pain (NP) is a chronic, persistent and intractable pain resulting from injuries or lesions of the somatosensory nervous system [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 [2]. Accordingly, NP substantially reduces the quality of life and increases the public health burden [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 [4, 5]. Astrocytes regulate NP through changes in glial signaling pathways, receptors and channel protein expression and glia-derived factors [6]. In addition, interactions mediated by IL-18/IL-18R between microglia and astrocytes have also been shown to promote NP [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 [8]. However, current studies of oligodendrocytes have primarily focused on the mechanisms of central nervous system (CNS) axonal myelination [9]. Microstructural changes in the myelin sheath contribute to aberrant neuronal firing by affecting nerve conduction and functions of neuronal circuits [10, 11]. NP, as resulting from abnormal spontaneous neuronal firing and peripheral nerve demyelination, represent key factors in the development of diabetic NP [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 [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 [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 [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 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 [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 [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 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₂ 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 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₂ FC of 0.2 and a maximum adjusted P-value (via the Bonferroni method) of 0.05 were considered as DEGs [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 [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 [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 4 [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 [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 [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. 1A.

Fig. 1

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. 1B-C, P > 0.05) or over any of the time periods sampled in the CON group (Fig. 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. 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. 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. 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. 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 [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. 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. 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. 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.

Fig. 2

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 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. 2A). As grounded on results from a previous study [25], the original subcluster was merged into 16 main subsets based on the most-significant DEGs (Fig. 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. 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. 2D). Notably, cells identified as OPCs, oligodendrocytes and microglia belonged to multiple clusters, suggesting that these cell types may share a heterogeneity (Fig. 2C).

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. 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 [25, 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. 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. 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. 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. 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.

Fig. 3

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

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. 2A).

We identified 12 oligodendrocyte subpopulations as based on UMAP visualizations (Fig. 3A) and used multiple markers to reveal subpopulation characteristics based on heat mapping of the top 20 DEGs in each subpopulation (Fig. 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).

We found that oligodendrocytes were widely expressed in different samples (Fig. 3D), but expression ratios and cell numbers within each subpopulation varied as a function of the stage examined (Fig. 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. 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. 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. 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. 3K-M). The differential marker genes Opalin, Ptprz1, Klk6 and Pmp22 for each pathway were analyzed as shown in Fig. 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.

Fig. 4

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

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. 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 [27]. Cluster analysis was used to separate the microglia into nine subgroups (Microglia-1 to Microglia-9), each demonstrating unique characteristics (Fig. 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. 4C-D), and their differential expressions were significantly different between subpopulations (Fig. 4E). Results from UMAP analysis showed that microglia within different stages formed distinct heterogeneous clusters (Fig. 4F-G). However, findings from the heat map of the top 10 DEGs also revealed considerable similarities among subpopulations (Fig. 4B).

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. 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. 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. 4J). Pathway expressions of immune-promoting factors and astrocyte development were found to be significantly downregulated in CCI mice (Fig. 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.

Fig. 5

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

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 [28]. From these analyses we found two distinct subpopulations of astrocytes, Astrocytes-1 and Astrocytes-2 (Fig. 5A). The top 20 specific markers for each subpopulation are shown in Fig. 5B. In addition, when investigating the expression of astrocytes in multiple samples and states, as based on specific subpopulation markers (Fig. 5C-F), we observed that cell numbers within the astrocyte subpopulations increased after CCI (Fig. 5D, F). While this increase was greater for Astrocytes-1 (Fig. 5E-F), differences between Astrocytes-1 and − 2 failed to achieve statistical significance with regard to astrocyte subpopulations (Fig. 5G, P = 0.259).

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. 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. 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.

Fig. 6

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

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. 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. 6A), results which were corroborated from the proposed time series as shown in Fig. 6B. Therefore, we evaluated the DEGs showing significant variation (Fig. 6C) and found that Cspg4 [29], Gria2 [30], Gnai1 [31], Sh3kbp1 [32], Tanc2 [33], Gal3st1 [34], and Arntl2 [35] all demonstrated increased levels of gene expression.

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. 6D). Notably, proliferative OPCs tended to be more frequently expressed during the S and G2/M phases (Fig. 6E-F). When separating OPC-1 from the proliferative OPCs, primarily expressed during the cytokinesis cycle, for pseudotime analysis (Fig. 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. 6H) and conducted an RNA velocity analysis on marker genes to characterize the G1, S and G2/M phases (Fig. 6I). With this analysis Zbtb20 [36], Ezh2 [37], Brca1 [38], Hdgf [39], Pex5l [40], Stmn1 [41], and Vegfa [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. 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 [43], Ntrk2 [44], Fads2 [45], Kcnmb4 [46], Rgcc [47], Slc1a1 [48], Dusp10 [49] and Car2 [50] were all significantly expressed in the CON group (Fig. 6L). In the CCI group, Cntn4, Adgrb3, Rgs7, Sox6, Slc1a1, Ttll7 and Car2 were significantly expressed (Fig. 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.

Fig. 7

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

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. 6J-K), we next evaluated the differentiation process of the oligodendrocyte subpopulations in the CON and CCI groups using a clustering analysis (Fig. 7A-B). The visualization of results for the different samples and sources assessed is shown in Fig. 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. 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. 7G).

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. 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. 7I). Based on these findings, we examined the marker gene changes in Oligo-2 stage along with pseudotemporal trajectories (Fig. 7J), with the result that the expressions of Apoe, Glul, Pex5l and Ptgds showed similar trends as those presented in Fig. 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.

Fig. 8

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

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., [25]it can be seen that the distribution of the spinal cord oligodendrocyte subpopulations was, in general, consistent in the two studies (Fig. 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. 8B), which was in accord with the results from the dendrogram and heat map analyses (Fig. 8C).

Similarly, when comparing the oligodendrocytes within spinal cord tissues between our data and that of Zeisel et al. [22]. , the comparisons of these data revealed that a large extent of integration connections were present for all oligodendrocytes (Fig. 8D). Accordingly, we further evaluated whether any notable relationship exists between cell subpopulations in these studies (Fig. 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. 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.

Fig. 9

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

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. 9A). Results of this assay demonstrated that there was an increase in spinal cord glial cell interactions after CCI.

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. 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. 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. 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.

Fig. 10

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

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. 10A-B). GO enrichment analysis of astrocytes indicated that pathways associated with cholesterol synthesis were downregulated after CCI (Fig. 10C).

We subsequently examined the differences in receptor–ligand interactions between microglia and oligodendrocytes (Fig. 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. 10G) and re-evaluated receptor–ligand interactions among glial cells in the CON and CCI groups (Fig. 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. 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. 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. 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. 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. 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 [51, 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 [53]. In the thalamus, oligodendrocytes perform their functions mainly by assisting astrocytes in transferring metabolites to synapses [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 [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 [56]. An absence or hypoplasia of the myelin sheath can lead to neuron exposure, which can then trigger ectopic firing and produce pain [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 [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 [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 [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 [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 [62, 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 [31]. Hiratsuka et al. [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 [65]. Interactions between NRP-1 and VEGFA enhance signaling and promote angiogenesis [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 [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 [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 [69], clemastine [70], GnbAC1 [71] and etazolate [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.