Unveiling the Temporal Dynamics and Molecular Regulation Profiles of Neutrophil Extracellular Traps Following Spinal Cord Injury

Introduction

Spinal cord injury (SCI) represents a devastating trauma to the central nervous system, characterized by the disruption of the spinal cord’s architecture and functionality as a result of external forces, resulting in impaired sensory, motor, and autonomic functions below the site of injury.^1–3 Globally, there are roughly 900,000 new instances of SCI annually, severely compromising patients’ quality of life and placing substantial economic burdens on both families and society.^4,5 The pathological process of SCI can be divided into two interrelated phases: primary injury, which is directly caused by the initial mechanical forces, and secondary injury cascade, which continues for days to months afterward.⁶ The secondary injury involves intricate pathophysiological mechanisms, including local inflammatory responses, oxidative stress, ferroptosis, and apoptosis, ultimately resulting in neurological degeneration and permanent neurological deficits.^7,8 These secondary injury reactions collectively lead to the expansion of the injury area and worsening of functional deficits, and are the main targets for current therapeutic interventions in SCI.^9,10

Immune-mediated inflammatory processes significantly contribute to the cascade of secondary damage following SCI. After SCI, neutrophils rapidly migrate from circulation into the injured spinal tissue, representing the earliest immune cell population to infiltrate the damaged area.¹¹ As critical effector cells of innate immunity, neutrophils participate in inflammatory responses primarily through phagocytosis, release of antimicrobial proteins, and formation of neutrophil extracellular traps (NETs). NETs are network-like structures released by neutrophils, mainly composed of DNA, histones, and various antimicrobial proteins.¹² Recent studies have revealed that NETs can exacerbate inflammatory responses through multiple mechanisms in process of SCI.¹³ The DNA and histone components released from NETs function as damage-associated molecular patterns that trigger microglial and astrocytic activation, subsequently enhancing pro-inflammatory cytokine production.¹⁴ Moreover, the proteolytic enzymes contained within NETs can directly degrade extracellular matrix structures and disrupt intercellular junction proteins, thereby compromising blood-spinal cord barrier function and exacerbating tissue edema and secondary tissue damage.¹⁵ Recent investigations have demonstrated the therapeutic potential of targeting NETs in SCI, with recombinant human DNase treatment shown to mitigate extracellular trap-mediated damage and improve long-term recovery in experimental models, while CD47-blocking antibodies have been found to interfere with NET formation and reduce spinal cord edema.^16,17 These NETs-mediated secondary injuries can amplify inflammatory responses and aggravate neuronal damage, but effective inhibition of NETs can significantly promote neural tissue repair and improve function.¹⁸ Therefore, NETs are not merely a consequence of inflammatory responses but serve as key drivers in the amplification and perpetuation of inflammation.¹⁹ However, the comprehensive molecular mechanisms of NETs in SCI have not yet been fully elucidated, particularly regarding the key genes associated with NETs formation and their regulatory networks, which remain to be systematically investigated.

In this study, we initially identify key genes linked to the formation and function of NETs in SCI through integrated differential expression analysis and weighted gene co-expression network analysis (WGCNA) methodology. Secondly, we used single-cell transcriptome analysis to explore the spatiotemporal expression patterns and interactions of these genes in neutrophils and other cell subsets, and determined the critical period of NETs formation within 1 to 7 days after SCI using a mouse SCI model. Subsequently, we also explored the NETs-related transcription factor (TF)-mRNA and miRNA-mRNA regulatory networks in SCI, and predicted potential therapeutic agents targeting NETs. These investigations will provide novel theoretical foundations for understanding the role of NETs in SCI and scientific evidence for the future development of precision therapeutic strategies targeting NETs.

Materials and Methods

Data Acquisition

For this study, we obtained five publicly available datasets from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). Two mouse SCI mRNA expression datasets (GSE42828 and GSE5296) were used for differential expression analysis. GSE42828 comprised 17 samples (4 Sham controls, 5 SCI-1d, 4 SCI-3d, and 4 SCI-7d), while GSE5296 included 17 samples (8 Sham controls, 3 SCI-1d, 3 SCI-3d, and 3 SCI-7d). The sex and age information for mice used in GSE42828 and GSE5296 datasets were not mentioned in the original database documentation. Moreover, we downloaded mouse SCI single-cell transcriptome data from GSE162610, which included 10 samples (3 Sham controls, 3 SCI-1d, 2 SCI-3d, and 2 SCI-7d) derived from 8–10-week-old female C57BL/6 mice. We also obtained human peripheral blood RNA expression profiles from SCI patients (GSE151371). This dataset comprised 48 samples (10 healthy controls and 38 SCI patients). Additionally, in this study, a set of 136 NETs-related genes (NRGs) was compiled from a series of previous studies (Table S1).^20–23

Data Processing and the Identification of NETs-Related DEGs

The expression data from GSE42828 and GSE5296 were downloaded from GEO database. The two datasets were then merged and the batch effect was corrected using the Combat function in the sva package in R. (v4.2.1). Principal component analysis (PCA) was conducted utilizing the FactoMineR package to visualize sample clustering after batch correction and confirm the effectiveness of batch effect removal. The limma package was utilized to conduct differential expression analysis, comparing the SCI-1d, SCI-3d, and SCI-7d groups against the Sham group. Differentially expressed genes (DEGs) were determined according to the following criteria: |log2(Fold Change)| > 1 and adjusted P-value < 0.05. Volcano plots were generated using the ggplot2 package to visualize the distribution of DEGs.

To identify differentially expressed NETs-related genes (DENRGs) in SCI, we intersected the DEGs from each time point with 136 NRGs collected from previous studies. Venn diagrams were created using the VennDiagram package to illustrate the overlap between DEGs at different time points and NRGs. The expression patterns of the 22 overlapping DENRGs (shared across all time points) were visualized using heatmaps generated by the pheatmap package.

Additionally, NETs scores for each sample were determined employing the gene set variation analysis (GSVA) algorithm with the 136 NRGs serving as the reference gene set to estimate the relative NETs-related gene expression patterns across different time points post-SCI.

Functional Enrichment Analysis

The Gene Set Enrichment Analysis (GSEA) was utilized to evaluate the distribution of genes within predefined gene sets in a gene list ranked by phenotypic correlation, aiming to uncover the potential impact of these gene sets on the phenotype. We obtained the gene set “c5.all.v7.2.symbols.gmt” from the Molecular Signatures Database (MSigDB). For our analysis, we specified the number of sample permutations as 1000 and established a significance threshold with a P-value of less than 0.05. To delve into the biological mechanisms and relevant signaling pathways involved in NETs formation in SCI, we converted gene symbols to Entrez gene IDs with “org.Mm.eg.db”. Subsequently, we performed Gene Ontology (GO) functional annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis on shared DENRGs using the ClusterProfiler package. Finally, we presented the enrichment analysis results visually by utilizing the ggplot2 and GOplot2 packages.

Protein-Protein Interaction (PPI) Network Construction and Hub Gene Screening

The 22 shared DENRGs were uploaded to the STRING database (http://string-db.org), with the threshold for the minimum interaction score was set at 0.4, disconnected nodes hidden, and “Mus musculus” selected as the organism. We utilized Cytoscape software (version 3.7.2) to visualize the PPI network and perform subsequent analysis. CytoHubba was used to screen the top 10 hub genes by degree value. Subsequently, key modules were discerned utilizing the Molecular Complex Detection (MCODE) plugin in Cytoscape, with the parameters set as a node score cut-off of 0.2, a degree cut-off of 2, a k-score of 2, and a maximum depth of 100. Spearman correlation analysis was performed on the 10 Hub genes.

WGCNA

To identify gene modules associated with SCI, we employed WGCNA on the gene expression data. Initially, we filtered the gene expression dataset to retain only genes exhibiting significant expression variability across samples. Subsequently, we conducted hierarchical clustering of the samples to identify and exclude outliers, thereby ensuring high data quality. For the construction of the network, we assessed the fit of the scale-free topology model and determined the optimal soft thresholding power to be β = 5. Subsequently, we computed the adjacency matrix and converted it into a Topological Overlap Matrix (TOM) to bolster the network’s robustness. Dynamic tree cutting was applied to cluster genes into modules, and highly similar modules were merged. To assess the relationship between modules and phenotypes, we calculated module eigengenes and performed Pearson correlation analysis with phenotypic data, comparing the Sham and SCI groups. Highly correlated genes were identified by setting thresholds for gene significance (GS > 0.5) and module membership (MM > 0.8). The results of modular trait association heatmaps and scatter plots of GS versus MM were visualized, which provided intuitive insights into the relationship between gene expression patterns and phenotypes in the context of SCI. Finally, the key DENRGs in SCI were determined by the intersection of the highly correlated genes in the blue module, the high-degree value Hub genes screened by PPI, and the NETs-related gene set.

Single-Cell Sequencing Data Analysis

The Seurat R package was used to perform downstream analysis on the public single-cell sequencing dataset GSE162610 following the methodology described in the original publication.²⁴ Single-cell RNA sequencing data were aligned to the mouse reference genome mm10 to maintain consistency with the original study and ensure reproducibility of results. To ensure the robustness of analyses, we retained genes expressed in at least 10 cells per sample. Normalization and scaling were performed using the SCTransform algorithm, with cell cycle effects controlled by regressing out G2M and S phase scores.^25,26 To mitigate potential batch effects, batch correction was applied prior to identifying 3000 highly variable genes for PCA. We chose the top 20 principal components to create a shared nearest neighbor (SNN) graph, which facilitated the unsupervised clustering of individual cells. Our analytical framework builds upon the established cellular annotations and gene expression profiles from the original study. We employed SingleR, leveraging established reference datasets to enhance classification accuracy.^27–29

To explore dynamic changes in cell type distribution and gene expression across time points, data processing was performed using the dplyr package, while ggplot2 was used for visualization. Bubble plots were generated to illustrate average gene expression levels and cell-type proportions. Additionally, neutrophil gene expression patterns were examined using the FeaturePlot function in Seurat, providing insights into their functional states and cellular heterogeneity.

Experimental Animals

A total of 70 female C57BL/6 mice (8 weeks old, 18–20g, Vital River Laboratory Animals Co., Beijing) were utilized in this study. All experimental animals were maintained in pathogen-free environments with regulated atmospheric conditions (ambient temperature, moisture levels), cyclical photoperiod alternations of equal duration, and unrestricted access to nutrition and hydration. Mice were group-housed (4–5 mice per cage) in standard laboratory cages throughout the experimental period. The entirety of animal research protocols adhered to established laboratory animal welfare guidelines and regulatory standards. All experimental procedures were conducted in the animal unit at Beijing MedCona Animal Experiment Center in Beijing, China, in accordance with protocols approved by the institutional animal welfare and ethical committee (Permit No. MDKN-2024-077).

Construction of Spinal Cord Injury Model

For mice in the SCI group, after anesthetizing with isoflurane (initial dose 4%, maintained at 2%) (RWD, R510-22, Guangdong, China), All mice underwent a midline incision followed by blunt dissection of the tissue overlying the T7-T9 vertebral laminae, after which a T8 vertebral laminectomy was performed. T8 spinal cord contusion model was prepared by the NYU Impactor-III (W.M.Keck, USA). Spinal cord contusion injury was performed by 5-g node falling freely from the height of 12.5 mm. After suturing the incision, the mice were manually assisted with bladder emptying twice daily and administered a cefuroxime injection once daily for three days post-surgery. For mice in the Sham group, only a T8 vertebral laminectomy was carried out, with the spinal cord remaining unharmed.

Isolation and Extraction of Primary Neutrophils

The mice were anesthetized through inhaling isoflurane, and subsequently, the femur and tibia were dissected. The tibial and femoral cavities were irrigated with PBS utilizing a 25G needle, and the resulting effluent was passed through a 70 µm filter mesh to generate a single-cell bone marrow suspension. Negative selection methodology employing the EasySep Mouse Neutrophil Enrichment Kit was utilized for the isolation and purification of neutrophil populations (19762; STEMCELL) according to the manufacturer’s instructions. The purified neutrophil was maintained in RPMI 1640 culture medium (Sigma Aldrich). For NET induction, neutrophils were incubated with the classical NET inducer Phorbol 12-Myristate 13-Acetate (PMA, 50 nM, Sigma-Aldrich) for 4 hours in a 37°C, 5% CO₂ incubator.

Immunofluorescence Staining

Following isoflurane anesthesia, experimental animals underwent cardiac perfusion with PBS followed by 4% paraformaldehyde (PFA). Spinal cord segments containing the injury epicenter were extracted and processed through a sequential protocol: fixation in 4% PFA for 24 hours, dehydration in 30% sucrose solution for 72 hours, and embedding in OCT compound. Cryosections were prepared at 10μm thickness. After blocking with 5% donkey serum, the tissue sections were incubated with primary antibodies at 4°C overnight, and subsequently incubated with secondary antibodies for a duration of 1 hour at 37°C. Thorough washing with PBS was performed before mounting with anti-fade medium containing DAPI. Specimens were examined using confocal microscopy (Leica TCS SP8). Primary antibodies utilized included rabbit anti-Histone H3 (citrulline R2+R8+R17) (1:100; Abcam, ab5103) and Goat anti-Myeloperoxidase / MPO (1:100; R&D Systems, AF3667). For quantitative analysis of NETs formation, MPO+CitH3+ double-positive cells were counted in the injury epicenter region. Observations were conducted on mid-sagittal sections of spinal cord specimens from each animal. Three representative images at 40× magnification were randomly selected from the injury epicenter on the mid-sagittal sections. Then, MPO+CitH3+ cells were manually counted by two independent investigators blinded to the experimental groups within these three standardized fields. The average count from the three images was calculated for each animal, and results were expressed as the number of MPO+CitH3+ cells per field. Cells were considered double-positive only when clear co-localization of both MPO (green) and CitH3 (red) signals was observed with nuclear or perinuclear staining patterns characteristic of NETs formation.

Western Blot

Total protein was extracted by thoroughly lysing 15mg of tissue obtained from the SCI epicenter using RIPA buffer supplemented with protease and phosphatase inhibitors. The concentration of the proteins was determined utilizing the BCA assay (Thermo Scientific). The protein specimens following denaturation underwent electrophoretic separation via SDS-PAGE (Epizyme) and were subsequently electroblotted onto PVDF membranes. These membranes were then subjected to blocking treatment using a 5% non-fat milk solution at room temperature, followed by overnight incubation with primary antibodies at 4°C. The following day, after extensive washing with TBST, membranes were incubated with HRP-conjugated secondary antibodies for 1 hour at room temperature. Protein bands were visualized using ECL reagent (Thermo Scientific) and detected using a chemiluminescence imaging system (SCG-W3000, Servicebio). For protein quantification, target band intensities were measured using ImageJ software (version 1.52a, NIH, USA). Briefly, Western blot images were converted to 8-bit grayscale, and rectangular areas of equal size were drawn around target band. After background subtraction, the integrated density values were recorded. All target protein signals were normalized to corresponding GAPDH internal control bands from the same membrane. Primary antibodies utilized in this study included rabbit anti-Histone H3 (citrulline R2+R8+R17) (1:1000; Abcam, ab5103) and rabbit anti-MPO (1:2000; Abcam, ab208670).

Quantitative Real-Time PCR (qPCR)

The expression of key DENRGs was validated using qPCR. Specifically, RNA was extracted from neutrophils or from homogenized spinal cord tissue containing the injury epicenter using cold Trizol reagent (Solarbio, China), followed by chloroform addition and centrifugation. The aqueous phase was collected, added isopropanol and centrifugated again. After washing by absolute ethyl alcohol, RNA was dissolved in DEPC water. The reverse transcription of RNA was carried out using the RT reagent Kit (Takara, China). The qPCR of cDNA was conducted by TB Green Premix Ex Taq II (Takara, China) according to the manufacturer’s instruction. Expression levels of each gene were normalized by GAPDH mRNA expression, and managed by 2^−ΔΔCt method to calculate the relative expression. The primers we utilized are shown in Table S2.

Construction of TF-mRNA and miRNA-mRNA Regulatory Networks

NetworkAnalyst (https://www.networkanalyst.ca/faces/home.xhtml) functions as a publicly available web-based visualization platform that facilitates thorough gene expression profiling and integrative meta-analysis. To construct the regulatory network of NETs-related TFs-target genes, we submitted 22 screened shared DENRGs to the platform. NETs-related TFs were predicted and generated from ENCODE database, and TF-DENRGS relation pairs were retained for further visualization and analysis based on Cytoscape software.

To enhance accuracy and reliability, three databases, including DIANA-micro T (https://diana.e-ce.uth.gr/home), miRWalk (http://mirwalk.umm.uni-heidelberg.de) and miRDB (https://mirdb.org) were used to jointly predict target miRNA upstream of DENRGs. The intersection of the predicted results from each database was considered to be the target miRNA. Ultimately, the NETs-related mRNA-miRNA regulatory network in SCI was presented utilizing Cytoscape software.

Drug Prediction and Molecular Docking Analysis

The Enrichr platform (https://maayanlab.cloud) is an accessible web-based database that aggregates genome-wide functionally-enriched information from a large number of gene sets. Its sub-website Drug Signatures database (DSigDB) supports the prediction of target-related potential small-molecule compounds acting on feature genes in specific disease profiles through drug-gene interaction relationships. Potential candidate drug compounds targeting NETs to treat SCI were predicted based on DSigDB via the Enrichr platform. Among the top ten, excluding 2 industrial compounds, the other eight therapeutic small molecule compounds are displayed in a heat map.

Additionally, the three-dimensional molecular representation of atorvastatin, which emerged as the top candidate, was extracted from the PubChem repository and archived in SDF format. Protein structure for the seven key DENRGs were procured from the PDB database (https://www.rcsb.org/). These protein structures are all saved in PDB format. For molecular docking analyses, we employed the web-based CB-Dock2 platform (http://clab.labshare.co.uk/cb-dock/php/).^30,31 As an innovative blind docking tool, CB-Dock2 effectively identifies potential binding regions on target proteins by implementing a curvature-dependent cavity detection algorithm to determine optimal center positions and dimensional parameters. The platform incorporates the docking engine AutoDock Vina (1.2.0), substantially enhancing docking accuracy. The Vina score below −5.0 kcal/mol indicates favorable interaction potential between protein targets and candidate compounds, while binding energy values less than −7.0 kcal/mol signify exceptionally strong molecular associations characterized by heightened binding stability.

Statistical Analysis

Statistical analyses were performed using GraphPad Prism 8.0, and the results are presented as mean ± standard deviation (SD). Sample sizes for each experiment were as follows: immunofluorescence analysis (n = 6 per group), Western blot analysis (n = 3 per group), qPCR analysis (n = 6 per group), and in vitro neutrophil experiments (n = 6 paired samples). Initial assessment included normality and lognormality examination for all datasets. The normality of data distribution was rigorously evaluated using Shapiro–Wilk test with a significance level of α = 0.05. In all study groups except for the in vitro neutrophil experiment validating Itgb2 expression, the data conformed to normal distribution assumptions. For datasets meeting normality criteria, parametric tests were employed for subsequent analyses, including paired t-test or one-way ANOVA followed by Tukey’s post hoc analysis. For datasets that do not meet the normality criteria, the nonparametric test was used. P < 0.05 were deemed indicative of statistical significance throughout the study.

Results

Identification of Differentially Expressed NETs-Related Genes in SCI

To study the role of NETs in SCI, we analyzed two public microarray datasets (GSE42828 and GSE5296) from the GEO database. After merging the datasets and correcting for batch effects, PCA showed that the two datasets were effectively integrated (Figure 1A and B). We calculated the NETs score for each sample in the integrated dataset using GSVA with the 136 NETs-related genes as the background gene set. The NETs score represents the relative enrichment of NETs-related gene expression patterns, providing a quantitative measure of NETs formation activity across different samples and time points. The result shows that the NETs scores were significantly elevated in all SCI groups compared to the Sham group, with the highest scores observed in the SCI-1d group, which to some extent revealed the temporal pattern of NETs formation after SCI (Figure 1C).

Figure 1 Identification of differentially expressed NETs-related genes in SCI. (A) PCA plot showing the distribution of samples before batch effect correction from two independent datasets (GSE42828 and GSE5296). (B) PCA plot showing the integrated data after batch effect correction, demonstrating the two datasets were effectively integrated. (C) The NETs score of each sample in different groups after SCI was calculated by the GSVA algorithm. NETs scores were significantly elevated in all SCI groups compared to the Sham group, with the highest scores observed in SCI-1d, followed by a progressive decrease at 3 days post-injury. (D-F) Volcano plots showing DEGs between SCI and Sham groups at 1 day (D), 3 days (E), and 7 days (F) post-injury. Purple dots represent upregulated genes, and green dots represent downregulated genes (|log2(Fold Change)| > 1 and adjusted P-value < 0.05). (G) Venn diagram depicting the overlap of DEGs identified at different time points (SCI-1d, SCI-3d, and SCI-7d) with the set of 136 NETs-related genes (NRGs). A total of 22 NETs-related genes were identified as common to the three groups. (H) Heatmap visualization of the expression patterns of 22 shared DENRGs across all groups. These genes were consistently upregulated at all three time points compared to the sham group. The color scale represents the normalized expression values, with red indicating upregulation and green indicating downregulation. *P < 0.05, ***P < 0.001. Abbreviations: SCI, spinal cord injury; PCA, principal component analysis; NETs, neutrophil extracellular trap; GSVA, gene set variation analysis; NRGs, NETs-related genes; DENRGs, differentially expressed NETs-related genes.

Differential expression analysis was performed between SCI groups (1d, 3d, and 7d) and the Sham group. A total of 846 DEGs were identified in the SCI-1d group, 1428 DEGs were identified in the SCI-3d group, and 924 DEGs were identified in the SCI-7d group (Figure 1D–F). To specifically identify NETs-related differentially expressed genes, we intersected the DEGs from each time point with 136 NRGs. The Venn diagram analysis revealed 22 shared DENRGs that were consistently differentially expressed across all three time points (1d, 3d, and 7d) post-SCI (Figure 1G). Notably, these 22 genes were predominantly upregulated compared to the Sham group, as visualized in the heatmap (Figure 1H). These findings demonstrate a significant and dynamic regulation of NRGs following SCI, suggesting a critical role for NETs in the acute phase of SCI.

Functional Enrichment Analysis of NETs-Related Genes in SCI

We also performed GSEA on the whole gene expression profiles at different time points after SCI. The analysis revealed significant enrichment of neutrophil-related pathways, including neutrophil migration, chemotaxis, and immune activation at all time points (Figure 2A–C). These findings suggest persistent neutrophil activity throughout the acute and subacute phases following SCI.

Figure 2 Functional enrichment analysis of NETs-related genes in SCI. (A–C) GSEA of gene expression data of SCI at different time points: 1 day (A), 3 days (B), and 7 days (C) after injury. The plots show significant enrichment of neutrophil migration, chemotaxis, and immune activation pathways at all time points. (D) Gene Ontology (GO) analysis of 22 shared DENRGs. The bar plot displays significantly enriched biological processes (BP, dark pink), cellular components (CC, medium pink), and molecular functions (MF, light pink). (E) Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of the 22 shared DENRGs. The bubble plot shows significantly enriched pathways, with bubble size corresponding to gene count and color intensity representing -log10(P value). Abbreviations: SCI, spinal cord injury; NETs, neutrophil extracellular traps; DENRGs, differentially expressed NETs-related genes; GSEA, gene set enrichment analysis; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

The purpose of the GSVA algorithm we introduced is to explore the temporal dynamics of NETs formation during neutrophil activation. The NETs score represents a composite quantitative measure derived from the coordinated expression levels of all 136 NETs-related genes, providing an assessment of the relative enrichment of NETs-related gene expression patterns within each sample. This scoring system enables comparative analysis of NETs formation activity between different experimental groups, with higher scores indicating greater NETs-related gene expression activity. As shown in Figure 1C, the NETs scores were significantly elevated in all SCI groups compared to the Sham group, with the highest scores observed in the SCI-1d group, confirming the temporal dynamics of NETs formation after SCI and supporting our subsequent functional enrichment analyses.

Additionally, GO analysis of the 22 shared DENRGs demonstrated enrichment in biological processes related to neutrophil function, including neutrophil chemotaxis, inflammatory response regulation, granulocyte migration, and leukocyte activation (Figure 2D). Enriched cellular components included membrane rafts or microdomain and receptor complex, while the enriched term of molecular function encompassed complement binding, integrin binding, and cytokine activity. KEGG pathway analysis showed significant enrichment in immune-related pathways, notably neutrophil extracellular trap formation, complement and coagulation cascades, and Toll-like receptor signaling (Figure 2E). The direct enrichment of the “neutrophil extracellular trap formation” pathway validates our study and confirms the involvement of these genes in NETs formation after SCI. These results collectively indicate that NRGs are actively involved in inflammatory processes that may influence secondary injury pathophysiology following SCI.

Establishment and Analysis of PPI Network

To investigate interactions among the 22 shared DENRGs, we constructed a PPI network based on STRING database (Figure 3A). The network comprised 21 nodes and 105 edges, indicating extensive interactions among these NRGs. MCODE algorithm revealed a highly interconnected functional module containing 14 nodes with a score of 11.07 (Figure 3B). CytoHubba analysis identified the top 10 hub genes: Itgam, Tlr4, Tlr2, Ccl2, Fcgr2b, Casp1, Itgb2, Cd44, Timp1, and Ccl3, with degree values ranging from 18 to 9 (Figure 3C and D). Correlation analysis demonstrated significant positive correlations among these hub genes (Figure 3E). Functional enrichment analysis showed significant involvement in inflammatory response, leukocyte activation, neutrophil migration, and regulation of TNF production (Figure 3F). Eight hub genes included in the gene expression profile of peripheral blood samples from SCI patients were extracted for validation. The results showed that all eight genes in the SCI group were upregulated compared with healthy controls. ITGB2, TLR2, TLR4, CASP1, ITGAM, and TIMP1 showed highly significant differences (p<0.05) (Figure 3G and H).

Figure 3 PPI network analysis and validation of hub NETs-related genes in SCI. (A) PPI network of the 22 shared DENRGs. The network consists of 21 nodes and 105 edges, demonstrating extensive interactions among these genes. (B) The most significant functional module identified by MCODE algorithm (score = 11.07) containing 14 highly interconnected nodes. (C) The top 10 hub genes extracted from the PPI network using CytoHubba plugin based on degree value. (D) Degree distribution plot of the top 10 hub genes (Itgam, Tlr4, Tlr2, Ccl2, Fcgr2b, Casp1, Itgb2, Cd44, Timp1, and Ccl3), with degree values ranging from 18 to 9. The color gradient represents the connectivity strength of each gene within the network. (E) Correlation matrix showing the significant positive correlations among the top 10 hub genes. Color intensity and asterisks indicate the strength and significance level of correlations. (F) Functional enrichment analysis of the identified hub genes, showing significant involvement in inflammatory response, leukocyte activation, neutrophil migration, and regulation of TNF production. (G and H) Validation of eight hub genes in peripheral blood samples from SCI patients compared to healthy controls. Box plots showing the expression of CASP1, CD44, FCGR2B, ITGAM (G) and ITGB2, TLR2, TLR4, TIMP1 (H) genes in SCI patients (purple) versus healthy controls (green). All eight genes were upregulated in the SCI group, with ITGB2, TLR2, TLR4, CASP1, ITGAM, and TIMP1 showing statistically significant differences. **P < 0.01, ***P < 0.001. Abbreviations: SCI, spinal cord injury; NETs, neutrophil extracellular traps; PPI, protein-protein interaction; MCODE, Molecular Complex Detection; DENRGs, differentially expressed NETs-related genes.

WGCNA Reveals SCI-Associated Gene Modules and Key NETs Genes

Given that GSVA analysis showed significantly elevated NETs scores across all SCI time points compared to the sham group, we performed WGCNA using SCI as a unified trait to identify gene modules universally associated with NETs formation following SCI. Figure 4A illustrates the determination of the optimal soft thresholding power (β = 5), selected by evaluating both the scale-free topology model fit index (left panel) and mean connectivity (right panel). Hierarchical clustering based on topological overlap generated a gene dendrogram with five distinct modules represented by different colors, including blue, yellow, red, turquoise and grey (Figure 4B). The module-trait relationship analysis revealed that the blue module exhibited the strongest positive correlation with SCI (Figure 4C). Further analysis of gene significance versus module membership in the blue module demonstrated a strong linear relationship between these parameters (cor = 0.9, p < 1e-200), with genes in the upper right quadrant (GS > 0.5, MM > 0.8) representing 412 that were highly associated with SCI (Figure 4D). Subsequently, the intersection of these 412 genes with the 136 NRGs revealed 18 overlapping genes (Figure 4E). Further intersection with the 10 PPI-derived hub genes calculated seven genes (Casp1, Ccl3, Fcgr2b, Itgam, Itgb2, Tlr2, and Tlr4) that represented key regulators in SCI-associated NETs formation (Figure 4F). Expression analysis confirmed that all seven key genes were significantly upregulated across all SCI time points compared to the sham group (Figure 4G). All identified key genes maintained elevated expression throughout the observation period, indicating their sustained involvement in the inflammatory response following SCI.

Figure 4 WGCNA reveals key co-expression modules and key NRGs associated with SCI. (A) Determination of the optimal soft thresholding power for network construction. Left panel: Scale-free topology model fit index (y-axis) as a function of soft-thresholding power (x-axis), with power 5 selected as optimal. Right panel: Mean connectivity (y-axis) decreases as the soft-thresholding power increases. (B) Gene dendrogram and module colors obtained by hierarchical clustering based on topological overlap, with branches representing genes and colors below denoting assigned modules. Five distinct co-expression modules were identified. (C) Module-trait relationships heatmap showing correlations between module eigengenes and phenotypic traits (Sham and SCI). The blue module exhibited the strongest positive correlation with SCI. (D) Scatter plot of gene significance for SCI trait versus module membership in the blue module. Genes above the horizontal red line (GS > 0.5) with high module membership (MM > 0.8) were considered highly associated with SCI, identifying 412 genes strongly correlated with the SCI phenotype. (E) Venn diagram showing the intersection between the 412 WGCNA-derived highly related genes for SCI and 136 NETs-related genes, revealing 18 overlapping genes. (F) Venn diagram illustrating the intersection of 18 WGCNA-derived NETs-related genes with the 10 hub genes identified from the PPI network, resulting in 7 key genes (Casp1, Ccl3, Fcgr2b, Itgam, Itgb2, Tlr2, and Tlr4) representing key regulators in SCI-associated NETs formation. (G) Box plots showing the expression patterns of the seven key NETs-related genes across experimental groups (Sham, SCI-1d, SCI-3d, and SCI-7d). All seven genes exhibited significant upregulation across all SCI time points compared to the sham group. ***P < 0.001. Abbreviations: SCI, spinal cord injury; WGCNA, weighted gene co-expression network analysis; NETs, neutrophil extracellular traps; PPI, protein-protein interaction; GS, gene significance; MM, module membership.

Single-Cell RNA Sequencing Reveals Cell-Specific Expression of Key NETs Genes in SCI

To explore the spatiotemporal expression characteristics of these seven key NRGs in SCI, we analyzed the single-cell RNA sequencing data of GSE162610. The UMAP visualization revealed 15 distinct cell populations within the spinal cord tissue following injury, with each cluster representing a different cell type (Figure 5A). Temporal analysis of cell distribution across different time points (uninjured, 1dpi, 3dpi, and 7dpi) demonstrated dynamic shifts in cellular composition, with notable alterations in immune cell populations, particularly neutrophils, during the acute phase post-injury (Figure 5B). The dot plot of canonical cell-type markers confirmed the accuracy of our cell type annotations, with each population exhibiting characteristic gene expression signatures (Figure 5C). Quantitative assessment of absolute cell numbers and relative cell proportions across different time points revealed a marked infiltration of neutrophils at 1dpi that progressively diminished at later timepoints (Figure 5D). Expression profiling of the seven key NETs-related genes (Tlr4, Tlr2, Itgb2, Itgam, Fcgr2b, Ccl3, and Casp1) across all identified cell types demonstrated their predominant expression in myeloid populations, with particularly robust expression in neutrophils, monocytes, and macrophages (Figure 5E). UMAP visualization of specific genes revealed distinctive expression patterns, with Ccl3 showing particularly pronounced expression in neutrophils, whereas Casp1 exhibited a more dispersed pattern. Spatial expression mapping of Fcgr2b, Itgam, Itgb2, Tlr2, and Tlr4 consistently demonstrated their enrichment in neutrophils and other myeloid populations, providing further evidence of their importance in neutrophil function and NETs formation following SCI (Figure 5F and G).

Figure 5 Single-cell RNA sequencing analysis reveals cell-specific expression patterns of key NETs-related genes in SCI. (A) UMAP visualization of 15 distinct cell populations identified from the single-cell RNA sequencing dataset (GSE162610), with each color representing a different cell type and cells clustered based on transcriptional similarity. (B) UMAP plots showing the distribution of cell populations across different time points post-SCI (uninjured, 1dpi, 3dpi, and 7dpi), revealing temporal changes in cellular composition after injury. (C) Dot plot showing the expression patterns of canonical cell type-specific marker genes across the 15 identified cell populations, confirming accurate cell type annotation. (D) Bar charts illustrating changes in absolute cell numbers (left) and relative cell proportions (right) across different time points post-SCI. Neutrophils exhibit marked infiltration at 1dpi that progressively decreases at later timepoints. (E) Dot plot showing the expression levels of the seven key NETs-related genes (Tlr4, Tlr2, Itgb2, Itgam, Fcgr2b, Ccl3, and Casp1) across all identified cell types. Dot size represents the percentage of cells expressing each gene within a cell type, while color intensity indicates the mean expression level. The analysis demonstrates predominant expression of these genes in neutrophils, monocytes, and macrophages. (F) UMAP feature plots showing the expression distribution of Casp1 (left) and Ccl3 (right) across all cells, with color intensity indicating expression level. Ccl3 shows particularly robust expression in neutrophils, while Casp1 exhibits a more dispersed expression pattern. (G) UMAP feature plots showing the expression patterns of Fcgr2b, Itgam, Itgb2, Tlr2, and Tlr4 across all cells. All five genes demonstrate strong expression in neutrophils and other myeloid populations, with varying intensity across different cell types. Abbreviations: SCI, spinal cord injury; NETs, neutrophil extracellular traps; UMAP, Uniform Manifold Approximation and Projection; dpi, days post-injury; OPC, oligodendrocyte progenitor.

Further temporal dynamic analysis focusing specifically on neutrophil populations revealed that all seven key genes showed highest expression at 1dpi, with the most concentrated neutrophil clustering observed at this timepoint (Figure 6A). Single-cell level quantitative analysis of neutrophil populations revealed distinct temporal expression patterns for the seven key NETs-related genes (Figure 6B). The analysis demonstrated that 1dpi showed the highest neutrophil infiltration (greatest number of expressing cells) combined with the most robust overall gene expression levels for most NRGs. Expression levels and neutrophil presence progressively decreased at 3dpi and 7dpi. These results confirmed that NETs-related gene activity basically peaked at 1dpi after SCI, coinciding with the peak of neutrophil infiltration. This temporal pattern also aligns with the highest NETs score observed in our transcriptomic analysis, potentially suggesting that NETs formation peaks during the acute phase of SCI.

Figure 6 Temporal expression dynamics of key NETs-related genes in neutrophil populations following SCI. (A) UMAP visualization showing the temporal expression patterns of the seven key NETs-related genes specifically in neutrophil populations across different time points post-SCI (uninjured, 1dpi, 3dpi, and 7dpi). Each row displays the expression pattern of a specific gene, with color intensity representing expression level (red indicating high expression, blue indicating low expression). (B) Single-cell level quantitative analysis of expression levels of the seven key NETs-related genes in neutrophils across different time points. Each dot represents an individual neutrophil expressing the corresponding gene. The analysis results reveal both the number of neutrophils present in each time group and the overall expression level of the corresponding gene expressed by neutrophils. The quantitative analysis demonstrates that 1dpi shows the highest neutrophil infiltration (as evidenced by the greatest number of data points) and the most robust overall gene expression levels for the seven key NETs-related genes. The progressive reduction in both neutrophil numbers and expression levels at 3dpi and 7dpi supports the temporal dynamics of NETs formation peaking in the acute phase post-SCI.

NETs Formation Peaks at Day 1 After SCI

To validate our transcriptomic findings in vivo, we examined established markers of NETs by establishing a mouse model of spinal cord contusion injury at T8 level. The schematic diagram of the whole experiment is shown in Figure 7A. Immunofluorescence analysis demonstrated prominent co-localization of MPO and CitH3 at the lesion epicenter in the SCI-1d group, with diminishing signals observed at 3 and 7 days post-injury (Figure 7B). Quantitative assessment revealed significantly elevated MPO+CitH3+ cell counts in SCI-1d specimens compared to SCI-3d, SCI-7d, and sham controls [F(2, 15) = 28,78, P = 7.33e-06] (Figure 7G).

Figure 7 NETs formation peaks at day 1 post-SCI with progressive attenuation at later time points. (A) Schematic diagram of the entire experimental design process. (B) Representative immunofluorescence images showing the co-localization of MPO (green) and CitH3 (red) at the lesion epicenter across different time points after SCI (Sham, SCI-1d, SCI-3d, and SCI-7d). The upper panel shows low-magnification images of the injury site (scale bars = 300 μm), while the lower panel displays high-magnification images of the boxed areas (scale bars = 50 μm). DAPI (blue) stains cell nuclei. (C and D) Representative Western blot images showing the protein expression of MPO (C) and CitH3 (D) in spinal cord tissue at different time points after SCI, with GAPDH serving as control. (E and F) Quantitative analysis of MPO (E) and CitH3 (F) protein expression levels normalized to GAPDH. Both markers exhibit peak expression at 1 day post-SCI, followed by significant decreases at days 3 and 7 (n = 3). (G) Quantification of MPO+CitH3+ cells per field across different experimental groups, showing significantly elevated NETs formation at day 1 post-SCI compared to Sham, SCI-3d, and SCI-7d groups (n = 6). Data are presented as mean ± SD. Statistical analysis was performed using one-way ANOVA followed by Tukey’s post hoc test. *P < 0.05, ***P < 0.001, ****P < 0.0001.

Western blot analysis further corroborated these findings. MPO protein expression exhibited peak levels at 1 day post-SCI, with subsequent significant decreases at 3 days and 7 days [F(3, 8) = 41.91, P = 3.08e-05] (Figure 7C and E) (Figure S1). Similarly, CitH3 expression was maximally upregulated at day 1, followed by progressive reductions at day 3 and day 7 [F(3, 8) = 199.6, P = 7.4e-08] (Figure 7D and F) (Figure S2). These results demonstrate that NETs formation peaks at day 1 after SCI and progressively decreases over time, corroborating our transcriptomic analysis. This finding establishes a critical temporal window during which NETs-mediated inflammatory responses are most pronounced following SCI.

Experimental Validation of Key Genes Expression by qPCR

To validate the expression of key NETs-related genes, we employed a two-tiered qPCR validation strategy. First, we confirmed the neutrophil-specific expression capacity of all seven key NRGs using in vitro isolated neutrophils. In vitro experiments, the expression of all seven identified key NRGs, including Tlr2 [t(5) = 12.87, P = 5.03e-05], Tlr4 [t(5) = 4.93, P = 0.004], Ccl3 [t(5) = 12.22, P = 6.5e-05], Itgam [t(5) = 8.55, P = 0.0004], Casp1 [t(5) = 3.56, P = 0.016], Itgb2 [W = 21.00, P = 0.0313], and Fcgr2b [t(5) = 8.18, P = 0.0004], was significantly upregulated in PMA-stimulated neutrophils compared with the control group (p < 0.05) (Figure 8A and B). This confirms that neutrophils possess the transcriptional machinery to express these genes upon activation.

Figure 8 Experimental validation of the expression of key NETs-related genes by qPCR. (A and B) qPCR analysis of key NETs-related genes in isolated neutrophils under in vitro conditions. (A) Expression of Tlr2, Tlr4, Ccl3, and Itgam, and (B) Casp1, Itgb2, and Fcgr2b in PMA-stimulated neutrophils compared to control groups. The expression of Tlr2, Tlr4, Ccl3, Itgam, Casp1, and Itgb2 and Fcgr2b was significantly upregulated in PMA-treated neutrophils (p < 0.05) (n = 6). (C–I) In vivo qPCR validation of key NRGs in spinal cord tissue at different time points after SCI (Sham, SCI-1d, SCI-3d, and SCI-7d). All genes showed significant upregulation compared to sham group across all time points (p < 0.05) (n = 6). Data are presented as mean ± SD. Statistical analysis was performed using paired t-test (A and B) or one-way ANOVA followed by Tukey’s post hoc test. *P < 0.05, **P < 0.01, ***P < 0.001. Abbreviations: NS, not significant; SCI, spinal cord injury; NETs, neutrophil extracellular traps; PMA, phorbol 12-myristate 13-acetate; qPCR, quantitative polymerase chain reaction; NRGs, NETs-related genes.

Subsequently, we examined the temporal expression dynamics of these genes in spinal cord tissue at different time points following SCI. While this approach captures composite gene expression from multiple cell types within the injury microenvironment, it provides essential validation of the overall temporal patterns identified in our bulk transcriptomic analysis. Importantly, our single-cell RNA sequencing data demonstrated that these genes are predominantly expressed in neutrophils during the acute phase, providing cellular context for interpreting the tissue-level expression patterns. The expression of Tlr2 [Welch’s F(3, 9.279) = 84.60, P = 4.71e-07] remained significantly elevated across all SCI time points compared with sham group (Figure 8C). The expression of Tlr4 [F(3, 20) = 46.57, P = 3.31e-09] and Fcgr2b [F(3, 20) = 56.24, P = 6.29e-10] peaked at 3d post-injury before declining at 7d (Figure 8D and I). The expression of Ccl3 [Welch’s F(3, 9.506) = 124.43, P = 6.17e-08], Itgam [F(3, 20) = 86.25, P = 1.29e-11] and Itgb2 [Welch’s F(3, 10.40) = 62.14, P = 6.18e-07] demonstrated an overall progressive increasing trend from 1d to 7d post-SCI, with sustained significant upregulation compared to sham group throughout all time points (Figure 8E, F and H). Casp1 [F(3, 20) = 82.00, P = 2.06e-11] exhibited highest expression at 1d post-SCI followed by significant reduction at later time points (Figure 8G). All genes showed significant upregulation compared to sham group across all time points (p < 0.05), confirming our bioinformatic findings and supporting the temporal dynamics of NETs-related genes after SCI.

Analysis of NETs-Related TF-mRNA and miRNA-mRNA Regulatory Network in SCI

To further study the transcriptional regulatory mechanisms controlling the expression of DENRGs, algorithms based on TF binding site and gene association position information provided by ENCODE were adopted to construct the regulatory network between TF and DENRGs. A total of 77 connections, containing 27 TFs and 14 DENRGs were uncovered (Figure 9A). CHD1 was screened as a potential important transcription factor regulating NETs in SCI (regulating 5 DENRGs). In addition, EP300, MAX, MXI1, RCOR1, NRF1, ETS1, and IRF4 can regulate 4 DENRGs respectively.

Figure 9 NETs-related transcriptional and post-transcriptional regulatory networks in SCI. (A) TF-mRNA regulatory network of 22 shared DENRGs in SCI. Yellow diamond nodes represent TFs, and cyan circular nodes represent DENRGs. Red arrows indicate regulatory relationships between TFs and target genes. The network comprises 77 connections containing 27 TFs and 14 DENRGs. (B) miRNA-mRNA regulatory network of DENRGs in SCI. Cyan circular nodes represent DENRGs and pink rectangular nodes represent miRNAs. The network consists of 91 nodes and 80 edges. The network visualization reveals complex post-transcriptional regulation of NETs-related genes. The construction of networks provides insights into the complex molecular mechanisms regulating NETs formation following SCI and highlight potential therapeutic targets.

We further explored the regulatory relationship between DENRGs and miRNA through the combined prediction of three miRNA online databases and construct the NETs-related mRNA-miRNA regulatory network to highlight the role of miRNA for this pathological mechanism in SCI. A total of 80 mRNA-miRNA pairs, including 78 miRNA and 13 mRNA were uncovered. Further, the visualization of this network, which comprised 91 nodes and 80 edges, was achieved by Cytoscape (Figure 9B).

Screening of Potential Drug Compounds

Analysis of the DSigDB database indicated that atorvastatin was an important potential therapeutic compound targeting NRGs in SCI (Figure 10A). Other compounds showing significant potential included cholecalciferol, simvastatin, hydrocortisone, and ropivacaine.

Figure 10 Prediction and molecular docking analysis of potential therapeutic compounds targeting NETs-related genes in SCI. (A) Bubble plot showing eight potential therapeutic compounds targeting NETs-related genes in SCI. Atorvastatin emerged as the most significant potential therapeutic agent, followed by cholecalciferol, simvastatin, hydrocortisone, and ropivacaine. (B-H) Molecular docking analyses showing the binding interactions between atorvastatin and the seven key NETs formation-related regulatory proteins. (B) Atorvastatin-ITGAM complex with a binding energy of −7.8 kcal/mol. (C) Atorvastatin-ITGB2 complex with a binding energy of −8.1 kcal/mol. (D) Atorvastatin-TLR2 complex with a binding energy of −9.3 kcal/mol. (E) Atorvastatin-TLR4 complex with a binding energy of −9.3 kcal/mol. (F) Atorvastatin-CASPAS1 complex with a binding energy of −6.7 kcal/mol. (G) Atorvastatin-CCL3 complex with a binding energy of −6.3 kcal/mol. (H) Atorvastatin-FCGR2B complex with a binding energy of −6.6 kcal/mol. Abbreviations: SCI, spinal cord injury; NETs, neutrophil extracellular traps.

Further Molecular docking analysis demonstrated good binding activity between atorvastatin and all seven key NETs formation-related regulatory proteins. In particular, atorvastatin exerted a strong binding activity on four key proteins, with binding energies all lower than −7.0 kcal/mol, including: TLR2 (−9.3 kcal/mol), TLR4 (−9.3 kcal/mol), ITGB2 (−8.1 kcal/mol), and ITGAM (−7.8 kcal/mol) (Figure 10B–H). The docking models showed atorvastatin occupied critical binding pockets within these proteins, suggesting a potential mechanism for inhibiting NETs formation and associated inflammatory processes in SCI.

Discussion

In this study, we employed an integrated bioinformatics approach to comprehensively identify and characterize the role of NETs-related genes in SCI. Through rigorous analysis of multiple transcriptomic datasets and experimental validation, we identified seven key NRGs (Casp1, Ccl3, Fcgr2b, Itgam, Itgb2, Tlr2, and Tlr4) that demonstrate significant temporal differential expression following SCI. The GSVA-based analysis found that the NETs score was highest at 1 day post-injury and then gradually decreased. Analysis of single-cell RNA sequencing further revealed that these key NRGs were predominantly expressed in neutrophils, with highest expression at 1 day post-injury, coinciding with maximal neutrophil infiltration. This prediction was subsequently validated through immunofluorescence and Western blot analyses of established NETs markers (MPO and CitH3), which demonstrated significantly elevated NET formation at the injury epicenter 1 day post-SCI, with progressive attenuation observed at 3 and 7 days post-injury. Furthermore, our drug prediction analysis using the DSigDB database identified several potential therapeutic compounds targeting NET formation in SCI, including atorvastatin, cholecalciferol, simvastatin, hydrocortisone, and ropivacaine. Molecular docking analysis of atorvastatin, as a promising therapeutic candidate, revealed strong binding activity to multiple key NETs formation-related regulatory proteins. These findings provide novel insights into the molecular mechanisms underlying NETs formation in SCI pathophysiology and highlight potential therapeutic targets for mitigating secondary inflammatory damage.

Our findings are consistent with and extend previous studies investigating NETs formation in SCI. Recently, Reid et al demonstrated that neutrophil-mediated extracellular trap formation peaks within the first day after SCI in a murine contusion model, which aligns remarkably well with our temporal analysis showing highest NETs scores at 1 day post-injury.¹⁷ While their study employed ELISA-based quantification of CitH3-DNA and MPO-DNA complexes, our quantitative Western blot and immunofluorescence approaches yielded consistent temporal patterns, which further corroborated our current conclusions. Importantly, Reid et al demonstrated that early administration of recombinant human DNase improved long-term functional recovery and tissue sparing, supporting the therapeutic potential of targeting NETs that we identified through our drug prediction analysis.

Earlier studies have established the pathological role of NETs in SCI but lacked comprehensive molecular characterization. Feng et al demonstrated that NETs exacerbate secondary injury via promoting neuroinflammation and blood-spinal cord barrier disruption in rat SCI models, and showed therapeutic efficacy of DNase1 treatment.¹⁵ However, their study primarily focused on the pathological effects of NETs without systematic identification of key regulatory genes. Our study advances this understanding by identifying seven key NRGs and their coordinated regulatory networks, providing specific molecular targets beyond the general anti-NET approaches used in previous studies. Similarly, Michel-Flutot et al characterized extracellular trap formation following cervical SCI, demonstrating NETs presence through histological analysis, but did not explore the underlying gene regulatory mechanisms that we have systematically characterized.¹³

Recent therapeutic investigations have explored various approaches to modulate NETs in SCI. Morishima et al showed that mesenchymal stem cell-derived exosomes alleviate SCI by regulating NET formation through exosomal miR-125a-3p, highlighting the role of miRNA regulation that complements our comprehensive miRNA-mRNA network analysis.³² Additionally, Diao et al demonstrated that CD47-blocking antibody interferes with NETs formation after SCI to reduce spinal cord edema, providing another therapeutic avenue.¹⁶ While these studies focused on specific interventions, our systematic drug prediction analysis identified multiple potential therapeutic compounds, particularly atorvastatin, which offers advantages over previous approaches due to its established safety profile and FDA approval status. Despite these advances in understanding NETs in SCI, previous studies have primarily focused on descriptive characterization or single therapeutic interventions. Our study provides the first comprehensive molecular framework by systematically identifying key regulatory genes, constructing regulatory networks, and predicting multiple therapeutic targets. This systematic approach bridges the gap between mechanistic understanding and therapeutic development, offering a more complete picture of NETs regulation in SCI than previously available. Furthermore, our identification of specific gene targets provides precision medicine opportunities that were not evident in earlier descriptive studies, potentially enabling more targeted and effective therapeutic interventions.

Our integrated analysis identified seven key NRGs (Casp1, Ccl3, Fcgr2b, Itgam, Itgb2, Tlr2, and Tlr4) that play crucial roles in the regulation of neutrophil function and NETs formation following SCI. Tlr2 and Tlr4, pattern recognition receptors predominantly expressed on neutrophils, function as primary sensors of damage-associated molecular patterns released after SCI.³³ Activation of these receptors initiates downstream signaling cascades that promote neutrophil activation, migration, and NETs formation through NF-κB and MAPK pathways.^16,34,35 In experimental SCI models, inhibition of Tlr4 has been shown to effectively reduce neutrophil infiltration and NET formation and improve functional recovery.¹⁶ Itgam and Itgb2, components of the Mac-1 integrin complex (CD11b/CD18), mediate neutrophil adhesion and migration to injury sites, while also contributing to reactive oxygen species generation necessary for NETs formation.^36,37 In the context of SCI, these integrins mediate neutrophil adhesion to the vascular endothelium and subsequent extravasation into the injured spinal cord, and their targeted blockade could limit post-injury neutrophil accumulation and further exacerbation of tissue damage. Although Fcgr2b, an inhibitory Fc-γ receptor, has been limitedly studied in SCI, it can regulate CNS immune complex-mediated inflammation as well as the threshold of neutrophil activation and the intensity of inflammatory response.³⁸ Its dysregulation may aggravate tissue damage through excessive release of NETs.³⁹ Ccl3 is a potent chemoattractant that coordinates the recruitment of neutrophils to the site of injury. It is particularly robustly expressed in neutrophils during the acute phase of SCI, consistent with its role in promoting neutrophil chemotaxis and activating inflammatory responses.⁴⁰ Additionally, Casp1 activation in the NLRP3 inflammasome pathway promotes IL-1β maturation and pyroptosis, thereby contributing to NET formation.⁴¹ After SCI, Casp1 activation leads to the production of IL-1β and IL-18, which become potent inflammatory cytokines that activate neutrophils and NET formation, further exacerbating the injury.⁴²

Based on our single-cell RNA sequencing analysis, we observed distinctive spatial expression patterns of the seven key NRGs across various cell populations in the injured spinal cord. Most notably, these genes exhibited predominant expression in neutrophils, with Tlr2, Itgb2, Itgam, Fcgr2b, and Ccl3 showing particularly robust expression in the neutrophil cluster. Temporal analysis revealed that expression of these genes in neutrophils peaked at 1 day post-injury, coinciding with maximal neutrophil infiltration and the highest NETs scores determined by GSVA analysis. This temporal concordance strongly suggests that the acute phase following SCI represents a critical window for NETs formation, potentially contributing to secondary inflammatory damage. The progressive reduction in both neutrophil numbers and gene expression levels at 3 and 7 days post-injury further supports this notion and aligns with our protein-level validation of NETs markers. Interestingly, while neutrophils exhibited the most pronounced expression of these genes, we also observed significant expression in monocytes and macrophages, suggesting potential intercellular cooperative mechanisms in NETs formation and regulation. Monocytes and macrophages may contribute to the inflammatory microenvironment by secreting cytokines and chemokines that modulate neutrophil recruitment, activation, and NETs release.^43,44 This cellular heterogeneity in the expression of NRGs highlights the complex cellular interplay during the inflammatory cascade following SCI. The spatial-temporal expression patterns observed in our single-cell analysis provide valuable insights into the complex cellular mechanisms that drive NETs formation after SCI, suggesting that therapeutic interventions targeting these pathways may be most effective when administered during the acute phase of injury, particularly within the first 24 hours.

Our findings regarding the temporal dynamics of NETs formation in SCI have critical clinical implications. The observation that NETs formation peaks at 1 day post-injury and subsequently decreases consistent with a progressive pathological process of secondary injury in SCI. This early peak coincides with the period of significant neutrophil infiltration and precedes the macrophage-dominated inflammatory response that typically manifests 3–7 days post-injury.^8,45 The 24-hour timepoint represents a critical window when the blood-spinal cord barrier disruption is the most severe, allowing neutrophil extravasation and subsequent NETs formation at the injury epicenter.⁴⁶ This temporal pattern suggests that NETs contribute significantly to the early phase of secondary injury cascade, potentially exacerbating tissue damage through the release of proteolytic enzymes, reactive oxygen species, and cytotoxic histones.^11,47,48 The identification of this temporal window offers valuable insights for developing targeted therapeutic interventions. Specifically, pharmacological inhibition of NETs formation within the initial 24 hours post-injury could significantly attenuate secondary inflammatory damage and improve functional outcomes. Furthermore, the rapid elevation of NETs markers during this timepoint suggests their potential utility as early diagnostic biomarkers for SCI severity assessment, potentially enabling early stratification of patients for appropriate interventions.

An important finding of our study is the distinct temporal expression patterns observed across different experimental approaches. While our single-cell analysis suggested that neutrophil-specific expression of key NRGs was most prominent at 1dpi, our bulk transcriptomic analysis (Figure 4G) and qPCR validation (Figure 8C–I) revealed more complex temporal dynamics, with several genes (Itgam, Itgb2, Fcgr2b, Tlr2, Tlr4) showing peak or sustained expression at 3dpi and 7dpi. These discrepancies likely reflect the fundamental differences between analytical approaches and sample compositions. The bulk RNA-seq data represents the collective gene expression of all cell types within the injured spinal cord tissue, where the signal from neutrophils may be diluted by other abundant cell populations such as microglia, macrophages, and astrocytes that also express these genes. In contrast, the single-cell analysis specifically captured neutrophil-intrinsic expression patterns, revealing the temporal dynamics within this specific cell population. Additionally, while neutrophil infiltration peaks at 1dpi, the sustained expression of NRGs at later timepoints in bulk tissue likely reflects the continued activation of resident and infiltrating myeloid cells, including macrophages and microglia, which also contribute to NETs-related processes. The qPCR validation using whole spinal cord tissue further supports this interpretation.

The integration of single-cell RNA sequencing with bulk transcriptomic analysis represents a powerful approach for understanding NETs formation in SCI. Single-cell analysis provides unprecedented resolution to dissect cell-type-specific gene expression patterns, allowing us to identify the precise cellular sources of NRGs and track their temporal dynamics without interference from other cell populations. This approach is particularly valuable for studying NETs formation, as neutrophils represent a relatively small and transient population within the injured spinal cord, making their specific contributions difficult to discern from bulk tissue analysis alone. The complementary use of bulk transcriptomic analysis provides essential context by capturing the overall inflammatory landscape and the collective contribution of multiple cell types to NETs-related processes. This dual approach allows for a more comprehensive understanding of how neutrophil-specific responses integrate with broader tissue-level inflammatory cascades. Moreover, the combination of these approaches with experimental validation through immunofluorescence, Western blot, and qPCR creates a robust framework for mechanistic understanding that bridges molecular discoveries with functional outcomes. This multi-layered analytical strategy is essential for developing targeted therapeutic interventions, as it identifies both the specific cellular targets (neutrophils) and the broader inflammatory context that must be considered for effective NETs modulation in SCI.

Our comprehensive transcriptional regulatory network analysis revealed intricate mechanisms regulating NETs formation in SCI. Notably, CHD1 emerged as a pivotal transcription factor regulating multiple NRGs, suggesting its central role in orchestrating neutrophil activation and NETs formation following SCI. CHD1, a chromatin remodeling protein with ATP-dependent helicase activity, likely controls access to regulatory regions of NRGs, facilitating coordinated transcriptional responses.^49,50 Additionally, transcription factors EP300, MAX, MXI1, RCOR1, and NRF1 demonstrated significant regulatory potential, collectively forming a complex transcriptional landscape that governs neutrophil function post-injury. The miRNA-mRNA network analysis identified extensive regulatory interactions, with numerous miRNAs targeting key NRGs, indicating sophisticated post-transcriptional control mechanisms. These miRNAs may specifically regulate neutrophil activation at different stages of SCI.^32,51,52 The exploration of these multi-layered regulatory networks significantly enhances our understanding of the molecular mechanisms underlying NETs formation in SCI and offers promising targets for therapeutic intervention.

The drug prediction analysis reveals promising clinical translational potential for several compounds targeting NETs formation in SCI. Notably, atorvastatin emerged as the top-ranked candidate with good binding activity to key NETs formation-related regulatory proteins. While previous studies have discovered that atorvastatin can improve functional recovery and reduces tissue damage following SCI, our current findings unveil a novel mechanism by which it may exert its neuroprotective effects through the regulation of NETs formation.^53–55 The pleiotropic effects of atorvastatin extend beyond its lipid-lowering properties, demonstrating significant anti-inflammation through inhibition of neutrophil activation and NETs formation. Recent evidence has shown that atorvastatin effectively suppresses NETs formation and reduces inflammation in septic lung injury via the ERK/NOX2 signaling pathway, improving survival rates.⁵⁶ And, high-dose atorvastatin treatment significantly reduces levels of NET-related proteins including CitH3, MPO, and neutrophil elastase.⁵⁷ Moreover, atorvastatin can attenuate TLR2/4-mediated NF-κB signaling and reduces integrin-dependent neutrophil adhesion and migration to the injury site, thereby limiting excessive NETs release and secondary tissue damage.^58,59 These effects were associated with reduced inflammatory factors and ameliorated tissue damage, potentially suggesting that this drug could repair SCI through an anti-NET mechanism. Other predicted compounds also showed significant promise. Cholecalciferol, through vitamin D receptor activation, may suppress neutrophil hyperactivation and NETs formation by modulating calcium signaling and inflammatory cytokine production.⁶⁰ Similarly, studies have shown that simvastatin can effectively prevent the formation of NETs by inhibiting the PAD4 or Mac-1 pathway.^61,62 Hydrocortisone and ropivacaine may provide anti-inflammatory and analgesic effects by attenuating neutrophil activation and NET-mediated damage.^63–65 These findings suggest that the predicted drugs could be used to target NETs formation, offering a practical therapeutic strategy to reduce secondary damage and enhance recovery after SCI.

This study integrates bioinformatics analysis with experimental validation to identify key NRGs in SCI. By combining multiple transcriptomic datasets, WGCNA and single-cell sequencing, we have established a comprehensive understanding of the regulatory mechanisms of NETs formation following SCI. The identification of seven key regulatory genes and their temporal expression patterns provide valuable insights into the molecular mechanisms of neutrophil-driven secondary injury after SCI. In addition, the drug prediction analysis offers practical translational value by identifying compounds that could potentially target NETs formation in SCI. However, several limitations warrant consideration. Firstly, this study predominantly relied on transcriptomic data and in vivo validation in murine models, which may not fully recapitulate the complexity of human SCI. Meanwhile, an important limitation is the potential influence of sex and age differences between datasets and experimental validation. Our experimental validation used 8-week-old female C57BL/6 mice. While this is largely consistent with the single-cell data coming from 8–10-week-old female mice, the sex and age information for the bulk RNA sequencing datasets is not available in the original database documentation. Given that age and sex have been shown to dramatically affect neutrophil function, future studies should systematically investigate the impact of these factors on NETs formation and regulation in SCI.⁶⁶ Secondly, in this study, our identification of NETs using MPO+CitH3+ co-localization, though representing a well-established and widely accepted methodology in the field, could potentially benefit from additional neutrophil-specific markers such as Ly6G to further enhance cellular specificity. Future studies incorporating multiple neutrophil-specific markers could provide additional validation of these findings. Thirdly, our qPCR validation was performed on whole spinal cord tissue rather than isolated neutrophils from the injured tissue, which limits the cellular specificity of our gene expression validation. While direct neutrophil isolation from injured spinal cord across multiple time points presents significant technical challenges due to low cell numbers and potential isolation-induced transcriptional changes, this approach would provide more definitive cell-specific validation. Our current validation strategy combining in vitro neutrophil studies with tissue-level analysis, supported by single-cell RNA sequencing data, provides reasonable evidence but represents a methodological compromise. Fourthly, our analysis primarily focused on the acute phase (1–7 days) post-injury, necessitating further investigation into long-term NETs dynamic modulation during chronic SCI. Furthermore, while atorvastatin showed promising potential as a therapeutic agent, its effect, mechanism and optimal administration regimen on NETs formation after SCI need to be further confirmed through continuous experiments. Notwithstanding these limitations, our current study provides a solid foundation for targeting NETs formation as a novel therapeutic strategy for SCI treatment.

Conclusion

Our study provides comprehensive insights into the temporal dynamics and molecular regulation of NETs in SCI. We identified seven key NETs-related genes (Casp1, Ccl3, Fcgr2b, Itgam, Itgb2, Tlr2, Tlr4) that play critical roles in NET formation following SCI. Temporal analysis revealed peak NET formation at day 1 post-injury with progressive attenuation at days 3 and 7, establishing a critical therapeutic window. Single-cell transcriptome analysis demonstrated predominant expression of these key genes in neutrophils during the acute phase, coinciding with maximal neutrophil infiltration. Regulatory network analysis identified CHD1 as an important transcription factor governing NET formation while revealing complex miRNA-mediated regulatory mechanisms. Drug prediction analysis identified atorvastatin as a promising therapeutic candidate, with molecular docking demonstrating strong binding to multiple NET formation-related proteins. These findings elucidate the molecular mechanisms underlying NETs formation in SCI and highlight potential therapeutic targets for mitigating Neutrophil-mediated secondary inflammatory damage, potentially improving functional outcomes in SCI patients.