Mitochondrial dysfunction and disulfidptosis co-regulate neuronal cell in neuropathic pain based on bioinformatics analysis

Abstract

Neuropathic pain (NP) affects approximately 6.9–10% of the world’s population and necessitates the development of novel treatments. Mitochondria are essential in the regulation of cell death. Neuroimmune mechanisms are implicated in various forms of cell death associated with NP. However, the specific involvement of mitochondrial dysfunction and disulfidptosis in NP remains uncertain. Further research is required to gain a better understanding of their combined contribution. Our comprehensive study employs a variety of bioinformatic analysis methods, including differential gene analysis, weighted gene co-expression network analysis, machine learning, functional enrichment analysis, immune infiltration, sub-cluster analysis, single-cell dimensionality reduction and cell-cell communication to gain insight into the molecular mechanisms behind these processes. Our study rationally defines a list of key gene sets for mitochondrial dysfunction and disulfidptosis. 6 hub mitochondrial genes and 3 disulfidptosis-related genes (DRGs) were found to be associated with NP. The key genes were predominantly expressed in neurons and were lowly expressed in the NP group compared to SHAM. In addition, our macrophages used the APP (Amyloid precursor protein)-CD74 (MHC class II invariant chain) pathway to interact with neurons. These results suggest that NP is interconnected with the mechanistic processes of mitochondrial dysfunction and disulfidptosis, which may contribute to clinically targeted therapies.

Introduction

Neuropathic pain (NP) refers to pain that arises from damage or dysfunction in the somatosensory nerve system. ¹ Risk factors associated with NP include conditions such as diabetes, nerve compression, traumatic events, oxidative stress, cancer, and inflammation.^2,3 In terms of therapeutics, the main classes of drugs target the α2δ subunits of calcium channels, sodium channels, and descending modulatory inhibitory pathways. ⁴ However, it is important to note that there is currently a need for more effective and safe pain treatments, leading to ongoing research efforts focused on targeted treatments and expanding treatment options to alleviate the global burden of pain.^5,6

Disulfidptosis is a newly discovered mechanism of cell death, representing an emerging form of programmed cell death. It is characterized by the abnormal accumulation of disulfides, which are chemical bonds formed between sulfur atoms within proteins.^7,8 The occurrence of disulfidptosis is triggered by the upregulation of solute carrier family 7 member 11 (SLC7A11) expression during glucose deprivation, initiating cell death through redox reactions and the formation of disulfide bonds.^8,9 Recent research has highlighted the diverse roles of mitochondria in cellular metabolism, cancer, immune system regulation, tissue homeostasis maintenance, mitochondrial quality control, wound healing, and adipose tissue function. ¹⁰ Mitochondria also play a crucial role in regulating various forms of cell death, including necroptosis, ferroptosis, and pyroptosis. ¹¹ Studies have indicated that neuroimmune mechanisms involved in neuropathic pain may be associated with multiple forms of cell death, such as necroptosis, pyroptosis, and ferroptosis.^12–15 Investigating the role of mitochondrial dysfunction and disulfidptosis in neuropathic pain could provide insights for the development of improved therapeutic strategies. However, whether mitochondrial dysfunction and disulfidptosis is involved in neuropathic pain is unclear and requires further study.

In our study, we utilized RNA sequencing data from rat model and single-cell data obtained from mouse model of NP to investigate the mechanistic connection between disulfidptosis and mitochondria in NP. Our analysis specifically examined the regulation of NP through gene expression and cell cluster expression, aiming to evaluate its diagnostic and therapeutic implications. This research serves as a fundamental basis for further exploration into the diagnostic and therapeutic possibilities of disulfidptosis mechanisms associated with mitochondria in individuals suffering from NP.

Materials and methods

Data source and workflow

We retrieved RNA sequencing data from the Gene Expression Omnibus (GEO) database. Specifically, we downloaded three datasets, namely GSE24982, ¹⁶ GSE63442, and GSE18803, ¹⁷ from GEO (https://www.ncbi.nlm.nih.gov/geo/) to create our entire dataset. To integrate the datasets, we utilized the R package inSilicoMerging. ¹⁸ Furthermore, we employed the removeBatchEffect function from the R package limma (version 3.42.2) to eliminate the batch effect, resulting in a matrix that is free from batch effects. ¹⁹ This integrated dataset consists of 38 samples from rattus norvegicus individuals diagnosed with neuropathic pain (NP) through L4 and L5 Dorsal Root Ganglion (DRG) surgery or sciatic nerve injury (SNI). We also included 38 samples from neurologically healthy individuals serving as SHAM samples. To clarify, the GSE24982 dataset consists of 20 pairs of ipsilateral L4 or L5 DRG tissue samples obtained from rats that underwent either sham surgery or Spinal Nerve Ligation (SNL) surgery. In the GSE63442 dataset, there are 6 samples from sham surgery and 6 samples from SNL surgery performed on rat DRG tissue. Lastly, the GSE18803 dataset includes expression data from 12 pairs of samples, comparing SNI samples to non-SNI samples. The scRNA-seq dataset GSE174430 for NP, which contains scRNA-seq datasets for 2 SNI and 1 SHAM sample, was downloaded from the GEO database. Our analysis was built upon a comprehensive dataset, enabling us to investigate multiple facets including weighted gene co-expression network analysis (WGCNA), differential gene expression analysis (DEGs), machine learning, immune infiltration, single-cell dimensionality reduction clustering, and cell communication. These analyses were conducted to gain crucial insights into the molecular mechanisms driving NP. It is worth mentioning that our study solely relied on publicly available data and did not involve any experimental work on humans or animals conducted by the authors.

Identification of disulfidptosis-related genes (DRGs) and mitochondria-related genes in NP

To identify the relevant DRGs, we conducted a comprehensive review of recent literature studies.^7,8,20 Based on the findings from this study, we initially identified 33 candidate genes associated with disulfidptosis. In addition, 1140 mouse genes encoding proteins that strongly support mitochondrial localization were acquired from the MitoCarta3.0 (https://www.broadinstitute.org/mitocarta). By combining the candidate DRGs identified from the literature review and the mouse Mito genes from MitoCarta3.0, we aimed to define the disulfidptosis/mitochondria-related gene set for exploring associations with NP disease.

Functional enrichment analysis of DRGs and mito: insights from KEGG and GO enrichment

Functional enrichment analysis is a valuable method for studying gene function and identifying relevant genomic information. In our research, we utilized various strategies to conduct functional enrichment analysis. To explore the pathways and molecular mechanisms associated with DRGs and Mito, we employed the Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis and Gene Ontology (GO) annotations. For gene set KEGG enrichment analysis, we leveraged the KEGG REST API available at (https://www.kegg.jp/kegg/rest/keggapi.html) to access the latest gene annotations of KEGG Pathway as the background. We mapped the genes to the background set and performed the enrichment analysis using the clusterProfiler R software package (version 3.14.3) to obtain gene set enrichment results. ²¹ For gene set GO enrichment analysis, we used the GO annotations of the genes from the org.Hs.eg.db R package (version 3.1.0) as the background. Similarly, we mapped the genes to the background set and conducted the enrichment analysis using the clusterProfiler R package (version 3.14.3) to obtain the gene set enrichment results. In both cases, p < .05 was considered statistically significant.

WGCNA: identifying NP-related gene networks

In order to identify biologically relevant gene networks and co-expressed modules associated with the disease, we utilized weighted gene co-expression network analysis (WGCNA) method. ²² Initially, we conducted gene expression profiling and calculated the Median Absolute Deviation (MAD) for each gene. Subsequently, we excluded the bottom 50% of genes with the smallest MAD values. Furthermore, we employed the “goodSamplesGenes” method from the R package WGCNA to remove outlier genes and samples. This preprocessing step allowed us to construct a scale-free co-expression network. Using WGCNA, we identified potential modules that were associated with different subgroups of NP and sham samples based on the integrated GEO data. We then focused on the significant module genes for further analysis, which ultimately led us to identify hub genes that play crucial roles in the disease.

Characterization of NP-related DEGs

For the identification of DEGs associated with NP, we employed the limma method, which is based on a generalized linear model. ²³ Specifically, we utilized the R package limma (version 3.40.6) for conducting the differential analysis based on the integrated GEO data. By applying appropriate cutoff criteria, specifically |logFC| > 1 and a p value <.05, we identified DEGs that exhibited significant changes in their expression levels. These criteria ensured that we selected DEGs with substantial fold changes and statistical significance. The resulting set of DEGs provides valuable insights into the alterations in gene expression associated with neuropathic pain and may contribute to a better understanding of the underlying mechanisms of this condition.

Analysis of immune cells

Analyzing the changes in immune cells after NP is essential for understanding disease progression and treatment response. The process of immune cell infiltration analysis involved the utilization of the ssGSEA algorithm provided in the R package GSVA. ²⁴ Specifically, we employed a set of 24 immune cell markers to calculate the degree of immune cell infiltration. This approach allowed us to assess the extent of immune cell presence in the analyzed samples. By applying the ssGSEA algorithm and utilizing the specific immune cell markers, ²⁵ we were able to quantify and evaluate the immune cell infiltration status in the context of our study. The correlation coefficients obtained using the “spearman” correlation method allowed us to explore potential associations between immune cell populations and the identified central genes. By conducting a comprehensive analysis of immune cell composition and their interactions with hub genes, our goal was to uncover valuable insights into the immune response and its involvement in the development of NP.

Machine learning for hub gene selection

In our study, we combined the Random Forest algorithm and LASSO regression to conduct gene selection and identify key genes. Random Forest is a highly effective machine learning algorithm that employs an ensemble of decision trees to achieve both accuracy and flexibility. For the implementation of Random Forest, we utilized the “randomForest” R package. ²⁶ Furthermore, we utilized the LASSO regression method to enhance prediction accuracy and feature selection. ²⁷ LASSO regression applies regularization techniques to select relevant variables and improve prediction precision. To perform the LASSO regression, we utilized the “glmnet” R package and employed the lasso-cox method. In order to optimize our model and evaluate its performance, we set the Lambda value to 0.09. Based on our findings, we identified 6 genes that exhibited strong predictive potential and ranked them accordingly. These selected genes are considered as key genes in our analysis.

Sub-cluster analysis

ConsensusClusterPlus, a class discovery tool with confidence assessments and item tracking developed by Wilkerson and Hayes in 2010, ²⁸ was utilized to perform cluster analysis. The analysis employed agglomerative pam clustering with a 1-Pearson correlation distance metric. The dataset consisted of 38 NP samples and 9827 genes. Resampling was performed on 80% of the samples, and the process was repeated 100 times. The optimal number of clusters was determined by examining the empirical cumulative distribution function plot.

scRNA-seq data processing and analysis

Using the R package Seurat ²⁹ (version 4.3.0.1), scRNA-seq data were converted to Seurat objects by the CreateSeuratObject() function and the Percentage of mitochondrial genes and ribosomes were calculated using the PercentageFeatureSet() function. The raw data was filtered according to the criteria for rejecting ineligible data: a minimum cutoff of 200 genes per cell was set, each gene was expressed in at least three cells, and 9788 cells passed the screening for subsequent analysis. The data were then normalized using the “LogNormalize” method. The FindVariableFeatures() function was used to extract the top 2000 genes with high coefficients of variation between cells and display the top 10 genes. To identify the most significant principal component (PC) values for cell clustering, we employed Dimheatmap, ElbowPlot, and clustree techniques. Subsequently, we utilized the FindClusters function with a resolution of 1 to determine the clusters. Cell clusters were categorized using t-Distributed Stochastic Neighbor Embedding (t-SNE) and uniform manifold approximation and projection (UMAP), ³⁰ and maker genes for each cluster were filtered using the FindAllMarkers() function with a cutoff value of logfc.threshold = 0.25; 1, Cell population expression ratio >0.25, adjusted p-value <.05. Cell cluster annotation was performed using the SingleR ³¹ (Version 1.10.0) package, with references loaded from the celldex package (Version 1.8.0) MouseRNAseqData.RData(). Cell differentiation trajectories were reconstructed via the monocle ³² (Version 2.26.0) package. Depending on the marker genes that differed between clusters, dimensionality reduction was performed by the reduceDimension function, where reduce_method = “ DDRTree” and max_components = 2. Characterized genes for different cell states were analyzed with |log2FC| > 1 and adjusted p-value <.05. In addition, the “CellChat” R package ³³ (1.6.0) was utilized to infer the cell-cell communication.

Statistical analysis and visualization of results

This study’s statistical analysis and visual plots were performed using R version 4.2.1, with visualizations created using the ggplot2 package (version 3.3.6). In addition to the packages above, several other bioinformatic analyses were conducted. For ROC analysis, the qROC package (version 1.18.0) was utilized. Differences between groups were assessed using the Wilcoxon rank sum test, with statistical significance set at p < .05 (*p < .05; **p < .01; ***p < .001; ****p < .0001).

Results

Batch effect removal and dataset integration

To enhance the reliability and stability of the NP dataset, we merged the GSE24982, GSE63442, and GSE18803 datasets to obtain batch-corrected data. The density plots (Figure 1(a) and (b)) clearly demonstrate significant differences in sample distributions among the datasets prior to the removal of the batch effect, indicating the presence of batch effects. However, after the removal of the batch effect, the data distributions between the datasets converged, suggesting successful elimination of the batch effect. Similarly, the UMAP plots (Figure 1(c) and (d)) illustrate that before the removal of the batch effect, the samples from each dataset exhibit distinct clustering, further confirming the existence of the batch effect. Conversely, after the removal of the batch effect, the samples from each dataset cluster together and intertwine, indicating a more effective elimination of the batch effect. Finally, the comparison box plots (Figure 1(e) and (f)) provide a comprehensive visualization of the samples before and after batch correction, demonstrating a satisfactory removal of batch effects. Overall, integrating three small NP samples to remove batch effects offers the benefits of reducing technical variability, improving data comparability, and enhancing the accuracy of analysis results.

Figure 1.

GEO data batch cancellation. (a) Density plots for the three datasets before de-batching. (b) Density plots for the three datasets after de-batching. (c) Uniform manifold approximation and projection (UMAP) analysis for samples before de-batching. (d) UMAP analysis datasets after de-batching. (e) Boxplots of gene expression of the dataset before de-batching for three selected GEO Databases. (f) Boxplots of gene expression of the integrated dataset after de-batching for three selected GEO databases. The GSE24982 dataset is represented by the red color, the GSE63442 dataset by the green color, and the GSE18803 dataset by the blue color. GSE: GEO series; GSM: GEO Sample.

Identification of DRGs and mito-related genes in NP

To identify key genes associated with NP, we conducted a differential gene expression analysis using the Limma package. The results of the analysis revealed significant differences in gene expression between the NP and sham sample groups. To identify the most significantly dysregulated genes, we constructed a volcano plot, which identified a total of 1694 upregulated genes and 3759 downregulated genes as candidate genes (|log2FC| > 1, p < .05) (Figure 2(a)). Additionally, heat maps were generated to visualize the overall distribution of the top 50 upregulated and downregulated genes (Figure 2(b)). These candidate genes were further investigated as potential key genes in the context of NP. Overall, we identified 5453 DEGs, which will undergo further analysis to gain insights into their roles and potential implications in the pathogenesis of NP.

Figure 2.

Recognizing DEGs and filtering DRGs in NP. DEGs are shown in (a) volcano plot and (b) heat map (|log2(FC)| > 1, p < .05). The heat map demonstrates the top 50 genes that are up- and down-regulated. The volcano plot demonstrates the DEGs between NP and SHAM samples. (c) Module-feature relationships in NP. (d) Hierarchical clustering of similarly expressed NP genes. Different module groups are color-coded. (e) Venn diagram showing the overlap between selected DRGs and mito as candidate genes in DEGs and WGCNA.

In the integrated dataset, we employed WGCNA to identify modules of genes exhibiting similar co-expression patterns. To determine appropriate soft thresholds, we assessed threshold power values (R^2) greater than 0.88 and a high average connectivity of 10. Ultimately, eight modules were retained, with the grey module representing genes that couldn’t be assigned to any specific module. To explore the correlation between modules and clinical characteristics, we calculated the correlations between module eigengenes and the NP and sham groups. Among these modules, the blue module (r = 0.64) and the brown module (r = 0.62) exhibited the strongest positive correlations with NP (Figure 2(c)). The clustering tree depicted the merging and separation of the modules (Figure 2(d)). We identified 2370 hub genes within the clinically significant modules, which displayed high connectivity. Specifically, there were 674 hub genes in the blue module and 1696 hub genes in the brown module. Further analysis was performed on all genes within these two modules to uncover their potential implications in NP.

In our study, we performed a comprehensive analysis to identify hub genes associated with NP. This analysis involved intersecting candidate gene sets derived from 1140 Mito-related genes, 33 DRGs, WGCNA, and DEGs. Venn diagram analysis revealed a total of 3759 DEGs and 2370 WGCNA trait genes, with 171 hub genes overlapping with Mito-related genes and 4 hub genes overlapping with the DRGs (Figure 2(e)).

Functional enrichment analysis provides insights into disulfidptosis and mitochondrial dysfunction

To gain deeper insights into the functional implications of disulfidptosis and mitochondrial dysfunction, we conducted enrichment analysis using the KEGG and GO datasets. The KEGG analysis of DRGs revealed associations with various cellular processes, disease infections, and organismal systems, including the Regulation of actin cytoskeleton, Pathogenic Escherichia coli infection, Salmonella infection, Focal adhesion, among others (Figure 3(a)). Additionally, the GO analysis of DRGs highlighted relevant biological processes (BP), cellular components (CC), and molecular functions (MF), such as hexose transmembrane transport, monosaccharide transmembrane transport, carbohydrate transmembrane transport, and carbohydrate transport (Figure 3(b)). Furthermore, the KEGG analysis of genes related to mitochondrial dysfunction showed associations with human diseases and energy metabolism of organisms, such as Metabolic pathways, Thermogenesis, Parkinson’s disease, Fatty acid degradation, and more (Figure 3(c)). Similarly, the GO analysis of mitochondrial-related genes revealed relevant terms including mitochondrion, mitochondrial part, and mitochondrial envelope (Figure 3(d)).

Figure 3.

Enrichment analysis in NP. (a) KEGG and (b) GO concentrated bubble chart in 33 DRGs. (c) KEGG and (d) GO concentrated bubble chart in 171 mito genes.

These enrichment analyses provide valuable insights into the functional implications of disulfidptosis and mitochondrial dysfunction, highlighting their involvement in various cellular processes, disease associations, and biological functions. These initial exploratory analyses provide supporting evidence for the definition of DRGs and their associated pathway mechanisms.

Evaluation of hub genes and validation of gene expression in NP

To assess the importance and reliability of hub Mito genes, we utilized machine learning algorithms, namely LASSO and Random Forest, for gene screening. Using the LASSO algorithm on the consolidation dataset, we identified the top 15 genes with high significance at a Lambda value of 0.09 (Figure 4(a)). Consistent results were obtained with the Random Forest algorithm, highlighting TOMM5, PDHB, MRPL39, SIRT3, MRPS15, and CPT1C as important genes (Figure 4(b)). These machine learning approaches provided a robust assessment of the significance and reliability of the hub Mito genes, emphasizing their relevance and influence in the study.

Figure 4.

Machine learning in NP. The most relevant hub Mito genes were selected using the NP dataset through (a) LASSO and (b) RandomForest regression. (c) Expression distributions of six hub mito genes between NP and SHAM groups. (d) ROC curves of six hub mito genes. (e) Expression correlations of six hub genes. (f) Expression distributions of four hub DRGs between NP and SHAM groups. (g) ROC curves of four hub DRGs.

To validate the expression of the identified genes in NP, we examined their differential expression in NP and sham data. The results demonstrated that all six Mito hub genes were differentially expressed and down-regulated in the NP samples, showing statistical significance (Figure 4(c)). Furthermore, these six genes exhibited promising diagnostic value (Figure 4(d)). The chord diagram of correlation analysis revealed a positive correlation among their expressions (Figure 4(e)). Additionally, among the four key genes identified in the Differentially Regulated Genes (DRGs), three genes showed statistically significant expression differences in the NP samples (Figure 4(f)) and displayed good diagnostic value (Figure 4(g)). Specifically, ACTB exhibited high expression in NP, while G6PD and SLC2A8 were expressed at lower levels.

These findings provide strong evidence supporting the dysregulation of these genes in the context of NP, and their differential expression patterns suggest their potential as diagnostic markers for neuropathic pain.

Immune cell composition in NP

To gain insights into the role of the immune system in the development of NP, we investigated the component of immune cells in the brain tissue. Utilizing the ssGSEA algorithm, we investigated the associations between the hub genes and immune cells. B cells, CD8 T cells, Eosinophils, Macrophages, NK cells, T helper cells, Tcm, Th17 cells and Th2 cells were highly expressed in NP samples compared to NK CD56 bright cells, pDC, T cells and Tgd in NP with low expression (Figure 5(a)).

Figure 5.

Analysis of immune infiltration in NP. (a) Comparison of 13 immune cell subtypes between NP and SHAM types. (b) Correlation analysis of 3 hub DRGs and 6 hub mito genes in 13 immune cells.

In our investigation, we explored the correlation between DRGs (ACTB, G6PD, and SLC2A8) and key Mito genes (TOMM5, PDHB, MRPL39, SIRT3, MRPS15 and CPT1C) with differentially expressed immune cells in NP. We observed certain commonalities in the correlations between these genes and immune cell populations.

Specifically, we found negative correlations between G6PD, SLC2A8, TOMM5, PDHB, MRPL39, SIRT3, and MRPS15 with B cells, Tcm cells, and Th2 cells (All r < −0.3, all p < .05) (Figure 5(b)). This suggests a consistent impact of disulfidptosis and mitochondrial dysregulation on the infiltration of these immune cell subsets. Furthermore, SLC2A8, TOMM5, and PDHB showed positive correlations with NK CD5 bright cells and Tgd cells (All r > 0.4, all p < .05) (Figure 5(b)). This indicates a potential association between these genes and the presence or activity of these specific immune cell populations.

These findings highlight the consistent effects of disulfidptosis and mitochondrial dysfunction on immune cell infiltration. The correlations observed between these genes and different immune cell populations provide insights into the potential interplay between neuropathic pain, immune responses, and the dysregulation of key genes involved in disulfidptosis and mitochondrial function.

Immune infiltration of sub-cluster analysis in NP

The study demonstrated maximizing the therapeutic effect by stratifying patient subgroups in a neuropathic pain trial. ³⁴ In our study, we conducted consistent unsupervised clustering on a dataset consisting of 38 NP (nasopharyngeal) samples from integrated GEO data. Based on the relative stability, we selected two clusters, namely cluster 1 (n = 22) and cluster 2 (n = 16), for further analysis (Figure 6(a)). We observed higher expression levels of macrophages, NK cells, and Th17 cells in cluster 2 compared to cluster 1, while Tgd (presumably referring to a specific cell type) showed higher expression in cluster 1 (Figure 6(b)). In our investigation, we examined the correlation between differentially expressed immune cells in clusters 1 and 2, and key mitochondrial genes (ACTB, G6PD, SLC2A8, TOMM5, PDHB, MRPL39, SIRT3, MRPS15, and CPT1C). Notably, hub mitochondrial genes such as G6PD and PDHB exhibited differential expression in these subgroups, with G6PD showing higher expression in cluster 2 and PDHB showing higher expression in cluster 1 (Figure 6(c)). Consequently, we further analyzed the correlation between the key genes G6PD and PDHB in the four immune cell types. Our findings indicated a positive correlation between G6PD and macrophages, NK cells, and Th17 cells, while a negative correlation was observed with Tgd. On the other hand, PDHB showed a positive correlation with Tgd and a negative correlation with macrophages and NK cells (All correlation coefficient |r| > .4, all p-value <.05) (Figure 6(d)). These findings provide insights into the differential expression patterns of immune cells and the association between key mitochondrial genes and DRGs and immune cell types in the context of NP samples.

Figure 6.

Sub-cluster analysis in NP. (a) Optimal unsupervised clustering analysis (k = 2). (b) Comparison of 4 immune cell subtypes between C1 and C2 cluster. (c) Expression differences of 3 hub DRGs and 6 hub mito genes between groups C1 and C2. (d) Correlation heat map of G6PD and PDHB in 4 immune cells. ***p < .001, ****p < .0001, ns: p > .05.

Single-cell quality control and dimension reduction clustering in NP cell subtypes

In the scRNA analysis of the GSE174430 dataset, a comprehensive pipeline was implemented to ensure data quality, reduce dimensionality, and identify distinct cell clusters. The gene features, counts, mitochondrial gene percentage, and ribosome percentage were assessed for each sample to evaluate sample quality. Adjustments were made to account for the proportion of mitochondrial and erythrocyte genes, ensuring reliable data representation (Figure 7(a)). A total of 9788 single-cell samples were used for downstream analysis. Correlations between per-sample counts and genes, mitochondrial gene percentage, and ribosomal RNA percentage were explored. Pearson correlation analysis of sequencing depth with mitochondrial percentage showed a correlation coefficient of −0.07, a number correlation coefficient of 0.03 with ribosome percentage, and a correlation coefficient of 0.82 with gene number (Figure 7(b)). This analysis provides reliable insights into potential biases or technical artifacts in the data. Highly variable genes among cells were identified using a volcano map, with the top 10 genes labeled (Figure 7(c)). This step allowed for the identification of genes that exhibited significant variation and potential functional relevance. The standard deviation of principal components (PCs) 1–30 was calculated using the ElbowPlot function (Figure 7(d)). This analysis helped determine the optimal number of PCs to retain for subsequent analyses, reducing the dimensionality of the dataset. The t-SNE algorithm was employed to visualize the dataset, reducing the high-dimensional data into a two-dimensional representation. This visualization facilitated the identification and visualization of 22 distinct cell clusters, providing valuable insights into the cellular heterogeneity present in the NP samples (Figure 7(e)). Overall, this scRNA analysis pipeline ensured data quality, reduced dimensionality, and enabled the identification and characterization of distinct cell clusters within the NP dataset.

Figure 7.

Quality control, dimension reduction, and clustering of single cells in NP: scRNA analysis of GSE174430 dataset. (a) Assessment of gene features, counts, mitochondrial gene percentage, and ribosome percentage for each sample to ensure sample quality. Adjustments were made to the proportion of mitochondrial and erythrocyte genes. (b) Exploration of correlations between per-sample counts and genes, mitochondrial gene percentage, and ribosomal RNA percentage. (c) Identification of highly variable genes among cells through a volcano map, with the top 10 genes labeled. (d) Calculation of standard deviation for principal components (PCs) 1–30 using the ElbowPlot function. (e) Visualization of 22 distinct cell clusters using t-SNE algorithm with a resolution of 1.

Seven cell clusters reveal major marker genes and cellular heterogeneity upon integration and characterization

The scRNA analysis of the NP dataset resulted in the identification and characterization of 22 distinct cell clusters, shedding light on the cellular heterogeneity within the sample. Based on the analysis, seven cell types were annotated as the primary markers for the 22 clusters. Each cell type was assigned a unique color, allowing for clear visualization and classification. The identified cell types included fibroblasts, dendritic cells (DC), induced pluripotent stem cells (iPS cells), osteoblasts, epithelial cells, macrophages, neutrophils, and bone marrow progenitor cells (BM Prog) (Figure 8(a)). A heatmap was generated to display the expression patterns of the top five marker genes for each cell cluster. Each row represented a gene, while each column represented a cell cluster. The top five marker genes were highlighted in red, indicating their significance in distinguishing the clusters (Figure 8(b)). UMAP and t-SNE dimensional plots were employed to visualize the 22 different cell clusters identified in the previous step (Figure 8(c) and (d)). These plots provided a comprehensive view of the cellular landscape, enabling the observation of distinct clusters and their relationships. In both the UMAP and t-SNE projections, the annotation of seven cells as primary labels for the clusters reinforces their significance in defining and characterizing the distinct cell populations, thereby providing valuable insights into the cellular composition of the NP dataset (Figure 8(e) and (f)). Overall, through integration and characterization of the 22 cell clusters, this scRNA analysis revealed primary marker genes and unveiled the cellular heterogeneity within the NP dataset, paving the way for further investigations and understanding of the underlying biology.

Figure 8.

Integration and characterization of 22 cell clusters reveals primary marker genes and cellular heterogeneity. (a) 7 cell types were annotated as the primary markers of the 22 clusters. Different cell types were colored with unique colors. The NP dataset cells may be classified into 22 clusters, which include fibroblasts, dendritic cell (DC), induced pluripotent stem cells (iPS cells), osteoblasts, epithelial cells, macrophages, neutrophils and Bone marrow progenitor cell (BM Prog). (b) Heatmap of the first 5 marker genes per cell cluster, with each row representing the gene and each column representing the cell cluster. The top 5 marker genes were labeled in red color. (c) UMAP and (d) tSNE dimensional plot helped to visualize 22 different cell clusters (b). 7 cells were annotated as primary labels of the clusters in (e) UMAP and (f) tSNE projection.

Differential expression and cellular distribution of hub mito genes and DRGs in neuronal and non-neuronal cell types in NP samples

In NP samples, we observed downregulation of two key DRGs (ACTB and SLC2A8) and five core Mito genes (TOMM5, PDHB, MRPL39, SIRT3 and MRPS15) in neuronal cells compared to non-neuronal cells (Figure 9(a)). This analysis aimed to identify genes that exhibit distinct expression patterns in neurons compared to non-neuronal cells. UMAP plots were used to visualize the expression distribution of these five central Mito genes and two DRGs in the clustered dataset (Figure 9(b)). The accessibility of each gene region was marked in red to provide insights into their spatial expression patterns across different cell types. Additionally, the expression levels of these key genes were investigated in different cell types under both SNI and SHAM controls (Figure 9(c)). This analysis aimed to understand the potential roles of these genes in NP development. By examining the differential expression and cellular distribution of these central Mito genes and DRGs in neuronal and non-neuronal cell types, this analysis provides valuable insights into their potential involvement in the pathological development of NP.

Figure 9.

Differential expression and cellular distribution of hub mito genes and DRGs in neuronal and non-neuronal cell types in NP samples. (a) Differential expression of DRGs and 5 hub Mito genes in neurons cells versus other cell types in NP samples. (b) UMAP shows the expression distribution of 5 hub mito genes and 2 DRGs in the clustered set. The accessibility of each gene region is indicated in red. (c) 5 hub mito genes and 2 DRGs expression levels in different cell types under SHAM and SNI controls.

Communication analysis in NP cells

Our analysis demonstrates the results of cell trajectory analysis and intercellular communication analysis of NP cells, focusing on their different modes of differentiation. Pseudotemporal analysis was performed using single-slice analysis to infer cell trajectories. The direction of cell differentiation is from dark blue on the right to light blue on the left (Figure 10(a)). This analysis helps in understanding the progression and differentiation patterns of NP cells over time. The trajectory map shows that DRGs (Figure 10(b)) and Mito (Figure 10(c)) genes manifest high expression in the middle and late stages of cell differentiation. It provides insights into how the expression of these genes changes as cells transition through different states during differentiation. The single-cell trajectory of neurons, defined by their clusters in NP (Figure 10(d)). This analysis helps in identifying the developmental trajectory and diversity of neuron subtypes in NP. A comparison of interaction weights and strengths between key cells in NP and SHAM control is provided (Figure 10(e)). This analysis reveals differences in intercellular communication networks between NP and control conditions. The heatmap shows a smaller number of potential ligand-receptor pairs in neuronal cells under NP than under SHAM conditions (Figure 10(f)). This may be related to the fact that key genes are less expressed in NP than in SHAM. Furthermore, we conducted a detailed analysis of specific ligand-receptor correlations in cellular communication under NP and SHAM conditions. Our findings revealed that neuronal cells exhibited strong interactions with macrophages primarily through the APP-CD74 ligand-receptor pathway. Additionally, neurons were found to engage in communication with Oligodendrocytes, Endothelial cells and neurons predominantly via the CADM1-CADM1 ligand-receptor pathway (Figure 11(a)). This investigation of ligand-receptor correlations provides valuable insights into the potential signaling interactions between cells in the context of NP and control conditions.

Figure 10.

Cell trajectory and intercellular communication analysis of distinct differentiation patterns in NP cells. (a) Pseudo-time analysis using monocle analysis to infer cell trajectory. Trajectory plot illustrating the expression patterns of (b) DRGs and (c) mito across cell states. (d) Single-cell trajectory of neurons defined by cluster from NP, with each point corresponding to a single cell and colors indicating neuron subtypes. (e) Comparison of interaction weights and strengths between key cells in NP and SHAM conditions. (f) Heatmaps displaying the number of potential ligand-receptor pairs in key cells in NP and SHAM conditions.

Figure 11.

Correlations between ligands and receptors in cellular communication: a comparison between NP and SHAM. (a) Analysis of ligand and receptor correlations for cellular communication in NP and SHAM conditions, respectively.

Discussion

Neuropathic pain presents a significant global health challenge, impacting around 6.9–10% of the global population, necessitating novel treatment approaches.^1,35 Mitochondria play a crucial role in regulating cell death, and various forms of cell death have been linked to neuroimmune mechanisms in neuropathic pain.^11,14 However, the involvement of mitochondrial dysfunction and disulfidptosis in neuropathic pain remains unclear. Further research is needed to elucidate their roles in this condition.

In our study, we used RNA-seq data and single-cell data from mouse models of NP to explore the mechanistic link between disulfidptosis and mitochondria in NP. Through DEGs analysis and WGCNA analysis, we identified 171 hub genes overlapping with Mito-related genes and 4 hub genes overlapping with DRGs. We further enriched 171 key MITO genes and 33 DRGs. Mito was mainly enriched in Metabolic pathways and mitochondrion-related pathways, and DRGs were mainly enriched in Regulation of actin cytoskeleton and transmembrane transport-related pathways. This result validates the reliability of the gene set definition. Previous research has demonstrated the presence of filamentous microtubules in neural spines and the significance of microtubule-modulating drugs in neuropathic pain (NP). ³⁶ Additionally, a study has proposed a connection between trigeminal neuropathic pain episodes and transmembrane proteins. ³⁷ Consequently, investigating the involvement of mitochondria and disulfidptosis in the mechanisms underlying the development of NP and exploring targeted therapeutic approaches offer valuable insights.

Machine Learning Validation of Mito’s Critical First 6 Genes, Including TOMM5, PDHB, MRPL39, SIRT3, MRPS15, and CPT1C. These genes are less expressed in NP than in the SHAM1 group and have good dispute value. Additionally, ACTB (belongs to DRGs) was highly expressed in NP, while G6PD and SLC2A8 showed low expression compared to the SHAM group. Translocase of Outer Mitochondrial Membrane 5 (TOMM5) is a protein-coding gene that encodes a subunit of the translocase of the outer mitochondrial membrane (TOM) complex. This complex is responsible for importing nuclear-encoded mitochondrial proteins into the mitochondria. ³⁸ Pyruvate dehydrogenase E1 subunit beta (PDHB) is an isoform of the pyruvate dehydrogenase complex and is a gene that encodes a lipid acylating protein in mitochondria. ³⁹ PDHB expression is upregulated early in peripheral nerve injury and promotes peripheral axon regeneration by regulating energy supply. ⁴⁰ Study identifies PDHB as a biomarker for predicting tumor immune response and immunotherapy. ⁴¹ PDHB was less expressed in the NP group compared to the SHAM1 group, indicating mitochondrial dysregulation, an important clinical therapeutic target. Mitochondrial Ribosomal Protein L39 (MRPL39) and Mitochondrial Ribosomal Protein S15 (MRPS15) are genes encoding ribosomal proteins localized in the mitochondria, involved in protein synthesis within the organelle.^42,43 SIRT3, a member of the sirtuin protein family encoded by the SIRT3 gene, is involved in the regulation of mitochondrial metabolism and the response to oxidative stress.^44,45 Research suggests that SIRT3-mediated CypD-K166 deacetylation alleviates neuropathic pain by improving mitochondrial dysfunction and suppressing oxidative stress. ⁴⁶ SIRT3 also attenuates neuropathic pain through a mechanism that deacetylates FoxO3a in the dorsal horn of the spinal cord in a diabetic model of rats. ⁴⁷ In addition, the reduction of spinal SIRT3 leading to neuropathic pain provides new insights into the mechanisms of chronic opioid therapy for neuropathic pain in HIV patients. ²⁰ CPT1C is a gene that encodes an isoform of the enzyme carnitine palmitoyltransferase 1, which is involved in the transport of long-chain fatty acids into the mitochondria for β-oxidation. ⁴⁸ In disulfidptosis, G6PD is involved in cellular energy production and protection against oxidative stress, while SLC2A8 is responsible for glucose transport and regulation of glucose metabolism.^49,50 Significantly, preliminary studies indicate that ACTB holds promise as a candidate biomarker and a potential therapeutic target for NP. ⁵¹ These genes are strongly associated with mitochondrial function and disulfidptosis, providing a solid foundation for investigating NP and exploring the mechanisms involving mitochondrial dysfunction and disulfidptosis-induced cell death.

Immune infiltration analyses highlighted the consistent effects of disulfidptosis and mitochondrial dysfunction on immune cell infiltration. Subtype analysis further exemplifying the association between the key gene PDHB and G6PD and immune cell type provided insights. Single-cell analysis revealed that the key genes were mainly expressed in Neurons and were able to be specific markers for cell clusters. The key genes were lowly expressed in the SNI group compared to the SHAM group, which validated the cellular level reason for the low expression of key genes in NP samples.

The proposed time-series analysis revealed that neuronal cells communicated closely with non-neuronal cells, and there was a clear trend of down-regulation in the occurrence of SNI, further suggesting that NP occurrence produces neuronal disruptive effects. Our findings suggest that neuronal cells exhibit strong interactions with macrophages mainly through the APP-CD74 ligand receptor pathway. It has been noted that the significant increase in CD74 expression in tubular epithelial cells is dependent on neighboring infiltrating macrophages, which is where necrotic injury occurs as a result of proximal tubular cell repair failure. ⁵² The study revealed that the depletion of macrophages in the dorsal root ganglia reduces the development and maintenance of nerve injury-induced mechanical hypersensitivity, which is a key factor in the onset and persistence of NP. ⁵³ In the immune infiltration, macrophages were highly expressed in NP, negatively correlated with TOMM5 and PDHB expression, and highly expressed and co-expressed with G6PD in the C2 subgroup analysis. These findings provide valuable insights into the immune response and molecular mechanisms underlying NP, with a specific focus on the involvement of macrophages in the APP-CD74 pathway, particularly in association with neurons. Additionally, these findings highlight the associations between macrophages and specific genes, such as TOMM5, PDHB, and G6PD, emphasizing their potential relevance in the context of NP.

There are certain limitations to our study. This paper primarily relies on NP mouse samples sourced from public databases, lacking experimental verification of the reliability of key genes and mechanisms. Additionally, further validation of the current results using human samples is necessary. Despite these limitations, our findings highlight the coexistence of factors mitochondrial dysfunction and disulfidptosis in the mechanism of NP occurrence, identify key gene targets, and lay the foundation for innovative therapeutic strategies for NP. Future research should delve into the mechanisms of action of these targets and conduct preclinical trials to evaluate targeted therapies.

In summary, our research has provided insights into the molecular mechanisms underlying mitochondrial dysfunction and disulfidptosis in NP. We have identified key genes and their involvement in disease progression, as well as the role of the macrophage-associated APP-CD74 pathway in neuronal cell communication. These findings contribute to our understanding of strategies for NP remission. While further validation studies are required, our results represent a significant advancement in NP research and hold promise for future therapeutic developments.