Ferroptosis-related genes participate in the microglia-induced neuroinflammation of spinal cord injury via NF-κB signaling: evidence from integrated single-cell and spatial transcriptomic analysis

Abstract

Background

Ferroptosis and immune responses are critical pathological events in spinal cord injury (SCI), whereas relative molecular and cellular mechanisms remain unclear.

Methods

Micro-array datasets (GSE45006, GSE69334), RNA sequencing (RNA-seq) dataset (GSE151371), spatial transcriptome datasets (GSE214349, GSE184369), and single cell RNA sequencing (scRNA-seq) datasets (GSE162610, GSE226286) were available from the Gene Expression Omnibus (GEO) database. Through weighted gene co-expression network analysis and differential expression analysis in GSE45006, we identified differentially expressed time- and immune-related genes (DETIRGs) associated with chronic SCI and differentially expressed ferroptosis- and immune-related genes (DEFIRGs), which were validated in GSE151371. Protein–protein interaction and microRNA-mRNA-transcription factor regulatory networks were constructed based on Search Tool for the Retrieval of Interacting Genes (STRING) and NetworkAnalyst, respectively, which were validated by chromatin immunoprecipitation followed by sequencing (ChIP-seq). Cell subclusters and unique features of microglia in SCI were identified by single-cell transcriptomic analysis, which were validated in GSE226286. Spatial expression patterns of DETIRGs and DEFIRGs were validated in brain injury (GSE214349) and SCI (GSE184369). Potential mechanisms underlying neuronal regeneration by neurotrophin-3 (NT3)-chitosan were revealed by transcriptomic analyses in GSE69334. Immune- and ferroptosis-related mechanisms of nanolayered double hydroxide loaded with NT3 (LDH-NT3) were investigated in vivo and in vitro.

Results

GBP2, TEC, UNC93B1, PLXNC1, NFATC1, IL10RB, and TLR8 were DETIRGs represented chronic SCI-specific genes and peripheral blood biomarkers. NFKB1 may regulate expression of CYBB and HMOX1 in a unique subcluster of M1 microglia within the middle SCI lesion, establishing links between microglial ferroptosis and neuroinflammation. Reduced inflammatory responses and microglial ferroptosis were potential effects of NT3-chitosan or LDH-NT3 on neuronal regeneration.

Conclusions

A novel subcluster of microglia exhibiting M1 polarization and ferroptosis phenotype was involved in SCI. These microglia may trigger neuroinflammation and induce neuronal degeneration within the middle site of SCI, which might be inhibited by NT3-chitosan or LDH-NT3.

Introduction

Spinal cord injury (SCI) can not only cause severe loss of sensory, motor, and autonomic function below the lesions [1], but also lead to huge financial burdens and despairing psychological stress on patients and society [2]. No effective clinical therapy exists for SCI because of the complicated pathological processes and molecular mechanisms. SCI can be classified into primary and secondary stages [3]. Primary injury is resulted from indirect or direct mechanical damages to the spinal cord, which may cause necrosis and demyelination of neuronal cells and axons. Secondary injury is a multifaceted cascade reaction after the primary injury, including spinal cord tissue edema, local hemorrhage, inflammatory response, oxygen free-radical injury, ionic imbalance, mitochondrial dysfunction, glutamate toxicity, necrosis, and programmed cell death (PCD).

Differing from other types of PCD (apoptosis, autophagy, and pyroptosis), ferroptosis is featured by dependence on ferric ion and accumulation of reactive oxygen species (ROS) [4]. Ferroptosis is involved in various central nervous system (CNS) degenerative diseases and trauma [5]. Ferroptosis is also critical in the pathological processes of SCI [6]. Hence, it is important to investigate the molecular and cellular mechanisms of ferroptosis in SCI. Besides, immune responses play vital roles in the secondary SCI, which can not only facilitate repair but also aggravate the injuries [7]. After SCI, biomarkers may release into the peripheral blood via the broken blood-spinal cord barrier (BSCB) [8]. During secondary SCI, various types of immune cells and cytokines in peripheral blood may enter the lesions via the broken BSCB, which play significant parts within the local microenvironmental immune responses [9]. Moreover, the immune states of peripheral blood can reflect the severe degree of SCI [10]. Hence, a pressing need arises for exploring the candidate signatures of blood and central immune features in SCI patients.

Microglia are the resident immune cells and core sensors of danger signals in the CNS. After stimulation, resident M2 type (anti-inflammatory) microglia may polarize to the M1 type (pro-inflammatory) [11], which can secrete pro-inflammatory factors [12]. However, ferroptosis and immune responses in microglia after SCI, as well as potential regulatory mechanisms remain a black box. Transcription factors (TFs) and microRNAs (miRNAs) are important regulators of differential gene expression in various diseases [13], especially in the inflammatory responses and ferroptosis post SCI [14, 15]. Constructing TF-miRNA-mRNA networks may reveal the potential regulatory mechanisms underlying the inflammatory responses and ferroptosis post SCI.

Here, we conducted a comprehensive bioinformatic analysis based on micro-array datasets (GSE45006, GSE69334), RNA sequencing (RNA-seq) dataset (GSE151371), spatial transcriptome datasets (GSE214349, GSE184369), and single cell RNA sequencing (scRNA-seq) datasets (GSE162610, GSE226286) to investigate potential biomarkers, miRNA-TF-gene regulatory networks, and unique cell populations involved in the ferroptosis and immune responses post SCI. According to our previous research, neurotrophin-3 (NT3)-loaded chitosan enabled an excellent micro-environment to promote new neurogenesis, nerve repair and growth in animals [8]. Further, nanolayered double hydroxide (LDH) is a candidate for constructing a favorable microenvironment for SCI recovery [16]. Here, the effects of LDH loaded with neurotrophic factor 3 (LDH-NT3) on microglial activation, immune responses, and ferroptosis in SCI were investigated in vivo and in vitro.

Methods

Data collection

This study was endorsed by the Ethics Committee of the Shanghai Tongji Hospital affiliated to Tongji University. The micro-array expression matrix GSE45006 of SCI rats was obtained from the GEO (gene expression omnibus) database (https://www.ncbi.nlm.nih.gov/geo/) [17]. Sequencing platform of GSE45006 is GPL1355 [Rat230_2] Affymetrix Rat Genome 230 2.0 Array (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL1355). In GSE45006, the thoracic spinal cord (T7) of rats was injured using aneurysm clip impact-compression injury model, the epicenter of the injured spinal cord was isolated, and RNA was obtained and processed. The global gene expression in acute, subacute and chronic stages of moderate to severe injury to the rat spinal cord was analyzed using micro-array gene chip approaches. The probes from the raw data in GSE45006 were converted to gene names, and the raw data were subjected to standardized pretreatment using the limma package in R. The Microarray data from GSE45006 were transformed to human genes using biomaRt package (Version: 2.46.0) in a certain homologous degree with “getLDS” function [18]. We divided the injured spinal cords into acute SCI group (1 day after SCI), subacute SCI group (7 days after SCI), and chronic SCI group (56 days after SCI), with four rats in each group, and the 4 intact spinal cords were defined as sham group (laminectomy from T6 to t8 without SCI). Besides, another micro-array data including 3 uninjured rats and 164 SCI rats were downloaded from the GEO database (Accession number: GSE69334) [8].

RNA-seq data of Homo sapiens peripheral blood mononuclear cells (PBMCs) form 58 samples (10 healthy individuals without CNS pathology history (healthy control, HC), 10 non-CNS trauma patients (trauma control, TC), and 38 patients with acute traumatic SCI)) were obtained using the Illumina’s HiSeq4000 sequencer, which were available in GEO database (Accession number: GSE151371) [19], which were used as the validation set. Samples with incomplete clinical information were excluded from this study.

ScRNA-seq data of normal and injured spinal cord at 1 dpi, 3 dpi, and 7 dpi in wild-type mice were obtained from GEO database (Accession number: GSE162610). Besides, another scRNA-seq data of whole brain cells from healthy control and 3-round microglia depletion-repopulation (3 × DR) mice were collected from GEO database (Accession number: GSE226286), which were used as the validation set [20]. Spatial transcriptome sequencing data of wild-type C57/BL6 mouse brain tissues were available in GEO database (Accession number: GSE214349) [21]. Besides, the spatial transcriptome sequencing data (Accession number: GSE184369) of wild-type C57/BL6 mouse spinal cord post SCI was used as the validation set [22].

Ferroptosis-related genes (FRGs) were determined using the PathCards database of ferroptosis pathways (http://pathcards.genecards.org/) [23]. 1793 immune-related genes (IRGs) were obtained from the ImmPort database (https://immport.niaid.nih.gov) [24]. The IRGs contained cytokines, cytokine receptors, and genes correlated with the B-cell antigen receptor and T-cell receptor signaling pathways, natural killer cell cytotoxicity, and antigen processing and presentation pathways. Lists of TFs were obtained from the Cistrome database (http://cistrome.org/) [25]. The flow chart and mechanism diagram of this study was constructed using Biorender (https://biorender.com/).

Micro-array data processing

Original micro-array data were read using affy package in R, followed by robust multi-array average (RMA) background correction, standardization, probe specific background correction, and summarizing probe set values in a expression measure [26]. Appropriate batch effects correction was conducted using the sva package [27].

RNA-seq data processing

Quality control (QC) of raw fastq files was performed using fastQC software and Trimmatic software to obtain clean data (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Then, hierarchical indexing for spliced alignment of transcripts 2 (HISAT2, Version: Ensemble GRCh38) was used in the mapping of clean data to the human reference genome, providing output alignments in Sequence Alignment/Map (SAM) format. The SAM files were converted to binary alignment/map (BAM) format using SAMtools (Version: 1.8) [28] and were analyzed by FeatureCounts [29] to estimate transcript count. The count was expressed as fragments per kilobase of transcript per million mapped reads (FPKM), and it was then converted to transcripts per million (TPM). Eventually, a gene expression matrix with gene symbol as row and sample ID as column was generated.

Differential expression analysis

According to previous studies, log₂ (fold change (FC)) was used to evaluate the significance and the change in expression of genes between different types of samples, and |log₂ FC|> 1.0 and false discovery rate (FDR) < 0.05 were considered significant in differential expression analysis [30]. In the micro-array datasets (GSE45006, GSE69334) and RNA-seq dataset (GSE151371), the “limma” R package was used to identify DEGs [31], and |log₂ (FC)|> 1.0 and FDR < 0.05 were the screening criteria for DEG selection [30]. By taking the intersection of DEGs and FRGs, differentially expressed ferroptosis-related genes (DEFRGs) were annotated. Likewise, differentially expressed immune-related genes (DEIRGs), were identified by taking the intersection of DEGs and IRGs. Differentially expressed ferroptosis- and immune-related genes (DEFIRGs) were identified by taking the intersection of DEFRGs and IRGs. Then, differentially expressed transcription factors (DETFs) were identified based on the same threshold using the “limma” R package. Heatmaps and volcano plots were constructed using the R packages “pheatmap” and “ggplot2” [32].

In the scRNA-seq dataset of GSE162610, the “FindAllMarkers” function with "wilcox" method was used to identify DEGs between normal control samples and SCI samples from top 2,000 variable genes. Only genes with |log₂Fold Change (FC)|> 0.5 and FDR < 0.05 were identified as DEGs.

Quantification of hallmark signaling pathways

50 hallmark pathways were obtained from the Molecular Signatures Database (MSigDB, version 7.1) (https://www.gsea-msigdb.org/gsea/msigdb/index.jsp) [33], and the absolute quantification of the pathways was investigated using gene set variation analysis (GSVA). Further, differentially activated pathways between uninjured samples and SCI samples were identified [33].

Functional enrichment analysis

Gene Oncology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were conducted to investigate the potential downstream signaling pathways of identified DEGs, DEFRGs, and DEIRGs [34]. Biological processes where these genes were mostly enriched were determined using the "cluster Profiler" R package [35] based on the cutoff value of P < 0.01 and FDR < 0.05.

Immune cell infiltration analysis

The Cell-type Identification By Estimating Relative Subsets Of RNA Transcripts (CIBERSORT) algorithm [36] was employed to investigate potential immune cell infiltrates in normal spinal cord and injured spinal cord. The gene expression matrix of GSE45006 was uploaded to CIBERSORT database (https://cibersort.stanford.edu/) to purify cellular subtype-specific gene expression. Only patient genes and immune cells with P < 0.05 were considered as statistically significant. Moreover, nonparametric tests were performed to determine differential infiltration degrees of 22 types of immunocytes. Eventually, correlation analysis was performed for predicting potential communication and mutual recognition among the 22 types of immune cells.

Weighted gene co-expression network analysis (WGCNA)

WGCNA (v1.61; https://cran.r-project.org/web/packages/WGCNA/index.html) is a powerful tool for identifying gene clusters or modules and constructing gene co-expression networks [37]. Therefore, the WGCNA integration algorithm R (v3.4.1) was utilized to explore highly relevant native gene modules associated with time points after SCI. There were ≥ 10 cutoff genes, cutting height = 0.85, Z‐ score ≥ 5, and stability-related stability correlation P ≤ 0.05 in our work. The connection of genes between the two was utilized to calculate the dataset of GSE45006, and genes with the correlation coefficient < 0.5 were excluded. The conservation status of WGCNA modules and the related features were investigated. Genes in the module that was positively associated with the time points post SCI were defined as time-related genes (TRGs). Further, differentially expressed time- and immune-related genes (DETIRGs) were identified by taking the intersection of DEIRGs and TRGs.

Analysis of TF-gene co-regulation and TF-miRNA interactions

The DEFIRGs and DETIRGs were mapped to the NetworkAnalyst (https://www.networkanalyst.ca/) to analyze gene regulation mechanisms in Homo sapiens [38]. Protein–Protein Interaction (PPI) analysis was conducted using the Search Tool for the Retrieval of Interacting Genes (STRING) database (https://string-db.org/). Gene-miRNA interaction analysis was conducted based on experimentally validated miRNA-gene interaction data from TarBase [39] and miRTarBase [40]. TF-gene interaction analysis was conducted based on the data from the ENCODE database (Encyclopedia of DNA Elements, https://www.encodeproject.org/). TF and hub gene target data were obtained from the ENCODE chromatin immunoprecipitation followed by sequencing (ChIP-seq) data. The predicted regulatory potential score < 1 and the peak intensity signal < 500 was set as cutoff value [41]. Finally, TF-miRNA-gene networks of the DEFIRGs and DETIRGs were constructed based on the NetworkAnalyst.

Prediction of potential small-molecule drugs

The Connectivity Map (CMAP) (https://clue.io) [42] is a commonly used database for investigating the correlations between gene expression patterns and small-molecule drugs. We uploaded the gene lists of DEFIRGs and DETIRGs of SCI to the CMAP database and chose the “Latest” version for further analysis. Compounds with negative scores might contribute to the therapy of the disease [43].

Spatial transcriptomic analysis

Intraventricular hemorrhage (IVH) is one of the most frequent complications after various acute traumatic brain injuries (TBIs) secondary to traumatic SCI. However, the pathophysiological mechanisms of how IVH stimulations post SCI lead to neurological dysfunction remain largely unknown [21]. Therefore, we analyzed the spatial transcriptome data (accession number: GSE214349) of wild-type C57/BL6 mouse brains from the needle-puncture group with blood injection on the seventh days (Injured group), no-surgery group (control group), and the needle-puncture group without blood injection (sham group), based on Spatial Transcript Omics DataBase (STOmics DB, https://db.cngb.org/stomics/) [44]. Furthermore, based on Comprehensive Repository of Spatial Transcriptomics (CROST) database (https://ngdc.cncb.ac.cn/crost) [45], we analyzed the spatial transcriptome data (accession number: GSE184369) of wild-type C57/BL6 mouse spinal cord post SCI to reveal the spatial transcriptomic atlases of the injured spinal cord and demonstrate the spatial expression features of identified DEFIRGs. Bayesian-based Space-clustering for spatial transcriptomics (BayesSpace) tool was used to perform the dimensionality reduction and clustering analysis [46]. Cell type annotation was performed using the Spatially-resolved Transcriptomics cell type identification tool (SPOTlight) [47]. CellChat algorithm was used to explore cell–cell communication patterns across different cell types [48].

Clinical correlation analysis

Based on the micro-array dataset GSE69334, integrated transcriptome analyses of spinal cord segments were performed at the lesion site, as well as caudal to the lesion regions and immediately rostral, over a period of 90 days (examined at 1, 3, 10, 20, 30, 60, and 90 dpi) post SCI. American spinal injury association (ASIA) grading system was used primarily to evaluate the neurological status of SCI through motor and sensory exams [49]. Here, the earliest ASIA Impairment Scale (AIS) grades were obtained between day 3 and 10 d after SCI to avoid producing substantially unstable AIS grades due to various factors such as strong sedation and spinal shock during the first 48 h. Likewise, we also avoided unnecessary biological and technical variance from 11 to 56 days post SCI. Specifically, 5 patients without AIS grades during the time window were removed from this study. Using nonparametric tests, differentially gene expression patterns across SCI samples with different AIS grades and lesion regions were identified.

SCI had profound impacts on circulating white blood cells (WBCs), causing alterations of WBC phenotypes and peripheral inflammatory responses in a dynamic cascade which may represent the pathological features of the evolving spinal lesions. Based on the peripheral blood RNA-seq dataset of GSE151371, DEG analysis was conducted across different groups (TC, HC, and SCI) as discussed previously [30]. Further, it should be emphasized that the larger the |log₂FC|, the greater the difference in gene expression. Therefore, the DEG analysis of GSE151371 was performed under two different conditions (|log₂FC|> 1 and FDR < 0.05 [30]; |log₂FC|> 2 and FDR < 0.05 [50]). After taking the intersection of different groups (HC vs. SCI, HC vs. TC, and TC vs. SCI), we extracted the intersecting DEGs for the subsequent analysis. Furthermore, univariate logistic regression was performed to identify the predictive value of the intersecting DEGs for SCI. Then, multivariate logistic regression was conducted to identify independent predictors of SCI, based on the intersecting DEGs with P < 0.05 in univariate logistic regression. Bidirectional stepwise regression was performed for constructing the regression model. For nomogram model construction, all samples in GSE151371 were randomly assigned to training set (n = 40) and validation set (n = 18) at a ratio of 7:3. A nomogram model was constructed using R software (version 4.3.1). Receiver operating characteristic (ROC) curves were developed and area under the curve (AUC) value was calculated to assess the model discrimination. Decision curves were constructed for evaluating clinical value based on the validation set.

External validation at single cell level

To validate the DEFIRG expression patterns of different microglia subclusters in CNS, scRNA-seq data of whole brain cells from control (2180 cells) and 3-round microglia depletion-repopulation (3 × DR) mice (1313 cells) were collected from GEO database (accession number: GSE226286) (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE226286) [20]. In dataset GSE226286, C57BL/6 J mice were administered a PLX5622-formulated diet for 2 weeks to deplete microglia in CNS, while the control mice were fed with an AIN-76A control diet (CD) for 3 weeks for repopulation [20]. Specifically, following three rounds of depletion–repopulation, every microglia has proliferated more than twenty times. Moreover, 3 × DR can shorten the telomere length of microglia, which showed an aged-like phenotype [51]. Because the cumulative number of microglia divisions is about 20 cycles in 2-month-old mice [52], every microglia may proliferate around 40 times, which was close to the Hayflick limit. Finally, the senescence index of 3 × DR microglia was 3.6-fold that of normal control microglia, echoing the aged-like phenotype of 3 × DR microglia [20]. A previous analysis revealed that 3 × DR microglia (aged -like microglia) exhibited ferroptosis phenotype in the CNS [53].

The dataset of GSE226286 was analyzed with Seurat 3.2.2 [54]. Using RunPCA with VariableFeatures, the top 20 principal components (PCs) of the dataset GSE226286 were identified, which were utilized for clustering and t-distributed stochastic neighbor embedding (TSNE)/uniform manifold approximation and projection (UMAP) visualization. For transcriptomic analysis of microglia, clusters exhibiting high expression of Cd79b, Camp, S100a4, S100a9, and Nkg7 were considered as non-microglia clusters and excluded from microglial scRNA-seq analysis. To determine different cell types in CNS, various cell-type-specific or enriched markers were used as previously described [55]. Specifically, these marker genes involved, but were not limited to, Hexb, Aif1, P2ry12, and Tmem119 for microglia; Aldh1l1, S100b, and Slc1a2 for astrocytes; Acta2 for vascular smooth muscle cells; Vtn for pericytes; Nkg7 for natural killer (NK) cells; Abcg2 for ECs; Dcx and Tubb3 for neurons; S100a9 for neutrophils; Mog, Mbp, and Mag for oligodendrocytes; Pdgfra for oligodendrocyte progenitor cells (OPCs); Ccdc153 for ependymocytes; and Ccr2 and Plac8 for monocytes. The expression level of our identified DEFIRGs in different cell types of 3 × DR and control mice were visualized using feature plots and violin plots.

Animal experiments and ethics statement

All experimental processes were approved by and conducted according to the standards of the Animal Welfare Committees of Tongji University in Shanghai, China. A total of 256 male and female C57BL/6 mice (aged 8–9 weeks) weighing 18–20 g (Jisijie Laboratory Animal Center, Shanghai, China) were utilized in this study. All animal operations were conducted with general anesthesia with isoflurane (1–2%). Animals were daily monitored to prevent infection, abnormal wound healing, and loss of weight. To assist micturition, animals’ bladders were squeezed twice a day, till the animals were executed.

Surgical procedure

Based on our previous work, completely transected and excised mouse models were constructed with female C57BL/6 mice [56]. The T8–9 dorsal spinal cord segments were exposed via laminectomy using a dissecting microscope (Nikon/Zeiss), which were then transected in the T8–9 level, and a 2.0 mm spinal segment was resected. The blade was repeatedly scraped along the ventral surface of spinal canal, and the residual fibers within the lesions were removed by aspiration, to make sure the complete transection was performed laterally and ventrally. Muscle and skin of the mice were closed in layers.

Bulk RNA-seq

The 1 cm long spinal cord segments located in lesion core, an average weight of 20 mg, at 15 min (n = 3), 1 day (n = 3), 3 days (n = 6), 1 week (n = 6), 2 weeks (n = 5), 4 weeks (n = 6), and 6 weeks (n = 4) post SCI, as well as spinal cord segments from uninjured mice (n = 6), were homogenized using 1 mL TRIzol (Invitrogen). Then, mRNA was obtained and purified using RNeasy Mini Kit (QIAGEN) as per instructions. The RNA qualities were evaluated using Agilent 2100. Subsequently, 1–3 μg qualified RNA was used to construct Illumina V2 RNA-seq library, followed by Hiseq 2000 SE50bp sequencing.

Extraction of spinal cord single cells

The 1 cm long spinal cord segments located in the lesion core were obtained at different time points post injury, which were dissociated in papain in Hibernate A-Ca (HE-Ca, BrainBits, LLC) at a concentration of 2 mg/ml for 10 min at 30 °C. Digested spinal cord tissues were triturated 10 times using fire-polished Pasteur pipettes. The single cell suspension was filtered using a 40 μm cell strainer, which was then centrifuged at 200g for 5 min at room temperature. The pellet was washed using PBS containing 0.04% bovine serum albumin. Viabilities of cells were evaluated using trypan blue staining. The single cells were diluted to a concentration of 1 × 10⁶/mL in PBS with 0.04% BSA. Libraries were constructed using the Chromium v3 Single Cell 3′ Library and Gel Bead Kit (10 × Genomics) as per instructions. QC was conducted using a DNA 1000 chip (Agilent Technologies), followed by sequencing using Illumina NovaSeq 6000.

Single cell RNA-seq data processing

After the procedures of 10 × Genomics Chromium (https://www.10xgenomics.com/instruments/chromium-x-series) [57], preprocessing of scRNA-seq data was performed. When data demultiplexing was completed, the sequencing results were categorized into two pair-ended reads fastq files. Then, these fastq files were trimmed to eliminate template switch oligo (TSO) sequence and polyA tail sequence. Moreover, clean reads were aligned with the GRCh38 (Version: 100) genome assembly and were quantified using Cell Ranger Software (Version 1.0.0) (http://10xgenomics.com/).

The quantified gene expression matrices were investigated using Seurat pipeline (Version: 3.2.2) [54]. Cells with more than 100,000 transcripts expression and less than 10% mitochondrial gene mapped, together with genes which expressed in more than two single cells were extracted for further analysis. Following the QC processes, single cells were merged into one Seurat object with the “IntegrateData” function, which was scaled and standardized using “ScaleData” function. In addition, the “vst” method was used to identify the top 2000 variable genes. For reducing dimensionality, principal component analysis (PCA) was performed, in which the first twenty PCs were merged as an input for UMAP analysis. Similarly, all control samples were also merged into one single Seurat object and processed by all analysis procedures aforementioned.

Annotation and characterization of cell types

DEGs in identified sub clusters were utilized as potential references to determine cell type of each unsupervised cluster, which were integrated with well accepted cell surface biomarkers from CellMarker [58]. Further, cellular feature plots, dot plots, violin plots, and heatmaps were generated for illustrating the markers of each cell type using SCANPY (Version: 1.7.1) and Seurat R package (Version: 3.2.2) with Python 3.6 [59].

Analysis of cellular communication

The core cellular communication patterns and ligand-receptor (LR) pairs across multiple cell types in SCI were identified using iTALK R package (Version: 0.1.0) (https://github.com/Coolgenome/iTALK) [60]. Based on the scarcity of Rattus norvegicus resources in present mainstream cellular communication algorithms, the top 2,000 variable genes were transformed to Homo sapiens genes using biomaRt R package (Version: 2.46.0) [18]. Subsequently, the standardized expression matrix of the genes was merged in the iTALK object using the “rawParse” function. Eventually, the top 200 LR pairs were identified using the edge bundle R package (https://github.com/garthtarr/edgebundleR) [61].

Trajectory analysis

Monocle package in R (Version: 2.18.0) was utilized to analyze specific cellular lineage pseudo-time trajectories and translational relationships across different cell clusters [62]. Genes utilized for ordering trajectories were obtained from the top DEGs across different cell clusters using the differential GeneTest function with default parameters, which were identified as cell differentiation-related genes (DRGs). Importantly, the DRGs may regulate different cellular chemotaxis fates. Following the selection of DRGs, single cells were projected into a low-dimensional space reduced from expression profiles and ordered in pseudo-time trajectories using Monocle′s reversed graph embedding algorithm. Pseudo-time of each cell type at different time points was identified using “reduceDimension” function with “DDRTress” method.

Preparation of LDH-NT3

As previously described, LDH was prepared with modifications [63], and NT3 was loaded into LDH using ion exchange intercalation with modifications [64]. The C57BL/6 mice were randomly divided into three groups. In LDH group (n = 76), LDH was implanted into the lesion sites. In LDH-NT3 treatment group (n = 65), LDH-NT3 was implanted into the lesion sites. In control group (n = 76), no treatment following the surgery was performed. Following the surgery, all mice were placed on sawdust bedding and kept warm. Water and food were provided ad libitum.

Spinal cord microglia culture

Microglial culture was prepared from the spinal cords of 1 day post-natal Sprague–Dawley rats [65]. Blood vessels and meninges were removed from the spinal cords, which was then minced and dissociated enzymatically using 0.25% Trypsin–EDTA at 37 °C for 20 min. Tissues were mechanically triturated in Dulbecco Modified Eagle Medium/Ham’s F12 (DMEM/F12; Gibco) containing 2% Penicillin/Streptomycin (P/S; Gibco) and 10% Fetal Bovine Serum (FBS; Gibco), which were then plated on the poly-L-lysine coated plates. Following cell culture for 3 weeks, the microglia were extracted by mild trypsinization for 20 min to 50 min (0.25% Trypsin–EDTA was diluted in 1:3). Subsequently, DMEM/F12 containing 2% P/S was added to the microglia [65]. Further, cells were tagged with PE-conjugated anti-CD11b + antibody (BD Biosciences, San Jose, CA, USA) followed by anti-PE antibody conjugated to the magnetic bead. Magnetically tagged CD11b + cells were isolated by utilizing MS columns as per the Miltenyi MACS protocol, which were then annotated as microglia [66]. Further, the isolate microglia were cultured in DMEM/F12 containing 2% P/S in 5% CO2 at 37 °C. The third passages of cells were utilized for the subsequent experimental processes.

Cell counting kit-8 (CCK8) assay

The cytotoxicity of blank LDH and LDH-NT3 to microglia was evaluated using the CCK8 assay as per the manufacturer's instructions. Firstly, the microglia were seeded into 24-well plates overnight. Wells with no cells but medium only were used as the blank control. Secondly, cells were exposed to different treatments. Specifically, the LDH: NT3 ratio of LDH-NT3 was set as 50,000:1 (w/w). The LDH concentrations for every group ranged from 10 μg/mL to 40 μg/mL, while the NT3 concentrations for every group ranged from 0.2 pg/mL to 0.8 pg/mL. 24, 48, or 72 h later, microglia were incubated using 50 µL CCK-8 solution at 37 °C for additional 3 h. The absorbance was detected at 450 nm using a microplate reader (Thermo Fisher Scientific Inc.). As compared with the blank control group (100%), the viability of microglia in these groups was evaluated, respectively.

Extraction of mRNA and RT-qPCR analysis

Total mRNA from microglia was extracted using TRIzol (Takara, Japan), which were reverse transcribed to cDNA with a PrimeScript RT Reagent Kit (Takara, Japan) and subjected to RT-qPCR analysis by using TB Green Premix Ex Taq (Takara, Japan). Glyceraldehyde-3-phosphate dehydrogenase (Gapdh) was used as a reference gene. The relative expression levels of target mRNAs in the experimental groups were quantified by Ct(2^−ΔΔCt) method and normalized to internal control Gapdh. The primers used for RT-qPCR were listed in supplementary Table S1.

Statistical analysis

In the present study, only two-tailed P value < 0.05 and FDR < 0.05 was suggested to be statistically significant. The statistical analysis processes were conducted with R version 4.3.1 software (Institute for Statistics and Mathematics, Vienna, Austria; www.r-project.org), Python version 3.6 software (https://www.python.org/), and Strawberry Perl version 5.30.0.1 software (https://www.perl.org/). For descriptive statistics, the mean ± standard deviation was used for continuous variables with a normal distribution. Moreover, the median (range) was used only for continuous variables with an abnormal distribution.

Results

DEG identification

Based on the micro-array data of GSE45006, 1,380 DEGs were identified between SCI and sham groups in rats (P < 0.05). The heatmap showed expression levels of the top 100 most significant DEGs between SCI group and sham group in rats (Fig. 1a). In the volcano plot, blue dots represented downregulated DEGs, and the red dots indicated upregulated genes (Fig. 1b). After transforming the DEGs to Homo sapiens genes, a total of 1,140 DEGs recorded in Homo sapiens were extracted for the subsequent analysis.

Fig. 1

Differentially expressed genes (DEGs) identification and immune cell infiltrates in SCI and sham samples. a Heatmap of the top 100 DEGs between SCI and sham samples. b Volcano plot of differentially expressed genes (DEGs) between SCI and sham samples. c Gene ontology (GO) enrichment analysis of DEGs in biological processes, cellular components, and molecular functions; Kyoto encyclopedia of genes and genomes (KEGG) analysis of DEGs. d Distribution of 22 kinds of immune cells in SCI and sham groups. e Infiltration degree of 22 kinds of immune cells in SCI and sham groups. f Correlation diagram between 22 kinds of immune cells

Functional enrichment analysis of DEGs

GO and KEGG functional enrichment analyses were performed based on the identified DEGs of GSE45006. The most involved processes or components in GO terms included establishment of localization, transport, regulation of biological quality (biological process), vesicle, plasma membrane part (cell component), transporter activity, inorganic cation transmembrane transporter activity, and channel activity (molecular function). Additionally, the most significantly enriched KEGG pathways mainly included cAMP signaling pathway, osteoclast differentiation, and retrograde endocannabinoid signaling (Fig. 1c).

Analysis of immune cell infiltration

Based on the micro-array data of GSE45006, the immune cell infiltration analysis was performed using the CIBERSORT algorithm. The stacked bar graphic indicated the main infiltrating immune cell types in SCI tissues were the resting dendritic cells and macrophages M2 (Fig. 1d). Noticeably, infiltration of macrophages M2 was increased in SCI tissue as compared with normal tissue, although not statistically significant (P > 0.05). Further, the infiltration of resting dendritic cells significantly increased in SCI tissues, whereas the infiltration of macrophages M0, resting mast cells, activated CD4 + memory T cells, and was reduced (all P < 0.05, Fig. 1e). The correlation heatmap showed a strong positive correlation between resting mast cells and macrophages M0 (R = 0.86), whereas resting dendritic cells exhibited a significantly negative correlation with resting CD4 + memory T cells (R = − 0.64, Fig. 1f). It suggested that varying infiltration degrees of immunocytes might function as a vital role in the neuroinflammation and immune responses in SCI.

WGCNA analysis to identify SCI time-related genes

Based on the micro-array data of GSE45006, β = 1 was selected as the suitable soft threshold to construct the scale-free network, and the scale-free topology model fit degree was − 0.91 (Fig. 2a). Additionally, the average connectivity was 3308.52 (Fig. 2b). Modules were hierarchically clustered according to the topological overlap matrix, and similar modules were merged in the clustering tree (Fig. 2c). Four co-expression modules were identified, which showed obviously different expression characteristics (Fig. 2d). It indicated that the purple module containing 519 genes was positively correlated with time points after SCI (R = 0.47, P < 0.05), while the light cyan module containing 508 genes was negatively correlated with time points after SCI (R = − 0.59, P < 0.05) (Fig. 2e). We selected the purple module to further investigate the potential upregulated genes in the chronic stage of SCI, which were defined as TRGs in this study.

Fig. 2

Weighted gene co-expression network analysis to identify differentially expressed time- and immune-related genes (DETIRGs). a Estimation of the scale independence index of the 1–30 soft threshold power (β = 1). b Determination of the mean connectivity of the 1–30 soft threshold power. c Identification and construction of gene co-expression modules that were illustrated as different colors in hierarchical clustering. d Clustering dendrogram of module feature genes. e Correlations between different modules and time points post SCI. Red represented a positive correlation, and green represented a negative correlation. f Venn diagram showing the 13 DETIRGs obtained by intersecting. g Gene ontology (GO) enrichment analysis of DETIRGs in biological processes, cellular components, and molecular functions. h Kyoto encyclopedia of genes and genomes (KEGG) analysis of DETIRGs

Functional enrichment analysis of DEIRGs

By taking the intersection of 1140 DEGs from GSE45006 and 1793 IRGs from ImmPort database, 119 DEIRGs were determined. Subsequently, we conducted GO and KEGG enrichment analyses specific to the DEIRGs (Fig. S1a). The most statistically significant GO terms were closely associated with regulation of response to stimulus, immune system process, cell surface receptor signaling pathway (biological process, Fig. S1b), extracellular region, plasma membrane part, whole membrane (cell component, Fig. S1c), signaling receptor binding, molecular transducer activity, and signaling receptor activity (molecular function, Fig. S1d). Besides, the top 3 most enriched KEGG pathways of the DEIRGs were cytokine-cytokine receptor interaction, natural killer cell mediated cytotoxicity, and tuberculosis (Fig. S1e).

Functional enrichment analysis of differentially expressed time- and immune-related genes (DETIRGs)

After overlapping 1140 DEGs from GSE45006, 519 TRGs identified by WGCNA, and 1793 IRGs from ImmPort database to take the intersection, 13 genes were identified as DETIRGs (Fig. 2f). Subsequently, GO and KEGG enrichment analyses were conducted based on the DETIRGs. The differential expression patterns of DETIRGs were illustrated by box plots between the rat spinal tissues of sham group and SCI group at different time points after SCI. Inevitably, all the 13 genes (T-lymphocyte activation antigen CD86 (CD86), C-X-C chemokine receptor type 4 (CXCR4), guanylate-binding protein 2 (GBP2), tyrosine-protein kinase Tec (TEC), nuclear factor of activated T-cells, cytoplasmic 1 (NFATC1), oxysterols receptor LXR-alpha (NR1H3), Plexin-C1 (PLXNC1), interleukin-10 receptor subunit beta (IL10RB), integrin alpha-L (ITGAL), macrophage metalloelastase (MMP12), toll-like receptor 7 (TLR7), protein unc-93 homolog B1 (UNC93B1), and Toll-like receptor 8 (TLR8)) were significantly upregulated in chronic stage of SCI (8 weeks after SCI), as compared with the sham group (all P < 0.05, Fig. S2).

The most statistically significant GO terms (Fig. 2g) included immune system process, immune response, regulation of immune response (biological process), receptor complex, lytic vacuole, lysosome (cell component), signaling receptor activity, molecular transducer activity, transmembrane signaling receptor activity (molecular function). The most significantly enriched KEGG pathways of the DETIRGs included toll-like receptor signaling pathway, human cytomegalovirus infection, and intestinal immune network for IgA production (Fig. 2h).

Construction of PPI and miRNA-mRNA-TF network of DETIRGs and prediction of potential inhibitors

The PPI network of the DETIRGs was constructed using the STRING database, which consisted of 23 nodes and 66 edges, and the average node degree was 5.74 (Fig. 3a). The hub genes within the PPI network can be grouped into 3 different clusters using k-means clustering analysis. Cluster 1 included TLR7, TLR8, NR1H3, UNC93B1, GBP2, and TEC, which were enriched in immune response-activating signal transduction (GO biological process), Golgi membrane (GO cellular component), and signaling receptor activity (GO molecular function, Fig. 3b); cluster 2 included CCL21 (C–C motif chemokine ligand 21), CD4 (T-cell surface glycoprotein CD4), CXCR4, ICAM1 (intercellular adhesion molecule 1), ICAM3 (intercellular adhesion molecule 3), IL10RA (Interleukin-10 receptor subunit alpha), IL10RB, IL22RA1 (interleukin-22 receptor subunit alpha-1), IL26 (Interleukin-26), ITGAL, ITIH4 (inter-alpha-trypsin inhibitor heavy chain H4), MMP12, NFATC1, and CD80 (T-lymphocyte activation antigen CD80), which were enriched in cell surface receptor signaling pathway (GO biological process), extracellular region (GO cellular component), and signaling receptor activity (GO molecular function, Fig. 3c); cluster 3 included SEMA7A (semaphorin-7A) and PLXNC1.

Fig. 3

Construction of protein–protein interaction (PPI) and miRNA-DETIRG-transcription factor (TF) regulatory network. a PPI network of DETIRGs was constructed based on the STRING database. b Gene ontology (GO) enrichment analysis of hub genes from cluster 1 in biological processes, cellular components, and molecular functions. c GO enrichment analysis of hub genes from cluster 2 in biological processes, cellular components, and molecular functions. d MiRNA-mRNA-TF regulatory network of DETIRGs was constructed based on the Networkanalyst database. e Connectivity Map (CMAP) analysis identified the top 10 most specific small-molecular inhibitors for DETIRGs. f Mode of action (MoA) and structural formula of the top 3 most specific drugs (staurosporine, zibotentan, and AKT-inhibitor-1–2). g ChIP-seq validation demonstrated the binding relationships between TFs and DETIRGs in the regulatory network

A miRNA-mRNA-TF interaction network of DETIRGs was constructed using NetworkAnalyst database, which included 36 nodes, 62 edges, and 13 seeds. CXCR4 might be regulated by TFs including STAT2 (signal transducer and activator of transcription 2), E2F1 (E2F transcription factor 1), USF1 (upstream transcription factor 1), MYC (MYC proto-oncogene), and STAT3 (signal transducer and activator of transcription 3). Additionally, TFs such as SPI1 (Spi-1 proto-oncogene) and SP1 (Sp1 transcription factor) may regulate the expression of DETIRGs including TRL7 and MMP12. CXCR4 and MMP12 were interacted with has-miR-9, while TLR8, TLR7, and CXCR4 were interacted with has-miR-9. (Fig. 3d). These potential interaction pairs took up the core position in the miRNA-mRNA-TF interaction network.

The 13 DETIRGs were imported into the CMAP web tool. After the signature query, we identified the top 10 small-molecular drugs exhibiting the highest negative enrichment scores with potential therapeutic effects via inhibiting DETIRGs (Fig. 3e, supplementary Table S2). Specifically, a τ score of − 90 indicated that only 10% of reference perturbations exhibited stronger connectivity to the query. Among the top 3 potential inhibitors (Fig. 3f), staurosporine with a median τ score less than −90 might be the most promising drug reversing abnormal immune responses by targeting DETIRGs.

ChIP-seq validation of TF-DETIRG interactions

To further validate the interaction mechanisms between DETIRGs and TFs, we utilized ChIP-Seq data reporting DNA binding sites for the TFs from Cistrome Data Browser (http://cistrome.org/db). We derived DETIRG lists which were potentially interacted with the predicted TFs by the identification that the TFs bind to these genes’ promoters. Importantly, we identified that, STAT2, STAT3, E2F1, USF1, and MYC had various DNA binding peaks in the promoter regions of CXCR4; SPI1 and SP1 had various DNA binding peaks in the promoter regions of TLR7 and MMP12 (Fig. 3g). This further demonstrated that these TFs may bind to and modulate DETIRG expression.

Identification of DEFRGs in SCI

After taking the intersection of 1140 DEGs from GSE45006 and 64 FRGs from PathCards database, 9 FRGs were identified as DEFRGs (Fig. S3a). The differences in expression of DEFRGs were shown by box plots between the rat spinal tissues of sham group and SCI group, which were statistically significant at P < 0.05. CYBB (cytochrome b-245 beta chain), HMOX1 (heme oxygenase 1), HSPB1 (heat shock protein family B (small) member 1), CPT1A (carnitine palmitoyltransferase 1A), SAT1 (spermidine/spermine N1-acetyltransferase 1), SLC39A8 (solute carrier family 39 member 8), and GCH1 (GTP cyclohydrolase 1) were differentially elevated in the rat spinal tissues of SCI group, in contrast to FDFT1 (farnesyl-diphosphate farnesyltransferase 1) and ACSL6 (acyl-CoA synthetase long chain family member 6) (all P < 0.05, Fig. S3b).

GO terms where the DEFRGs mostly enriched included cellular response to oxidative stress, response to oxidative stress, cofactor metabolic process (biological process, Fig. S3c), endomembrane system, endoplasmic reticulum part, endoplasmic reticulum (cell component, Fig. S3d), protein dimerization activity, identical protein binding, and protein homodimerization activity (molecular function, Fig. S3e). The KEGG pathways where the DEFRGs mostly enriched were ferroptosis, metabolic pathways, and fatty acid degradation (Fig. S3f).

Construction of PPI and miRNA-mRNA-TF network of differentially expressed ferroptosis- and immune-related genes (DEFIRGs) and prediction of potential inhibitors

By taking the intersection of 1,140 DEGs from GSE45006, 64 FRGs from PathCards database, and 1,793 IRGs from ImmPort database, 2 DEFIRGs (CYBB and HMOX1) were selected and defined as DEFIRGs (Fig. 4a). Furthermore, we identified that the expression of CYBB and HMOX1 was significantly upregulated after SCI, with persistent overexpression during the acute (1 day), subacute (7–14 days), and chronic (8 weeks) phases of SCI, compared to the sham group (all P < 0.05, Fig. S4a, b).

Fig. 4

Construction of protein–protein interaction (PPI) and miRNA-DEFIRG-transcription factor (TF) regulatory network. a Venn diagram showing the 2 DEFIRGs obtained by intersecting; PPI network of DEFIRGs was constructed based on the STRING database. b Gene ontology (GO) enrichment analysis of hub genes from cluster 1 in biological processes, cellular components, and molecular functions. c GO enrichment analysis of hub genes from cluster 2 in biological processes, cellular components, and molecular functions. d MiRNA-mRNA-TF regulatory network of DEFIRGs was constructed based on the Networkanalyst database. e Connectivity Map (CMAP) analysis identified the top 10 most specific small-molecular inhibitors for DEFIRGs. f Mode of action (MoA) and structural formula of the top 3 most specific drugs (benzanthrone, piperacetazine, and westcort). g ChIP-seq validation demonstrated the binding relationships between TFs and DEFIRGs in the regulatory network

The PPI network of the DEFIRGs was constructed using the STRING database to identify potential interactions across the overlapping genes that contained 12 nodes and 52 edges, and the average node degree was 8.67. The hub genes involved in the PPI network can be divided into 3 different clusters by k-means clustering analysis. Cluster 1 included CYBB, neutrophil cytosol factor 1 (NCF1), neutrophil cytosol factor 2 (NCF2), neutrophil cytosol factor 4 (NCF4), and Ras-related C3 botulinum toxin substrate 2 (RAC2), which were enriched in respiratory burst (GO biological process), phagocytic vesicle (GO cellular component), and anion binding (GO molecular function, Fig. 4b); cluster 2 included cytochrome b-245 light chain (CYBA), NADPH oxidase 1 (NOX1), NADPH oxidase 4 (NOX4), NADPH oxidase activator 1 (NOXA4), and NADPH oxidase organizer 1 (NOXO1), which were enriched in reactive oxygen species metabolic process (GO biological process), plasma membrane part (GO cellular component), and oxidoreductase activity, acting on NAD(P)H, oxygen as receptor (GO molecular function, Fig. 4c); cluster 3 included HMOX1 and biliverdin reductase A (BLVRA), which were enriched in heme catabolic process (GO biological process).

The NetworkAnalyst database was utilized to construct the miRNA-mRNA-TF interaction network of DEFIRGs, which involved 35 nodes, 35 edges, and 2 seeds. As illustrated in Fig. 4d, CYBB and HMOX1 could be regulated by two common TFs: nuclear factor kappa B subunit 1 (NFKB1) and SPI1. HMOX1 was interacted with three miRNAs, including has-miR-873, has-miR-196a, and has-miR-16. CYBB was interacted with six miRNAs, including has-miR-186, has-miR-138, has-miR-369-3p, has-miR-551b, has-miR425, and has-miR-653. Nevertheless, these results need further experimental validation.

The 2 DEFIRGs were imported into the CMAP web tool, and we identified the top 10 small-molecular drugs with the highest negative enrichment scores with potential therapeutic effects via inhibiting expression of DEFIRGs (Fig. 4e, supplementary Table S3). Among the top 3 potential inhibitors (Fig. 4f), benzanthrone with a median τ score less than −90 might be the most promising drug reversing abnormal immune responses by targeting DEFIRGs.

ChIP-seq validation of TF-DEFIRG interactions

We derived DEFIRG lists from Cistrome Data Browser (http://cistrome.org/db), which had potential interactions with TFs in the regulatory network. Various DNA binding peaks of NFKB1 and SPI1 were identified in the promoter regions of CYBB and HMOX1 (Fig. 4g), which further demonstrated the interaction mechanisms of core TF-DEFIRG pairs within the regulatory network.

Gene expression landscapes of cells in SCI revealed by scRNA-Seq analysis in GSE162610

Following strict QC, scRNA-seq dataset GSE162610 of 66,428 cells collected from 7 SCI samples and 3 uninjured samples were integrated and analyzed (Fig. 5a). As shown in the UMAP plots, a total of 39 clusters and 15 cell types, namely astrocyte, dendritic cell, Div-myeloid cell, endothelial cell, ependymal cell, fibroblast, lymphocyte, macrophage, microglia, monocyte, neuron, neutrophil, oligodendrocyte, OPC, and pericyte were clearly identified based on gene expression profiles and specific cell markers [55] (Fig. 5b, c). Further, the bar plot showed the average count and specific proportion of the cell types within spinal cord tissues of 7 SCI samples and 3 uninjured samples (Fig. 5d). Specifically, the percentage of the absolute number of microglia, macrophage, and dendritic cell in SCI samples was obviously higher than that in uninjured samples.

Fig. 5

Comprehensive analysis of Single-cell sequencing (scRNA-seq) and spatial transcriptome analysis. a Dimensionality reduction analysis of SCI samples (1 day, 3 days, and 7 days after SCI) and uninjured samples. b Unsupervised clustering identified 39 clusters based on 66,428 single cells. c Unsupervised clustering and dimension reduction analysis annotated 15 main cell types (astrocyte, dendritic cell, myeloid-derived cells, endothelial cell, ependymal cell, fibroblast, lymphocyte, macrophage, microglia, monocyte, neuron, neutrophil, oligodendrocyte, oligodendrocyte progenitor cell, and pericyte) based on known cell markers. d Bar plot showed the cell number and percentage of the 15 cell types in all SCI and uninjured samples. e Venn diagram showing the intersection of microglial markers, DEFIRGs, DETIRGs, DEFIRG-related TFs, and DETIRG-related TFs. f Cellular feature plots of DEFIRGs in microglia. g Cellular feature plots of overlapping DETIRGs in microglia. h Cellular feature plots of overlapping TFs in microglia. i Normalization and clustering of the cell types with similar gene expression. j Hematoxylin and eosin (H&E)-stained sections and spatial gene expression in traumatic brain injury mice models

The most highly variable genes of each cell type were identified as typical cell markers, and we took the intersection of the cell markers, DEFIRGs, DETIRGs, DEFIRG-related TFs, and DETIRG-related TFs. Importantly, CYBB (47.9%), HMOX1 (75.5%), NFKB1 (48.5%), SPI1 (87.3%), CD86 (54.5%), CXCR4 (49.0%), IL10RB (74.2%), TLR7 (34.7%), and UNC93B1 (88.6%) were mainly expressed in microglia; STAT3 (47.5%) and ITGAL (41.4%) were mainly expressed in neutrophils; SP1 (29.2%) and NR1H3 (26.0%) were mainly expressed in pericytes; GBP2 (26.0%) was mainly expressed in endothelial cells; STAT2 (31.3%) was mainly expressed in astrocytes; MYC (26.6%) was mainly expressed in fibroblasts (supplementary Table S4). The Venn plot showed the intersection relationship across microglia markers, DEFIRGs, DEFIRG-related TFs, DETIRGs, and DETIRG-related TFs (Fig. 5e).

As significant immune cells, microglia undergo a series of ferroptosis-related changes after SCI that increase the susceptibility to neural degeneration [67]. Therefore, we focused on microglia to investigate ferroptosis and immune microenvironment after SCI, and validate our previous findings. A total of 19,687 microglia annotated by canonical markers and dimension reduction were extracted for the subsequent analysis. Cellular feature plots showed high expression level and clear cell distribution of DEFIRGs (CYBB and HMOX1) (Fig. 5f), DETIRGs (CXCR4, TLR7, UNC93B1, CD86, IL10RB) (Fig. 5g), and TFs of DEFIRGs and DETIRGs (NFKB1 and SPI1) (Fig. 5h) in microglia.

Hematoxylin and eosin (H&E)-stained biopsies and spatial gene expression patterns

Based on the spatial transcriptomic dataset of GSE214349, a total of 19 original frozen brain sections of five IVH-stimulated mice were used as the input data. After performing spatial transcriptomic sequencing analysis, the H&E-stained brain sections were evaluated to examine spatial gene expression of overlapping DEFIRGs and DETIRGs in secondary TBI. The SME algorithm [12] was used for normalizing and clustering the cell types with similar gene expression characteristics (Fig. 5i). Specifically, five main cell types were annotated, including astrocyte (blue cluster), mural cell (yellow cluster), neuron (green cluster), microglia (red cluster), and Schwann cell (purple cluster). Spatial characteristics of expression alterations of DEFIRGs (CYBB and HMOX1), DETIRGs (CXCR4, TLR7, UNC93B1, CD86, IL10RB), and TFs of DEFIRGs and DETIRGs (NFKB1 and SPI1) were revealed in the biopsies from injured group, control group, and sham group. It highlighted that expression levels of DEFIRGs in periventricular tissues, the epicenters or middle sites of secondary injury lesions, was increased after IVH stimulation, compared to the control group and sham group (Fig. 5j).

Based on the spatial transcriptomic dataset of GSE184369, we investigated the cell type distribution and specific gene expression post SCI at the spatial transcriptome level. A total of 41,26 spots were identified in 3 SCI mice samples (including 1,390, 1,409, and 1,327 spots in VISDS000257, VISDS000262, and VISDS000271, respectively; Fig. S5a). Dimensionality reduction and clustering identified 8 cell types based on known marker genes. We compared the spots of microglia with Hmox1 and Cybb, and it showed that the spatial expression position of Hmox1 and Cybb was similar to that of the microglia. We examined the spatial location of microglia and gene expression of Hmox1 and Cybb in the red box to highlight co-localization features (Fig. S5a). Cell–cell communication analysis revealed the complex interactions between microglia and other cell types, in which microglia occupied the central site with high strength of interactions (Fig. S5b).

Microglia lineages’ characteristics in SCI

Based on the scRNA-seq dataset GSE162610, the microglia were merged into a Seurat object, on which the cell type annotation was conducted and these cells were segregated into 6 major cell subclusters (Microglia 1–6). Figure 6a showed cell proportion of different microglia subclusters in SCI samples and uninjured samples. Importantly, the expression of CYBB, HMOX1, CXCR4, TLR7, UNC93B1, CD86, IL10RB, NFKB1, and SPI1 was higher in SCI samples than that in uninjured samples. The UMAP plot and bar plot showed the total quantity distribution of the 6 microglia subclusters (Fig. 6b, c). Specifically, Microglia2 and Microglia3 showed relatively higher proportion in SCI patients than that in uninjured samples, while Microglia1 was lower. Cleveland's dot plot showed the expression level of several key marker genes (JUN (Jun proto-oncogene) [68], CTSD (cathepsin D) [69], IFITM3 (interferon induced transmembrane protein 3) [70], and C1QL3 (complement component 1, q subcomponent-like 3)) [71] in 6 microglia subclusters, and relative cellular proportion across different samples (Fig. 6d). It showed that high-level expression of CTSD was significantly apparent in all microglia clusters, while C1QL3 was lowly expressed. Importantly, CTSD was critical in protection from neurodegeneration and maintaining cell homeostasis. After SCI, lysosomal membrane permeabilization allowed leakage of CTSD in cytosol, suppressing lysosomal activity and inhibiting autophagy flux, which may induce neuronal cell death [72].

Fig. 6

Identification of the microglia clusters, cellular communication analysis, cell cycle analysis, and pseudo-time analysis. a Dimensionality reduction analysis of microglia from SCI samples (1 day, 3 days, and 7 days after SCI) and uninjured samples. b Unsupervised clustering identified 6 microglia subclusters (Microglia1-6). c Bar plot showed the cell number and percentage of the 6 microglia subclusters in all SCI and uninjured samples. d The Cleveland's dot plot showed the expression levels of 4 classical marker genes of microglia; the stacked bar plot displayed the proportion of 6 microglia clusters in different samples (7 SCI samples and 3 uninjured samples). e Heatmaps of the most variable markers of 6 microglia subclusters. f The cell cross-talk network showed the intersected cellular communication among 6 microglia clusters. g Cell cycle analysis of the 6 microglia subclusters. h The heat map showed the biological features and expression levels of 31 signaling pathways in 6 microglia clusters by gene set variation analysis (GSVA). i The cell trajectory state plot illustrated the developmental trajectory of microglia with pseudo-time, 5 cell states shown in different colors, and all cell subpopulations illustrated in the trajectory. j Violin diagrams showing the expression level of DEFIRGs and potential upstream TFs in 6 microglia subclusters. k Violin diagrams showing the expression level of M1 polarization markers in 6 microglia subclusters

The heatmaps showed expression of the top 30 variable genes and the top 5 marker genes in each microglia subcluster respectively, and these microglia subpopulations exhibited decidedly different gene expression patterns (Fig. 6e).

There was obviously strong intercellular communication across 6 microglia clusters, in which Microglia2 occupied the central site of cellular communication and mostly sent out and received signals to other microglia subclusters, and SPP1 (secreted phosphoprotein 1) and APOE (apolipoprotein E) were the key cellular communication genes (Fig. 6f). Cell cycle analysis indicated that these microglia subclusters mainly composed of cells in the S and G2M phases with enriched expression of cell cycle-associated signatures, which showed high activity of proliferation (Fig. 6g).

To further identify the potential functional activities of these microglia subclusters, activation level of GO, KEGG, and hallmark pathways in each microglia cluster was shown in the heatmap. Importantly, various pro-inflammatory pathways, such as GO positive regulation of inflammatory response, GO inflammatory response, hallmark inflammatory response, GO cytokine production involved in inflammatory response, GO inflammatory cell apoptotic process, and reactome inflammasomes were notably activated in Microglia3 subcluster, which may form a pro-inflammatory immune microenvironment after SCI, and thus can aggregate neuroinflammation and induce neural apoptosis and degeneration (Fig. 6h).

Pseudo-time analyses recapitulated microglia differentiation trajectories

Based on the scRNA-seq dataset GSE162610, pseudo-time analysis was performed to reconstruct cell-developmental trajectories of all cells of the 6 microglia subclusters using Monocle 2 in R. Importantly, these microglia bifurcated into decidedly different developmental states at each time point, and branch point 1 and point 2 clearly separated the developmental trajectory into 5 cell states (Fig. 6i). The differentiation of microglia subclusters followed four directions, and the overall differentiation of microglia occurred from cell state 1, 5, and 2 to cell state 3 and 4 (Fig. 6j). Microglia1-3 subclusters showed early occurrence, which were present on two branches and differentiated later via taking the second branch point as a critical turning point. Importantly, the DEFIRGs (CYBB, HMOX1) and DEFIRG-related TFs (NFKB1, SPI1) were highly-expressed and co-expressed in Microglia3 subcluster, indicating that ferroptosis occurred in this unique microglia subcluster after SCI. In addition, it revealed that activating expression of ferroptosis-related biomarkers co-occurred with the high expression of M1-associated biomarkers in Microglia3 subcluster (Fig. 6k): secreted markers (IL-1β, IFN-γ), surface markers (CD86), and intracellular/transcription factor markers (STAT1) [73]. Furthermore, the expression trends of CYBB and HMOX1 were in consistent with the developmental trajectory of Microglia3, which showed a persistent activation status during the pathological processes of SCI. Taken together, these findings indicated that ferroptosis and M1 polarization were both promoted in Microglia3 subcluster, which exhibited pro-inflammatory capacities. According to previous evidence [74], we speculated that, high resistance of M1 polarized microglia to ferroptosis may contribute to the presence of Microglia3 in SCI mice. Given the strong pro-inflammatory effects of ferroptotic cells [75], the modulation of Microglia3 may be of vital importance to control inflammatory and immune responses post SCI.

External validation of the differential expression of DEFIRGs and DETIRGs on transcriptome level

In dataset GSE151371, a total of 8,987 DEGs were identified between 38 SCI patients and 10 HC samples, and 3,327 DEGs were identified between 38 SCI patients and 10 TC samples (Fig. S6a), based on the cutoff value of |log₂FC|> 1, P < 0.05. By taking the intersection of these DEGs (SCI vs. HC, SCI vs. TC), DEFIRGs, and DETIRGs, we identified that 7 DETIRGs (GBP2, TEC, UNC93B1, PLXNC1, NFATC1, IL10RB, and TLR8) and 2 DEFIRGs (CYBB, HMOX1) (Fig. S6b) were still significantly differentially expressed in Homo sapiens between SCI patients and HCs. Specifically, NFATC1, IL10RB, and TLR8 were the DETIRGs that were differentially expressed in SCI samples, as compared with HCs and TCs. Further, based on the threshold of |log₂FC|> 2 and P < 0.05, 590 DEGs were identified in the blood samples between the HC and SCI groups, and 198 DEGs were identified in blood samples between the SCI and TC groups (Fig. S6c). However, no intersection was observed between these DEGs and DETIRGs or DEFIRGs (Fig. S6d). According to previous studies, the threshold of |log₂FC|> 1 and P < 0.05 has been widely used as a general principle to identify DEGs [76,77,78]. Therefore, the external validation in GSE151371 further confirmed that 7 DETIRGs (GBP2, TEC, UNC93B1, PLXNC1, NFATC1, IL10RB, and TLR8) and 2 DEFIRGs (CYBB, HMOX1) were differentially expressed in human blood samples, potentially functioning as peripheral blood biomarkers post SCI.

Univariate logistic regression showed that (Fig. S6e), UNC93B1 (odds ratio (OR) = 10.890, 95% confidence interval (CI)  2.620–45.270, P < 0.001), NFATC1 (OR = 0.080, 95% CI 0.020–0.330, P < 0.001), PLXNC1 (OR = 10.720, 95% CI 2.360–48.710, P = 0.002), TLR8 (OR 95.150, 95% CI 9.310–972.670, P < 0.001), IL10RB (OR = 20.840, 95% CI 4.180–103.920, P < 0.001), and HMOX1 (OR = 15.450, 95% CI 3.190–74.700, P < 0.001) were predictors for SCI. Multivariate logistic regression showed that (Fig. S6f), TLR8 (OR = 89.170, 95% CI 4.690–1695.630, P = 0.003) and HMOX1 (OR = 29.880, 95% CI 1.760–506.960, P = 0.019) were independent predictors for SCI. A nomogram including the two independent predictors was established (Fig. S6g). In this model, the points on the axes of independent risk factors (peripheral blood RNA-seq expression level of TLR8 and HMOX1) corresponded to different risk scores. The AUC of the ROC curve of the nomogram established by the training set was 0.963 (95% CI 0.919–1.000), and the AUC of the ROC curve of the validation set was 0.963 (95% CI 0.860–1.000) (Fig. S6h). Decision curves (Fig. S6i) demonstrated excellent performance of the nomogram over a wide range of threshold probabilities. Besides, expression level of the top 100 most significant DEGs between SCI patients and HC samples (Fig. S6j) and the top 100 DEGs between SCI patients and TC samples (Fig. S6k) were shown in the heatmaps, respectively.

External validation of DEFIRG expression in aged-like microglia exhibiting ferroptosis phenotype at single-cell resolution

Based on the validation scRNA-seq dataset GSE226286, the UMAP plot showed that, gene expression features of 3 × DR microglia were different from the control group (Fig. S7a). Importantly, we identified Cybb was highly-expressed in 3 × DR microglia compared to microglia in control mice, whereas Hmox1 was highly expressed in both 3 × DR microglia and control microglia (Fig. S7b). We further investigated the transcriptome between brain tissues of 3 × DR and age-matched control mice at single-cell resolution. A total of 16 cell types (17 populations) were identified in brain tissues of 3 × DR and control mice (Fig. S7c, d). Specifically, 3 × DR microglia and control microglia located in different clouds, whereas other cell populations did not exhibit such separation patterns, further demonstrating that microglia were the major cell population influenced by 3 × DR treatment directly. Importantly, high expression of Cybb and Hmox1 was only identified in 3 × DR microglia/macrophages, and the unique expression patterns were further demonstrated in the violin plots (Fig. S7e).

Region-specific ferroptosis-related alterations after SCI and potential anti-ferroptosis activity of NT3 by inhibiting DEFIRG expression

Based on the micro-array dataset GSE69334, a total of 1409 DEGs were identified across different spinal cord regions (caudal, middle, rostral, and uninjured sites) in 164 SCI and 3 normal Rattus norvegicus with |log₂FC|> 1, FDR < 0.05 set as cutoff value, which were shown in the heatmap and volcano plot (Fig. 7a). Further, 23 DETFs were identified across different spinal cord regions based on the threshold of |log₂FC|> 1 and FDR < 0.05, which were shown in the heatmap and volcano plot (Fig. 7b). Activation levels of 50 hallmark pathways across different spinal cord regions were illustrated in the heatmap (|log₂FC|> 0.1, FDR < 0.05) (Fig. 7c). Specifically, hallmark interferon alpha response, hallmark interferon gamma response, hallmark TNFA signaling via NFKB, hallmark IL2 STAT5 signaling, hallmark IL6 JAK STAT3 signaling, hallmark apoptosis, and hallmark inflammatory response were significantly activated in the middle region of injured spinal cord, indicating that pro-inflammatory and neural degeneration events may happen in the middle regions. Besides, hallmark angiogenesis was also significantly activated in the middle SCI regions, which was needed for proper tissue repair, indicating dramatic injury repair and reconstitution events occurred in middle regions after SCI. Importantly, NT3-chitosan treatment promoted angiogenesis, and it also inhibited the inflammatory pathways, which were widely accepted as the key to repressing secondary injuries post SCI.

Fig. 7

Region-specific ferroptosis-related and immune-related alterations after SCI based on GSE693334. a Heat map and volcano plot of differentially expressed genes (DEGs) in different clusters (Time, Treatment, and Location). b Heat map and volcano plot of 23 differentially expressed transcription factors (DETFs) in different clusters (Time, Treatment, and Location). c Heat map illustrated the differential regulation of 50 hallmark signaling pathways in different clusters (Time, Treatment, and Location). d Box plots showing differential expression level of differentially expressed ferroptosis- and immune-related genes (DEFIRGs) and related transcription factor (TF) across different SCI lesion sites and uninjured spinal tissues; box plot showing differential expression level of Hmox1 (heme oxygenase 1) across different treatment groups and uninjured group. e Mechanism map of our scientific hypothesis. NFKB1 (nuclear factor kappa B subunit 1) may regulate expression of ferroptosis-related markers CYBB (cytochrome b-245 beta chain) and HMOX1 that were specifically expressed in a unique subcluster of microglia with M1 polarization, which may trigger neuroinflammation and induce neuronal degeneration within the middle site of SCI

Interestingly, we identified that the DEFIRGs (Cybb, Hmox1) and their TF (Nfkb1) were all significantly upregulated in the middle regions of injured spinal cord, compared with caudal regions and rostral regions of injured spinal cord as well as uninjured spinal cord (all P < 0.05, Fig. 7d). Furthermore, Hmox1 expression was significantly inhibited after NT3-chitosan treatment, as compared with lesion control group and uninjured group (all P < 0.05, Fig. 7d). We speculated that NT3-chitosan may suppressed ferroptosis of microglia, and the anti-inflammatory effects elicited by NT3-chitosan may also facilitate the nerve regeneration.

Taken together, upon stimulation, the resident anti-inflammatory or M2 phenotype microglia may differentiate into a unique microglia subcluster exhibiting pro-inflammatory phenotype that secretes the pro-inflammatory cytokines (IL-1β, FN-γ), which can form a neuroinflammatory microenvironment. Besides, within the unique microglia subpopulation, TF (NFKB1) may promote the expression of ferroptosis markers (CYBB, HMOX1), and enhance ferroptosis and various pro-inflammatory pathways in this microglia subcluster, and this effect could be reversed by NT3-chitosan treatment (Fig. 7e). Further, rupture of microglia after ferroptosis can release pro-inflammatory damage-associated molecular patterns (DAMPs), which may enhance the innate immune system transformed from an anti-inflammatory status to a pro-inflammatory state, and thus contribute to neural damage and degeneration [79].

ScRNA-seq and bulk RNA-seq revealed temporal transcriptomic alterations reflecting crucial pathological events in SCI mice models

The 10 mm long spinal cord segments encompassing the lesion were dissected and subjected to scRNA-seq. A total of 59,558 single cells were sequenced, and 12 cell types were identified according to specific cell markers (Fig. 8a). The microglia cluster was the most abundant cell type identified in spinal cord post injury based on our scRNA-seq data. Specifically, 8 microglia subclusters were identified by dimension reduction analysis (Fig. 8b). It showed that, microglia subcluster 1 and 4 were the two main relatively homogenous microglia subpopulations identified in uninjured samples (Fig. 8c). Almost all microglial cells appeared to be transferred into the microglia subcluster 5 in four hours post injury (Fig. 8c). Microglia subcluster 5 may be transferred into subcluster 2 in one day post SCI. Infiltration degree of microglia subcluster 2 appeared to increase in one day post SCI, which were the second most abundant microglial subcluster within the lesion (Fig. 8c). Importantly, microglia subcluster 2 in one day post SCI appeared to be pro-inflammatory, because these cells decreased expression level of the M2 microglial markers and promoted expression level of M1 microglial markers (Fig. 8d). In 3-day post-SCI, the microglia subcluster 2 was still the major microglia subpopulation, while microglia subcluster 6 peaked and microglia subcluster 5 eventually disappeared (Fig. 8c). Additionally, infiltration of microglia subcluster 3 increased, and microglia subcluster 8 began to emerge (Fig. 8c). By 7 days after SCI, infiltration of the whole microglia population reached the peak, and the main subclusters consisted of clusters 2/3 and 7/8 (Fig. 8c). By 14 days post SCI, it showed a drop in the whole microglia population within the lesion site, which was consistent with the increase in microglia subcluster 2 and decreases in microglia subclusters 3, 7, and 8 (Fig. 8c). By 38 days after SCI, the whole microglia population relative to total spinal cord cells again increased, with microglia subclusters 2, 7, and 8 becoming the main microglia subclusters (Fig. 8c). Interestingly, even by 38 days after SCI, microglia population were still different from those identified within the uninjured samples, indicating persistent changes induced by microglia in the immune microenvironment post-SCI. Among the 8 identified microglia subclusters, subcluster 2 was the most heterogeneous and, which dynamically changed at different time points post injury to promote inflammatory responses.

Fig. 8

Anti-inflammatory and anti-ferroptosis effects of nanolayered double hydroxide loaded with neurotrophic factor 3 (LDH-NT3). a Unsupervised clustering of 59,558 spinal cord cells sequenced from 39 mice, color-coding annotated 12 main types of cells according to expression of specific cell markers. b Unsupervised clustering identified 8 microglia subclusters in spinal cord injury. c Line chart showing temporal alterations of the relative infiltration degrees of the 8 microglia subclusters. d Bar plot showing the temporal changes of relative infiltration degrees of M1 type microglia, M2 type microglia, and other microglia across different time points after spinal cord injury. e Box plots showing differential expression level of Nfkb1 (nuclear factor kappa B subunit 1), Hmox1 (heme oxygenase 1), and Cybb (cytochrome b-245 beta chain) across different time points after SCI (all P < 0.001) based on bulk RNA sequencing data. f mRNA analysis of M1 markers in spinal cord microglia. The data were presented as the mean ± SD *P < 0.05; **P < 0.01; ***P < 0.001. g RT-PCR analysis of Hmox1 in 3 weeks, 4 weeks, and 5 weeks post operation of SCI mice

Besides, spinal cord segments from 39 samples were also subjected to bulk RNA-seq, and the temporal expression characteristic of DEFIRGs and their TF appeared to be present, indicative of unique pathological events happening at different time points after SCI (Fig. 8e). In SCI mice, expression of Nfkb1 and Cybb increased after injury, which remained at high expression levels till 42 days after SCI (all P < 0.001). In addition, expression of Hmox1 quickly increased after SCI, reached the peak value in 1 day to 3 days post-SCI, and decreased gradually, whereas its expression remained above normal control levels thereafter (P < 0.001), which also fit well with the time scale of the rise of microglia subcluster 2. The pathophysiological staging indicated that earlier (before 3 days post SCI) anti-inflammatory and anti-ferroptosis interventions may be the key to inhibit downstream chain responses resulting in adverse clinical outcomes. Importantly, the second wave of microglia activation 14 days post SCI should also be attended, because it may play an important part in the long-term continued damages to the neurons, further obstructing neuronal regeneration. Microglial ferroptosis should be a prominent target for the next phase of therapeutic development to treat SCI.

LDH-NT3 inhibited the expression of M1 markers in spinal cord microglia and repressed ferroptosis signature in SCI mice models

The molecular effects of LDH-NT3 on spinal cord microglia were further investigated. As shown in Fig. 8f, LDH-NT3 treatment in spinal cord microglia induced a decreased expression of the M1 markers tumor necrosis factor-α (TNF-α) and inducible nitric oxide synthase (iNOS), as compared with LDH treatment, NT3 treatment, and normal control groups, based on RT-PCR results for microglia (all P < 0.05). The fold changes in Hmox1 in SCI mice in 3, 4, and 5 weeks post-operation in each group were presented in Fig. 8g, indicating that the expression of Hmox1 was significantly lower in the LDH-NT3 treatment group than in the LDH treatment group. Expression of Hmox1 in LDH-NT3 treatment group even nearly recovered to the control level in 4 weeks after operation.

Discussion

Here, we determined key biomarkers related to ferroptosis and immune microenvironment after SCI temporally and spatially, providing clues for subsequent studies at the transcriptional and sing-cell levels. WGCNA was performed to identify the hub gene module that was positively related to the time points after SCI. We identified 2 DEFIRGs (CYBB and HMOX1) and 13 DETIRGs (CD86, CXCR4, GBP2, TEC, NFATC1, NR1H3, PLXNC1, IL10RB, ITGAL, MMP12, TLR7, TLR8, and UNC93B1) based on known FRGs and IRGs respectively, then constructed PPI and microRNA-mRNA-TF networks and conducted functional enrichment analyses. The key mRNA-TF interaction pairs were all validated by ChIP-seq data. To make our findings clinically meaningful, various publicly available RNA-seq and scRNA-seq data from SCI patients or SCI mice models were used to validate our results. Based on the CMAP database, we also predicted the top 10 small-molecular drugs targeting the 2 DEFIRGs and 13 DETIRGs respectively, which may reverse the expression level of these genes.

Specifically, DEFIRGs including CYBB and HMOX1 were drivers of ferroptosis and both were up-regulated, indicating ferroptosis was naturally activated post SCI [80]. As a part of cytochrome B, CYBB may induce oxidative stress and ferroptosis after SCI [81]. Overexpression of HMOX1 can increase the ROS level and iron concentration, thereby promoting ferroptosis [15]. CYBB was the upstream molecule of hypoxia-inducible factor 1 (HIF-1) pathway, in which HIF-1α was the effector molecule and HMOX1 was the downstream product [80]. Hypoxia and ischemia of injured spinal cord may lead to HIF-1α accumulated and activated the HIF-1 pathway, which may upregulate HMOX1 and exacerbate oxidative stress via releasing ferrous iron [81]. As for DETIRGs, CD86 was one of the costimulatory receptors [82], which was a known M1 macrophage-related marker, and the M1 macrophages can secrete proinflammatory factors and aggravate neuroinflammation [83]. The chemokine receptor CXCR4 was widely expressed in immune cells and CNS cells [84], which can regulate migration of resting leukocytes in response to stromal cell-derived factor-1 (SDF1) [6,7,8,9]. Attraction of CXCR4 + macrophages was an important programmed response to SCI, which happened between 7 and 14 days after SCI [85]. GBP2 was a member of 65–67 kDa GTPases family, which was a marker for IFN responsiveness, being induced to transcription induction by IFN-γ in a few minutes [86]. A previous study showed GBP2 co-localized with ionized calcium binding adaptor molecule 1 (IBA1), and highlighted the vital role of GBP2 in microglia activation and neuronal apoptosis after CNS trauma [87]. TEC was a member of Tec family of nonreceptor tyrosine kinases, which was implicated in the modulation of cell proliferation, differentiation, survival, and cytokines expression in neutrophils following various stimulus [88]. TEC was also a critical mediator of the inflammatory responses in macrophages [89]. NFATC1 was a Ca2 + -responsive member of the nuclear factors of activated T cells (NFATs), a family of TFs playing an important role in the immune responses [90]. Nuclear receptor NR1H3 was an important modulator of macrophages and lipid homeostasis [91, 92], especially in the inflammatory responses of macrophages [93]. PLXNC1 encoded a member of the plexin family, which was related to the inflammatory responses and immune reactions [94]. A previous study showed PLXNC1 directly affected transmigration of leukocytes to enhance acute inflammatory responses [95]. IL10 exerted its effects via interleukin 10 receptor alpha (IL10RA) and IL10 receptor beta (IL10RB), which can suppress proinflammatory cytokine production, decrease macrophage activity, inhibit monocyte differentiation to dendritic cells, and downregulate costimulatory molecules on antigen-presenting cells [96]. ITGAL encoded an integrin component of lymphocyte function-associated antigen-1 (LFA-1), which provided a critical signal for the links between T cells and other cells and played a vital part in immune and inflammatory responses [97, 98]. MMP12 was known as macrophage metalloelastase and was implicated in the breakdown of extracellular matrix, which was highly expressed in chronic SCI [99]. Endosomal Toll-like receptors (TLRs) including TLR7 and TLR8 have been demonstrated to be expressed in neurons, which acted as negative regulators of neurite growth and mediators of neural apoptosis [100]. As an endoplasmic reticulum (ER)-resident transmembrane protein, UNC93B1 was widely expressed in CNS, which regulated TLR trafficking from the ER in immune cells, including dendritic cells and bone-marrow-derived macrophages [101]. Modulation of UNC93B1 was critical for regulating UNC93B1-dependent TLR function in CNS [102].

Further, we queried whether the transcriptomic changes of DEFIRGs in peripheral blood have diagnostic values for SCI. The overlapping DETIRGs identified by taking the intersection of training set (GSE45006) and validation set from blood RNA-seq dataset (GSE151371), were GBP2, TEC, UNC93B1, PLXNC1, NFATC1, IL10RB, and TLR8. These 7 genes represented chronic SCI-specific genes, which were also peripheral blood biomarkers for SCI. These findings indicated that early SCI gene expression were distinctly different chronic-stage SCI gene expression patterns, and aberrant peripheral immune responses may exist in chronic SCI patients. Therefore, these validated DETIRGs may inform possible and promising candidates to investigate the mechanisms related to SCI progression at late stage and the development of therapeutic strategies to slow or halt progression of SCI. Additionally, HMOX1 and CYBB were the key DEFIRGs in PPI and microRNA-mRNA-TF regulation networks, and NFKB1 was the common TF that can directly regulate the expression of these DEFIRGs.

The microRNA-mRNA-TF regulation networks indicated that three miRNAs (has-miR-873, has-miR-196a, and has-miR-16) may inhibit the expression of HMOX1, and six miRNAs (has-miR-186, has-miR-138, has-miR-369-3p, has-miR-551b, has-miR425, and has-miR-653) may repress the expression of CYBB. A previous study showed that, downregulation of miR-16-5p may upregulate expression of Apelin-13 to activate ERK1/2 pathway, thus alleviating neuronal apoptosis and inflammation responses after SCI [103]. Upregulation of miR-186-5p has been demonstrated to improve neurological outcomes after SCI and can repress neuroinflammation via Wnt5a-, TLR3-, or CXCL13-involved signaling pathways after SCI [104]. However, downregulation of miR-138 after SCI might have deleterious effects on neurons by activating pro-apoptotic factors, especially on spinal neurons [105]. Taken together, these identified miRNAs have been demonstrated to affect the processes of neuronal regeneration after SCI by regulating neuroinflammation or neuronal apoptosis, whereas their effects on DEFIRG expression still needed to be validated.

CYBB and HMOX1 were persistently high-expressed after SCI till the chronic stage, which were also differentially expressed in peripheral blood cells between SCI samples and HCs. Our results informed an alternative approach to consider that circulating blood cells represented “sensors” of SCI–induced molecular changes in CNS immune responses and ferroptosis. If associations can be demonstrated between the transcriptomes of the SCI lesions and that of the peripheral blood, novel biomarkers for molecular diagnostics, measurements of outcomes, and possible therapies for SCI can be established.

Results of functional enrichment analysis indicated that DEFIRGs were mainly enriched in cellular response to oxidative stress and various cytokine biosynthetic and metabolic processes, while DETIRGs were mainly associated with immune system process and inflammation responses. These findings indicated that several pro-inflammatory pathways might be the common downstreams of DEFIRGs and DETIRGs.

HMOX1 played a key role in ferroptosis pathway, functioning as a double-edged sword in ferroptosis mechanisms. Specifically, it can promote the anti-oxidant ability of cells to resist ferroptosis [106]. Additionally, upregulation of HMOX1 may promote ferrous ion accumulation to trigger ferroptosis [107]. As a part of cytochrome B, overexpression of CYBB may lead to oxidative stress after SCI [81]. NFKB1 played a prominent part in modulating inflammation, innate and adaptive immunity, cell proliferation, cell differentiation, and programmed cell death [108]. NFKB1 has been demonstrated to regulate the neuroinflammation [109] and cell death in the CNS injuries [110], which was associated with the potential mechanism of SCI [111]. NFKB played a pro-inflammatory and pro-ferroptosis role in SCI [112, 113].

To further explore and validate the functional mechanisms and cellular location of the DEFIRGs and DETIRGs at the single cell level, scRNA-seq dataset GSE162610 of 66,428 single cells was merged in the UMAP analysis, with 39 clusters and 15 cell types clearly annotated. Interestingly, DEFIRGs (CYBB, HMOX1), DEFIRG-associated TFs (NFKB1, SPI1), and various DETIRGs (CD86, CXCR4, IL10RB, TLR7, and UNC93B1) were mainly high-expressed in microglia, which further highlighted the key role of microglia in ferroptosis and immune-related pathological mechanisms after SCI. Further, a total of 19,687 microglia were categorized into 6 subclusters based on their distinct gene expression patterns. Extensive intercellular communication was identified among these microglial subclusters within the SCI lesions. Either through microglia-secreted factors or direct cell-to-cell contact among these cells, microglia functioned as active and key regulators of immune microenvironment.

Interestingly, CYBB, HMOX1, NFKB1, and various M1 type polarization markers were high-expressed and co-localized in a unique microglia subcluster (Microglia3) based on the scRNA-seq analysis of GSE162610, and this unique microglia subcluster showed early-occurrence and differentiated latter after SCI. Based on the pseudo-time analysis, there were still a considerable number of these unique microglia until the late stage of differentiation. It suggested that ferroptosis and M1 polarization were both enhanced in Microglia3 subcluster, which had potential pro-inflammatory capacities, possibly because of M1 type microglia may be better able to survive in ferroptosis-related processes. For the validation dataset GSE226286, the aged-like microglia (iron-enriched cells exhibiting ferroptosis phenotype) obtained from 3 × DR mice, showed high expression of CYBB and HMOX1. Taken together, our findings indicated that, upregulation of NFKB1 triggered by inflammatory stimulus could promote the expression levels of CYBB and HMOX1 and abundance of M1 type microglia after SCI.

Inflammasomes activation in microglia may induce metabolic shift from oxidative phosphorylation to glycolysis, which can contribute to the polarization of microglia to M1 type, eventually leading to neuroinflammation and neurodegeneration [114]. Additionally, iron enrichment can induce microglia take a pro-inflammatory and ferroptosis phenotype, and myelin debris in SCI lesions was the major source of microglial iron content [115]. On the other hand, M2 phenotype microglia was more sensitive to ferroptosis, and inhibition of M2 type microglia can attenuate neuronal inflammation.

To reveal the underlying ferroptosis programs that were changed in spatially specific manners post SCI, we comprehensively investigated the transcriptomes of different injured spinal cord segments (middle, caudal, and rostral sites) in dataset GSE69334. It indicated that CYBB, HMOX1, and NFKB1 were all significantly overexpressed in the middle sites of injured spinal cord, compared to caudal and rostral sites of injured spinal cord and uninjured spinal cord tissues. Besides, various inflammatory pathways, such as hallmark interferon alpha response, hallmark interferon gamma response, hallmark TNFA signaling via NFKB, hallmark IL2 STAT5 signaling, hallmark IL6 JAK STAT3 signaling, and hallmark inflammatory response were also significantly upregulated in the middle sites of injured spinal cord.

Currently, various biomaterials have been designed and have exhibited advantages for SCI recovery [116,117,118]. Several biomaterials, including fiber-collagen scaffolds, PLGA nanoparticles (NPs), and chitosan have been demonstrated to deliver factors/drugs and promote regeneration of neurons following transplantation in the lesions [119,120,121]. The nanocomposites are attractive for SCI therapy due to their low toxicity and promising recovery effects. Nevertheless, most of the beneficial effects are derived from loaded drugs or factors, and few of the biomaterials have immunoregulation capacities. LDH was a kind of clay showing anion exchange capacities [122]. It had various biological capacities, such as pH sensitivity, excellent biocompatibility, safe biodegradation, anti-inflammatory abilities, and high loading efficiency of genes/factors/drugs. LDH NPs have been widely investigated in biomedical engineering, including regulation of immune responses and regulation of stem cell fate [123]. Additionally, its anionic and anion exchange capacities endowed LDH with the properties to gently load vulnerable drugs/factors without weakening their effects. Our findings indicated that NT3-chitosan promoted angiogenesis and inhibited inflammatory pathways and ferroptosis markers in the middle sites of SCI lesions, creating an excellent environment to generate newly born nerves and facilitate functional recovery after SCI, which was consistent with our previous work [8]. Furthermore, LDH-NT3 inhibited the expression of M1 markers in spinal cord microglia. As a carrier, LDH loaded with NT3 showed better anti-inflammatory and anti-ferroptosis effects than LDH and NT3 itself. Therefore, LDH-NT3 can be utilized to construct a favorable microenvironment for functional recovery after SCI.

There were still several limitations in this study. Firstly, there exist differences among mice, rats, and humans in immune system activation and response to SCI. It is important to understand the potential limitations of extrapolating data from mice/rats to humans. Caution is required when extrapolating our results from mouse and rat studies to the clinic. Secondly, drawing inference from disparate tissue types can introduce discrepancies into the results. The mechanisms of secondary injury of TBI and SCI are heterogeneous, and the classification of TBI and SCI as more homogenous disease entities has contributed to the failure of pharmacotherapy trials. Although we performed spatial analysis using mice brain and spinal cord samples and transcriptomic analysis using peripheral blood samples in our validation analysis, data from post-mortem human SCI tissue samples may provide more reliable results.

Conclusions

Our findings highlighted the significance of a unique M1 type microglia subcluster exhibiting ferroptosis phenotype in injured spinal cord, which appeared early and last a certain continuous period after SCI. NFKB1 may regulate expression of ferroptosis-related markers CYBB and HMOX1 that were specifically expressed in these microglia, which may trigger neuroinflammation and induce neuronal degeneration within the middle site of SCI, and these effects may be reversed by NT3-chitosan and LDH-NT3 treatment. These hub genes could serve as both therapy targets in CNS and diagnostic markers in peripheral blood for SCI patients. These conclusions may interpret the immune- and ferroptosis-related molecular and cellular mechanisms underlying SCI and provide novel therapeutic strategies.