Introduction
Traumatic brain injury (TBI) represents a critical global public health concern, characterized by structural and functional disruption of the brain due to external mechanical forces. Affecting approximately 69 million individuals globally each year,1-3 TBI creates a staggering socioeconomic burden, with profound multidimensional consequences that extend beyond acute clinical manifestations. Although the acute phase of brain trauma is well documented, the long-term pathological outcomes, especially those following mild injuries, remain largely underappreciated in clinical settings. Recent clinical and epidemiological data are redefining TBI not merely as a single event but as an evolving pathology that drives sustained cognitive and motor impairments, thereby severely degrading patient well-being over time.3 Notably, growing data implicate TBI as a risk factor for delayed neurodegenerative processes, underscoring its dual nature as both an acute injury and a chronic disease entity. Current TBI management remains predominantly supportive, focusing on intracranial pressure modulation, oxygenation, and surgical interventions.3-5 However, effective pharmacological treatments for TBI remain limited.6
TBI represents a progressive disorder characterized by a primary mechanical insult and a subsequent secondary injury phase involving excitotoxicity, oxidative stress, and neuroinflammatory cascades. Microglia, the primary immune effectors of the central nervous system, exhibit stage-dependent functional plasticity. As microglia polarize toward a pro-inflammatory state, they may undergo metabolic reprogramming characterized by reduced oxidative phosphorylation and increased reliance on glycolysis. Identifying the molecular drivers of this “immunometabolic shift” is critical for developing neuroprotective interventions.7-14 Therefore, this study aimed to integrate single-cell and bulk transcriptomic datasets to identify microglia-associated regulatory genes (MRGs) and characterize their potential roles in post-TBI neuroinflammatory and immunometabolic dysregulation.
Materials and methods
Gene data and clinical information collection
TBI can cause structural damage and functional impairment in the neocortex, the outer layer of the brain, via direct impact or secondary injury mechanisms, and the highly vulnerable hippocampus is also susceptible to neuronal apoptosis and neurogenesis inhibition induced by shear injury, ischemia, and hypoxia; accordingly, individuals with TBI in clinical settings frequently exhibit advanced cognitive dysfunction as well as core deficits in memory and spatial navigation. To explore the commonly affected genes in murine hippocampal and neocortical tissue samples, the following datasets were selected for subsequent analysis.
All transcriptomic datasets analyzed in this study were retrieved from the public Gene Expression Omnibus (GEO) repository (https://www.ncbi.nlm.nih.gov/geo/ ). The training set was derived from the GSE58485 dataset (GPL1261), which included 3 TBI samples (GSM1412408-GSM1412410, representing a moderate-to-severe focal brain injury model at 3 days post-injury) and 3 control samples (GSM1412411-GSM1412413, uninjured neocortical tissues).15 We confirm that all six selected samples were exclusively derived from untreated, wild-type adult male C57BL/6 mice. These samples were utilized for differential expression analysis, Gene Set Enrichment Analysis (GSEA), and immune cell infiltration analysis.
To corroborate our findings, an independent validation cohort was established using transcriptomic profiles derived from the Clariom D Mouse/Murine Assay microarray (GSE242025). This dataset features five pairs of ipsilateral injured (utilizing a moderate-to-severe controlled cortical impact model at 7 days post-injury) and contralateral comparator cerebral tissues derived from 8- to 12-week-old male and female C57BL/6 mice, thereby ensuring a paired statistical design.16 The data were employed for validation of expression levels. Additionally, the single-cell dataset GSE101901 (GPL17021) was downloaded, containing transcriptomic profiles from 3 TBI samples (derived from a mild-to-moderate concussive brain injury model at 24 h post-injury) and 3 control samples (hippocampal tissues).17 This dataset was used for single-cell analysis.
Single-cell quality control
Quality control and downstream single-cell processing were performed using the ‘Seurat’ suite (version 5.1.0) in R.18 To ensure matrix robustness, strict filtering parameters were applied: transcripts expressed in fewer than 3 individual cells were discarded; droplets exhibiting feature counts below 200 or exceeding 2,500, as well as total UMI counts over 6,000, were excluded; and elements with a mitochondrial read ratio surpassing 20% were similarly filtered out. Following filtration, the remaining high-quality cellular matrices underwent global-scaling normalization. We subsequently isolated transcripts demonstrating substantial cell-to-cell expression variance to define the highly variable genes (HVGs) for downstream processing. HVGs were visualized in a scatter plot, with the top 10 genes exhibiting the highest variation. The data were normalized again through the ScaleData function. HVGs were then evaluated by principal component analysis. Appropriate dimensionality was determined by evaluating JackStraw statistical thresholds and identifying the distinct inflection point via elbow plot heuristics. Cells were clustered using a graph-based clustering approach (via the FindNeighbors and FindClusters functions implementing the Louvain algorithm) at a resolution of 0.5. Subsequently, the uniform manifold approximation and projection method was employed for nonlinear dimensionality reduction and two-dimensional visualization.
Cluster-specific marker identification was conducted by calculating the area under the receiver operating characteristic curve, applying stringent cutoff criteria including a minimum cellular detection rate of 10%, a log-fold change exceeding 0.25, and an adjusted significance level below 0.05. Then, cluster-specific differentially expressed genes were compared with marker genes from the literature and the CellMarker database (http://117.50.127.228/CellMarker/ ) to identify TBI marker genes (TMGs).17 Furthermore, the cells were annotated based on TMGs. The sample-level proportional abundances of these cell types were evaluated between the TBI (n = 3 individual mice) and control (n = 3 individual mice) groups in the GSE101901 dataset. Changes in cell proportions were reported strictly as descriptive trends to identify key cells, avoiding pseudoreplication from single-cell-level statistical testing.
Pseudotime analysis
Key cells were re-clustered using the R package “Seurat” (v 5.1.0) based on HVGs.18 To reconstruct the dynamic transcriptional-state progression of relevant cellular clusters, we deployed the Monocle computational pipeline (v 2.28.0) to infer pseudotemporal trajectories.19 Specifically, ordering genes were identified using the Monocle 2 function differentialGeneTest with a model formula based on cell subclusters (i.e., fullModelFormulaStr = “~Cluster”). Genes meeting a Benjamini–Hochberg (BH)-adjusted P-value < 0.05 were subsequently selected for trajectory construction via the setOrderingFilter function. Finally, dimensionality reduction was performed using the reduceDimension function with the DDRTree algorithm, followed by the orderCells function to align cells along the pseudotime trajectory. Crucially, these trajectories were inferred on the merged dataset (comprising both TBI and control samples) to enable direct cross-condition mapping. To establish the biological directionality, the pseudotime root was designated as the state predominantly populated by control cells, representing the homeostatic baseline prior to injury. No formal sensitivity analyses regarding alternative root choices or algorithmic parameters were performed. Subsequently, the pseudotemporal trajectory of each key cell type was identified, and expression trends of TMGs in the inferred pseudotemporal trajectory were also detected.
Cell-cell communication analysis
Intercellular signaling networks and ligand-receptor interactomes were computationally inferred using CellChat.20 The analysis was performed on the merged single-cell RNA sequencing (scRNA-seq) dataset, and inferred interactions with a communication strength greater than 0.1 were retained for visualization and descriptive interpretation. The predicted ligand-receptor pairs among different cell types were visualized using interaction networks and bubble plots. Because the analysis was conducted on the merged dataset, it was not used to test between-group differences in communication strength.
Identification of microglia-associated gene modules using hdWGCNA
We adopted an hdWGCNA strategy specifically to unearth modular transcriptional networks driving microglial behavior. The analysis utilized scRNA-seq data from all samples. To mitigate the extreme sparsity inherent to single-cell matrices, microglial transcriptomes were pooled into robust metacells—each constructed from 10 k-nearest neighbors—prior to network assembly.21
We iteratively tested multiple integer soft powers and selected the minimum soft thresholding power of 6 (which achieved a scale-free topology fit index of approximately 0.8) to guarantee that the resulting co-expression matrix adhered to scale-free distribution principles. A weighted co-expression network was subsequently constructed with the ConstructNetwork function based on the selected soft threshold. Transcripts were dynamically grouped into functional modules relying on average linkage hierarchical clustering applied to a topological overlap-based dissimilarity matrix. Pairwise correlations between modules were evaluated via Spearman correlation analysis (P < 0.05).
Module eigengenes (MEs) were calculated as the first principal component of each module’s gene-expression matrix. Batch effects were corrected using the R package “harmony” (v 1.2.1).22 Gene connectivity within modules was quantified by calculating the module membership (kME), defined as the correlation between individual genes and their corresponding MEs, using the ModuleConnectivity function. To select the most representative core regulators within each co-expression network and minimize background transcriptional noise, the top 10 genes with the highest kME values in each module were selected, collectively defining MRGs for downstream analyses.
Differential expression analysis and candidate gene screening
For both bulk transcriptomic datasets (GSE58485 and GSE242025), the pre-processed, base-2 logarithm (log2)-transformed expression matrices were retrieved directly from the GEO database. Background correction and quantile normalization had been executed according to the standard platform protocols prior to matrix deposition. Probe-to-gene mapping was performed utilizing the corresponding platform annotation files (GPL1261 and GPL20258, respectively), and redundant probes mapping to a single gene were aggregated by calculating the mean expression value. Because all samples within each respective dataset were generated in a single experiment, no additional batch-effect correction algorithms were required. Bulk transcriptomic variances between the independent TBI (n = 3 independent animal tissues) and control (n = 3 independent animal tissues) groups in the GSE58485 dataset were statistically evaluated utilizing the linear modeling capabilities of the “limma” package (v 3.58.1).23 An independent design matrix was applied, and redundant probes were averaged after mapping to specific genes. Significance boundaries for nominally differentially expressed genes were defined using an exploratory threshold of a raw P-value < 0.05 coupled with an absolute log2 fold-change > 0.5. We additionally calculated the false discovery rate (FDR) using the BH procedure to account for multiple testing. To further visualize the results, volcano plots were generated using the R package “ggplot2” (v 3.5.1),24 and the top 10 upregulated and top 10 downregulated nominal DEGs were labeled. Additionally, a heatmap visualizing the expression patterns of these top 20 dysregulated genes was generated using the R package “ComplexHeatmap” (v 2.16.0).25
Furthermore, to identify MRGs differentially expressed in TBI, the intersection of DEGs and MRGs was determined. The intersecting genes were designated as candidate genes in TBI. The R package “VennDiagram” (v 1.7.3) was employed to visualize the intersections.26
Functional enrichment
To interrogate the biological relevance of the candidate targets, we deployed the clusterProfiler algorithm to map significant ontological categories (Gene Ontology; GO) and biochemical cascades (Kyoto Encyclopedia of Genes and Genomes; KEGG), using an adjusted P-value threshold of < 0.05.27 Specifically, GO enrichment analysis categorized genes into three functional domains: biological processes (BP), cellular components (CC), and molecular functions (MF), while KEGG analysis identified significant biological pathways. The pathways were ranked according to the counts from high to low.
Protein-protein interaction (PPI) network construction
A PPI network was constructed using the STRING database, with interactions retained at a minimum medium-confidence score of 0.4.28 After obtaining the interaction information, isolated targets were removed. The resulting interaction topology was imported into the Cytoscape graphical environment for advanced network modularization and visualization.29 Based on topological analysis, candidate genes with the highest degree centrality within the network were identified as hub genes for downstream analysis.
Expression validation
To cross-validate the expression levels of the hub genes, we analyzed the GSE242025 validation dataset, which features a strictly paired design comprising ipsilateral injured and contralateral comparator tissues from the same animals (n = 5 paired mice, 10 tissue samples in total). Statistical differences in the extracted expression values of individual hub genes were evaluated using standard paired Student’s t-tests (or paired Wilcoxon signed-rank tests for nonnormally distributed data). This paired approach directly matched the ipsilateral and contralateral samples by individual animal ID to account for within-subject baseline variance, with exact P-values reported. Hub genes that exhibited significant nominal differential expression (P < 0.05) with consistent trends across datasets were identified as key genes.
GSEA
Next, GSEA was employed to explore global biological functions and pathway alterations in TBI. Given that the GSE58485 training dataset is limited to 3 TBI and 3 control samples, genes were preranked based on the moderated t-statistic derived from the limma group-level comparison (TBI versus control). GSEA was then conducted using the R package “clusterProfiler” (v 4.8.3) with KEGG canonical pathway gene sets from the C2 collection of the Molecular Signatures Database (MSigDB) (https://www.gsea-msigdb.org/gsea/msigdb ), mapped to Mus musculus.27 Enrichment was considered significant with |normalized enrichment score (NES)| > 1, FDR < 0.25, and P < 0.05.
To quantify dynamic fluctuations in hallmark pathway activities across different tissue conditions, we executed gene set variation analysis (GSVA) methodologies mapped against the standardized mouse hallmark signatures provided by MSigDB.30 Sample-level GSVA scores were calculated for each tissue. Subsequently, differential analysis of the group-level GSVA scores between the independent TBI (n = 3) and control (n = 3) groups was conducted using the “limma” package (v 3.58.1),23 applying an independent design matrix. Both nominal P-values and BH-FDR were reported to identify pathways showing significant alterations.
Immune infiltration analysis
The proportional abundance of leukocyte subpopulations in the GSE58485 dataset, comprising three TBI and three control samples, was estimated using CIBERSORT with a mouse-specific ImmuCC signature matrix.31 Given the limited sample size, between-group differences in estimated immune-cell proportions were summarized descriptively rather than subjected to confirmatory hypothesis testing. Spearman correlation coefficients were calculated using the R package “psych” (v 2.4.3) to explore associations among immune-cell signatures and between key-gene expression and immune-cell signatures.32 These coefficients were treated as exploratory effect-size estimates without confirmatory P-value thresholds or confidence-interval interpretation. The results were visualized using the R package “ggplot2” (v 3.5.1).24
In silico prediction of transcription factor binding sites (TFBS)
To explore the transcriptional relationships within the Ccl4/Lgals3/Egr1 expression pattern, we performed motif scanning on the promoter regions of target genes. The specific DNA-binding consensus sequence matrix for the transcription factor EGR1 was sourced from the JASPAR database (Matrix ID: MA0162.1). Regulatory genomic architectures—capturing the region spanning 2 kilobases upstream to 0.5 kilobases downstream of the respective transcription start sites—were securely extracted based on the mouse reference genome assembly mm10 (GRCm38). The FIMO (Find Individual Motif Occurrences) tool from the MEME Suite (v 5.5.9) was utilized to scan these promoter sequences for potential EGR1 motif occurrences. A prespecified threshold of P < 0.0001 was applied to identify predicted motif occurrences.
Drug prediction
To explore potential therapeutic agents targeting TBI, homology conversion of mouse genes to human orthologs was performed using the R package “homologene” (v 1.4.68.19.3.27) (https://CRAN.R-project.org/package=homologene ). Key genes were converted to their human counterparts. Drug-gene interactions were subsequently queried from the Drug Gene Interaction Database (DGIdb) (https://www.dgidb.org ). Subsequently, a gene-drug interaction network was constructed based on the retrieved data, and network visualization was implemented using Cytoscape software (v 3.10.2). Interactions were filtered to include only direct pharmacological relationships documented in DGIdb.
Real-time quantitative polymerase chain reaction (RT-qPCR) experimental verification
To validate key gene expression in vivo, 10 six-to-eight-week-old male C57BL/6 mice (22 ± 2 g) were obtained from GemPharmatech (Chengdu, China). The mice were randomly assigned by simple random drawing to two experimental groups, Sham and TBI (n = 5 per group). The animals were maintained in a standard specific pathogen-free environment featuring a 12-h light-dark cycle at a constant temperature of 22 ± 1 °C and provided unrestricted access to food and water. All experimental protocols adhered to the National Institutes of Health (NIH) Guide for the Care and Use of Laboratory Animals and were approved by the Ethics Committee of the General Hospital of Western Theater Command (Approval No. 2023ky119-1). All animal experiments were designed and reported in strict accordance with the ARRIVE 2.0 guidelines. No formal a priori statistical power calculation was performed. The sample size of five mice per group was empirically determined and was consistent with that used in the referenced independent validation cohort (GSE242025).16 To minimize experimental bias, investigators performing tissue collection and RT-qPCR data analysis were strictly blinded to group allocation. The operating surgeon was not blinded because the surgical procedures inherently differed between groups. This specific paradigm, utilizing the parameters detailed below, represents a moderate-to-severe focal TBI model.33 We selected this severe controlled cortical impact model because it generates highly reproducible contusions and robust secondary neuroinflammatory cascades, providing an optimal pathological environment to study pronounced microglial activation and macrophage infiltration. Briefly, mice were anesthetized via intraperitoneal injection of 2.5% tribromoethanol (250 mg/kg; purchased from Sigma-Aldrich, USA). During the perioperative period, deep anesthesia was verified by the absence of the pedal withdrawal reflex, and body temperature was maintained at 37 °C using a homeothermic heating pad. Routine postoperative pharmacological analgesia was intentionally omitted to prevent severe confounding effects from systemic anti-inflammatory agents on the microglial transcriptomic and immunometabolic profiles, which was explicitly approved by the institutional ethics committee. Following anesthesia induction, mice were secured in a stereotaxic frame. Following anesthesia, a unilateral 3-mm cranial window was carefully drilled over the right parietal cortex, meticulously avoiding any disruption to the underlying dural integrity. Sham-operated mice underwent identical anesthesia and craniotomy without the subsequent cortical impact. For the TBI group, the mechanical insult was precisely delivered via an automated electromagnetic impactor, driving a 2.5-mm flat-tipped rod perpendicularly into the exposed cortex at 3.5 m/s with a tissue deformation depth of 1.0 mm and a dwell duration of 100 milliseconds. The incision was sutured before recovery. At 7 days post-injury, mice were euthanized by an intraperitoneal overdose of sodium pentobarbital (150 mg/kg). Deep anesthesia and death were confirmed by the complete cessation of respiration and heartbeat, alongside the absence of pedal and corneal reflexes. Following this confirmation, cervical dislocation was performed as a secondary physical method to ensure death before the ipsilateral pericontusional cortex 2 mm from the lesion core was microdissected. Total RNA was extracted using TRIzol reagent (Vazyme, China), and RNA concentration was evaluated using a fluorescence detection instrument. cDNA was synthesized using the HP All-in-one qRT Master Mix II (YunGen Biosciences, China) for subsequent qPCR assays using standard protocols. The qPCR assay was performed with a CFX Connect Thermal Cycler (Bio-Rad, USA) using the 2× Universal Blue SYBR Green qPCR Master Mix (Servicebio, China). The thermocycling conditions consisted of an initial denaturation at 95 °C for 2 min, followed by 40 cycles of 95 °C for 10 s and 60 °C for 30 s. Technical triplicates and no-template controls were performed for all samples. Primer amplification efficiencies (90–110%) and the presence of a single melt curve peak were validated for all targets. Gapdh was selected as the endogenous control due to its widespread establishment as a reliable reference gene in murine TBI models. Its suitability in our specific cohort was verified retrospectively, as the raw cycle threshold (Ct) values of Gapdh exhibited no significant statistical variance between the Sham and TBI groups (P > 0.05). The relative quantification of mRNAs was computed using the 2−ΔΔCT method. The sequences of all primers can be found in Table 1. Once the RT-qPCR results were generated, they were first exported to Excel and subsequently imported into GraphPad Prism 10 (https://www.graphpad.com/) for statistical analysis and visualization. Statistical significance for qPCR data was determined using the two-sided independent Student’s t test (for normally distributed data evaluated by the Shapiro-Wilk test) or the Mann-Whitney U test, with exact P-values reported. Results are presented as mean ± standard deviation, and no animals were excluded from the final analysis.
| Gene | Forward primer (5′–3′) | Reverse primer (5′–3′) | Target accession number |
|---|
| Ccl4 | TTCCTGCTGTTTCTCTTACACCT | CTGTCTGCCTCTTTTGGTCAG | NM_013652 |
| Lgals3 | GTTGCGGTCAACGATGCTC | TGATCCCCAGTTGGCTGATTT | NM_010705 |
| Egr1 | CCGAGCGAACAACCCTATGA | AGGCTGAAAAGGGGTTCAGG | NM_007913 |
| Gapdh | TGTGTCCGTCGTGGATCTGA | GAGTTGCTGTTGAAGTCGCA | NM_008084 |
General statistical considerations
All statistical analyses were evaluated using two-sided tests unless otherwise specified. Individual animals or independent tissue samples were strictly treated as the biological statistical units; individual cells in scRNA-seq were not treated as independent biological replicates. For the GSE242025 dataset, which utilized ipsilateral and contralateral tissues from the same animals, the paired design was accounted for in the statistical analyses. When conducting multiple comparisons across genes or pathways, the BH procedure was applied to control the FDR, with both nominal P-values and BH-adjusted P-values reported where applicable. A predefined threshold of an adjusted P < 0.05 was considered statistically significant, except for prespecified exploratory analyses using nominal P-values. All graphical error bars represent the standard deviation or standard error of the mean, and the exact sample sizes (n) and specific statistical tests are explicitly annotated in the respective figure legends and method subsections. A comprehensive summary of all statistical analyses, including sample sizes, specific tests, and correction methods, is provided in Supplementary Table. 1.
Results
Single-cell transcriptomic profiling of the TBI ecosystem
Following stringent quality control protocols (Supplementary Fig. 1), our final analytical scRNA-seq matrix consisted of 17,597 genes across 7,173 high-quality cells. Subsequently, 2,000 HVGs were selected, including Ttr, Hbb-bs, Hbb-bt, and Ccl12. By leveraging the top 17 statistically significant principal components, the cellular population was objectively partitioned into 14 distinct transcriptional clusters. Subsequently, seven major cell types were annotated: oligodendrocytes, astrocytes, microglia, neurons, endothelial cells, oligodendrocyte precursor cells, and mural cells based on the TMGs (Fig. 1). Comparative analysis between TBI and control samples highlighted apparent descriptive shifts in the proportional abundance of astrocytes, microglia, and neurons between the TBI and control groups at the sample level. While astrocytes showed an apparent descriptive increase in abundance, microglia were prioritized for downstream network analysis because of their established role as innate immune cells in post-TBI neuroinflammation and their relevance to the study objective.
Pseudotemporal dynamics of key cell types
To delineate the pseudotemporal trajectories of astrocytes, microglia, and neurons in the TBI microenvironment, pseudotime analysis was performed. After re-clustering based on HVGs and dimensionality reduction, astrocytes resolved into six subclusters spanning 11 pseudotime states, microglia into five subclusters across 11 states, and neurons into seven subclusters over nine states. Trajectory reconstruction revealed distinct pseudotemporal paths for each cell type. Comparative analysis between TBI and control groups highlighted predicted dynamic state transitions in TBI-associated astrocytes and neurons (Supplementary Fig. 2), as well as microglia (Fig. 2).
Gene expression dynamics along pseudotime were quantified. For astrocytes, Cst3 and Slc1a2 expression peaked at intermediate pseudotime points. Microglia displayed progressive upregulation of Jun, Fos, and P2ry12. Neurons exhibited declining Calm2 expression.
Intercellular communication networks
CellChat-based network quantification of the merged scRNA-seq dataset suggested that endothelial cells had relatively high inferred communication centrality, with prominent predicted interactions with astrocytes. Within the immune-related communication network, microglia also showed substantial inferred signaling connectivity with other cell types. Predicted bidirectional communication patterns were identified among oligodendrocytes, astrocytes, microglia, neurons, endothelial cells, oligodendrocyte precursor cells, and mural cells, indicating that each cell type could act as both a predicted signal sender and receiver. Representative ligand–receptor pairs meeting the predefined CellChat criteria included Tnf–Tnfrsf1a, Sema6a–Plxna2, and Hbegf–Egfr. In particular, the Tnf–Tnfrsf1a pair was identified as a predicted microglia-associated interaction with astrocytes and neurons. These results represent computationally inferred communication patterns and do not demonstrate actual ligand secretion, receptor activation, downstream inflammation or apoptosis, or increased signaling after TBI. Detailed inferred cellular interaction networks and ligand–receptor interaction weights are shown in Supplementary Fig. 3.
Identification of microglia-associated gene modules via hdWGCNA
To delineate co-expression networks underlying microglial responses in TBI, hdWGCNA was performed on scRNA-seq data using a soft-thresholding power of β = 6 (which yielded a scale-free topology fit index of approximately 0.8) (Supplementary Fig. 4). Subsequently, nine distinct co-expression modules were identified. Spearman correlation analysis revealed correlations among distinct modules (P < 0.05). Furthermore, 90 MRGs were identified across the nine modules (Supplementary Table. 2).
Identification of candidate genes in TBI
To further investigate the transcriptional alterations in TBI, differential expression analysis was conducted between TBI samples and control samples in the training set. Bulk exploratory analysis yielded a total of 2,195 nominally differentially expressed transcripts following trauma (based on an exploratory threshold of raw P < 0.05 and |log2FC| > 0.5). However, after applying BH multiple testing correction, 89 genes passed the rigorous adjusted P < 0.05 threshold (Fig. 3a and Supplementary Table. 3). Additionally, a heatmap showing the expression patterns of the top 20 differentially expressed genes (10 upregulated and 10 downregulated) between TBI and control samples was generated (Supplementary Fig. 5). By mapping the global injury-induced differential signatures against our isolated microglial modules, we pinpointed 43 core transcripts occupying this structural intersection (Fig. 3c). Because the MRGs were derived from microglial hdWGCNA modules, this intersection prioritized genes associated with microglial co-expression patterns. However, bulk-tissue differential expression cannot establish the cell-type-specific origin of these changes. These candidate genes could provide insight into potential therapeutic approaches for TBI.
Key genes involved in a variety of biological functions in TBI
Functional annotation of candidate genes revealed significant enrichment in BP and signaling pathways associated with TBI pathophysiology. GO analysis identified 218 enriched terms (adjusted P < 0.05), including 202 BP, 14 MF, and 2 CC terms (Supplementary Table. 4). Key BP terms were linked to immune regulation and cellular stress responses, such as “leukocyte migration”, “mononuclear cell migration”, and “extrinsic apoptotic signaling pathway”. MF terms predominantly included “extracellular matrix binding” and “chemokine activity”. CC terms highlighted subcellular localization to the Golgi apparatus subcompartment and the centriole.
KEGG pathway analysis identified nine enriched pathways. Cytokine–cytokine receptor interaction and antigen processing and presentation were among the prominently enriched pathways, together with pathways related to prion disease, efferocytosis, infectious disease responses, and lysosomal function (Supplementary Fig. 6). These enrichment results indicate pathway-level associations and do not establish the corresponding cellular or molecular mechanisms.
Identification of key genes in TBI
To elucidate the functional interplay of candidate genes at the protein level, a PPI network was constructed using the STRING database (confidence score ≥ 0.4). Among the 43 intersection genes, D17H6S56E-5 was excluded due to a lack of interaction data in the STRING database. Among the remaining 42 candidate genes, 10 genes were identified as isolated nodes without structural edges. The final retained functional network comprised 32 active nodes and 79 interacting edges. Topological analysis revealed 10 genes (Ctss, Mki67, Ccl3, Birc5, Lgals3, Il1a, Ccl4, Egr1, Rrm2, and Nusap1) as hub genes with the highest degree centrality.
Cross-dataset validation identified Ccl4, Lgals3, and Egr1 as candidate genes showing consistent directional expression changes. Ccl4 and Lgals3 were upregulated, whereas Egr1 was downregulated in the validation dataset, and these expression directions were further supported by RT-qPCR (Fig. 3d–f).
To explore their potential temporal relationships, the expression patterns of Ccl4, Lgals3, and Egr1 were retrospectively mapped onto the inferred microglial pseudotime trajectory. Egr1 expression decreased during earlier portions of the inferred trajectory, whereas Ccl4 and Lgals3 increased during later portions (Supplementary Fig. 2). Because pseudotime represents a computational ordering of cross-sectional data, these patterns do not establish chronological precedence or causal regulation. They support a hypothesis-generating association among the three genes that requires direct experimental validation.
Global biological pathway alterations in TBI
GSEA was employed to identify coordinated global pathway alterations in TBI. Preranked GSEA based on the limma moderated t-statistic (TBI versus control) identified significantly enriched pathways (|NES| > 1, FDR < 0.25, nominal P < 0.05). The top five significantly upregulated and downregulated pathways in TBI are presented in Figure 4a and Figure 4b. The complete list of all significantly enriched GSEA pathways can be found in Supplementary Table. 5. GSVA further delineated hallmark pathway activity differences between TBI and control groups. At the group comparison level, angiogenesis was markedly upregulated (adjusted P < 0.05), while heme metabolism and oxidative phosphorylation pathways were suppressed (adjusted P < 0.05).
Immune infiltration analysis in TBI
Deconvolution of the GSE58485 samples suggested descriptive differences in leukocyte signatures between the TBI and control groups (Fig. 5a). Given the limited sample size of three mice per group, these estimates were not subjected to confirmatory hypothesis testing. The estimated M0 macrophage signature appeared higher in the TBI samples and was therefore reported as an exploratory descriptive trend. Among the pairwise correlations between immune cell signatures, mast cells and memory B cells showed the highest correlation coefficient.
Exploratory Spearman correlation analysis was performed to assess associations between key gene expression and immune cell signatures. Lgals3 showed a strong positive descriptive correlation with the M0 macrophage signature (r = 0.85), whereas Egr1 showed a strong negative descriptive correlation (r = −0.76). Ccl4 showed a moderate positive descriptive correlation with the same signature (r = 0.50) (Fig. 5b). Because these estimates were based on only six samples, they are highly uncertain and should be interpreted solely as exploratory effect-size estimates rather than confirmatory statistical evidence. In addition, the ImmuCC-derived M0 macrophage signature is a computationally inferred state and does not establish the presence, abundance, or influx of a discrete in vivo M0 macrophage population. It may instead reflect unpolarized or transitioning monocyte/macrophage-like transcriptional patterns. Overall, these findings suggest candidate associations between Ccl4, Lgals3, Egr1, and immune cell signatures after TBI, but do not demonstrate immune cell recruitment or direct immune regulation.
Prediction of TFBS
To assess the structural plausibility of EGR1-mediated regulation, we performed in silico motif scanning of the Ccl4 and Lgals3 promoter regions. FIMO identified four predicted EGR1 motif occurrences in the Ccl4 promoter and six predicted occurrences at three unique locations in the Lgals3 promoter, all meeting the prespecified threshold of P < 0.0001 (Fig. 6a and Supplementary Table. 6). These computational findings support a candidate regulatory hypothesis but do not demonstrate EGR1 binding, target-gene status, or transcriptional repression. Direct validation using approaches such as ChIP-qPCR or ChIP-seq and Egr1 perturbation is required.
Drug prediction of key genes
Finally, the DGIdb search yielded database-derived candidate gene–drug associations requiring experimental validation. Specifically, candidate associations were identified for Ccl4 (cyclosporine, clodronate disodium, and epoetin alfa), Lgals3 (e.g., olitigaltin, davanat, and belapectin), and Egr1 (brivoligide and genipin). Detailed parameters for each candidate association, including the specific interaction types (e.g., inhibitor, antagonist) and source databases, are provided in Supplementary Table. 7.
Discussion
This study integrated single-cell and bulk transcriptomic data to identify Ccl4, Lgals3, and Egr1 as candidate microglia-associated genes. Pseudotime, motif-scanning, and expression analyses collectively support a correlative and hypothesis-generating Ccl4/Lgals3/Egr1 expression pattern associated with inflammatory and immunometabolic changes after TBI. These findings provide hypothesis-generating insights into the immunometabolic processes associated with TBI. Rather than representing an isolated anatomical disruption, brain trauma initiates secondary molecular cascades involving neuroinflammation, metabolic dysfunction, and cell death. Because these secondary processes develop over hours to days, they may provide a potential therapeutic window for limiting further tissue damage.3,4,14 Previous studies have emphasized that microglial transcriptional states are context-dependent and functionally diverse.12,34 By combining predictive trajectory modeling with metacell-based modular analysis, our study highlights the candidate Ccl4/Lgals3/Egr1 expression pattern as a potential contributor to microglia-associated neuroinflammatory processes. These expression changes may serve as candidate biomarkers and may be associated with remodeling of the lesion microenvironment, although their regulatory roles require direct experimental validation.
Functioning as a zinc-finger transcription factor, EGR1 has been implicated in neural and inflammatory responses across several experimental contexts.35-37 Altered Egr1 expression has also been observed in TBI-associated microglial transcriptional responses.38 Although elevated EGR1 activity has been associated with neuroinflammation and neuronal injury in experimental Parkinson’s disease and intracerebral hemorrhage models,39,40 its precise regulatory role in the TBI microenvironment remains context-dependent.
Within the analyzed datasets, Egr1 expression decreased while Ccl4 and Lgals3 expression increased. The predicted EGR1 motif occurrences in both promoters provide an in silico basis for a candidate regulatory relationship. Pseudotime mapping placed the decrease in Egr1 earlier along the inferred trajectory than the increases in Ccl4 and Lgals3; however, pseudotime does not establish chronological or causal ordering. These observations therefore generate the hypothesis that reduced EGR1 activity may be associated with altered Ccl4 and Lgals3 expression after TBI. Direct evidence of promoter occupancy, transcriptional repression, and downstream functional effects is lacking and requires Egr1 perturbation and ChIP-based validation. The accompanying immune and metabolic pathway changes should likewise be interpreted as associations rather than consequences of Egr1 loss.
Ccl4 was upregulated after TBI and may represent one component of the microglial inflammatory response.41 In the CIBERSORT analysis, Ccl4 showed a moderate positive descriptive correlation with the M0 macrophage signature (r = 0.50). Given the small training cohort (n = 6), this association is exploratory and does not establish that Ccl4 drives macrophage recruitment. Nevertheless, it is consistent with the established chemotactic functions of CCL4 and supports further investigation of its relationship with infiltrating or transitioning monocyte-derived macrophages after TBI.12,41 The descriptive correlations of Lgals3 and Egr1 with the same signature likewise support candidate associations rather than a confirmed coordinated mechanism.
Beyond its canonical proinflammatory effects, CCL4 has been linked to oxidative stress and endothelial inflammation in aging-related vascular dysfunction.42 Mechanistically, CCL4 and related ligands secreted by activated microglia engage the neuronal CCR5 receptor to initiate downstream cascades, triggering mTORC1 pathway activation and impairing autophagic clearance of aggregation-prone proteins.43,44 CCL4 and related ligands have been implicated in CCR5-associated inflammatory and autophagy pathways in other neurological and inflammatory conditions.42-46 These previous findings provide biological context for the present associations but do not establish an equivalent mechanism in acute TBI. The DGIdb query identified cyclosporine and other candidate compounds associated with Ccl4. Cyclosporine has previously been investigated in TBI studies.47,48 However, the predicted Ccl4–drug associations identified here require direct experimental validation. The precise temporal and context-dependent roles of CCL4 after TBI remain to be definitively determined in future longitudinal studies.
Concurrently, this candidate expression pattern involves the induction of Lgals3, which has been implicated as a potential modulator of neuroinflammatory responses in various neurological conditions. Encoding galectin-3, a β-galactoside-binding lectin implicated in TLR4- and TREM2-mediated pattern recognition, Lgals3 serves as both a phenotypic marker of disease-associated microglia in neuroinflammatory disease and a biomarker of central nervous system inflammation.49-53
Our analyses associated Lgals3 expression with lysosomal and ribosomal pathways and identified a descriptive correlation between Lgals3 and the M0 macrophage signature. These associations do not establish that Lgals3 fixes microglia in a deleterious state, promotes leukocyte trafficking, or drives downstream cellular injury. Previous studies have implicated galectin-3 in inflammatory signaling and microglial responses in several neurological conditions,49-56 providing biological context for further investigation. According to the DGIdb database,57 the compounds associated with Lgals3 should be regarded as database-generated candidates requiring direct pharmacological and experimental validation.
CellChat analysis identified predicted microglial ligand–receptor patterns, including the Tnf–Tnfrsf1a and Hbegf–Egfr axes. These computationally inferred interactions are biologically plausible in light of previous evidence that activated microglia may influence astrocytic and neuronal responses after TBI.7,12,50,52,58 However, because the CellChat analysis was performed on merged scRNA-seq data and was not experimentally validated, it does not establish enhanced intercellular signaling, astrocyte polarization, neuronal injury, or a causal role of the Ccl4/Lgals3/Egr1 expression pattern. These findings therefore provide a hypothesis-generating framework for future investigation of microglial–neural crosstalk after TBI.11,34
Limitations
Several limitations should be considered. First, the integration of in silico motif scanning, pseudotime analysis, CellChat inference, bulk transcriptomics, and RT-qPCR provides correlative and hypothesis-generating evidence only. It does not establish direct EGR1 binding, transcriptional regulation of Ccl4 or Lgals3, microglial-state transitions, or downstream immune and metabolic effects. These relationships require direct validation using Egr1 perturbation models, ChIP-qPCR or ChIP-seq, and protein-level assays such as Western blotting or immunofluorescence. Second, the analyses integrated murine datasets generated using different injury models, brain regions, post-injury time points, sexes, and transcriptomic platforms. These differences may limit cross-dataset comparability and direct translation to human TBI. In addition, the contralateral tissues used as comparators in GSE242025 may have been influenced by systemic neuroinflammatory crosstalk or diaschisis and may not represent completely uninjured brain tissue. Third, although the candidate genes were initially identified from microglial single-cell modules, downstream validation relied on bulk transcriptomic data and whole-tissue RT-qPCR. Contributions from other cell types, including reactive astrocytes and endothelial cells, therefore cannot be excluded. Fourth, the small GSE58485 training cohort (n = 6) substantially limited the precision and statistical power of the immune-cell proportion and gene–cell correlation analyses. These analyses were therefore interpreted as exploratory and descriptive rather than confirmatory. In addition, the ImmuCC-derived “M0 macrophage” signature is a computational construct that may not reliably distinguish resident microglia from infiltrating monocyte-derived macrophages in murine brain tissue. Fifth, although Ccl4, Lgals3, and Egr1 showed nominal differential expression, they did not meet the FDR < 0.05 threshold in the small training dataset, with adjusted P-values of approximately 0.147, 0.075, and 0.132, respectively. Their selection therefore depended on consistency across datasets and RT-qPCR validation and should be interpreted cautiously. Finally, the RT-qPCR experiments included only male mice, preventing assessment of sex-related differences in post-TBI immune responses. The available datasets also did not permit full resolution of the spatial and temporal dynamics of the candidate Ccl4/Lgals3/Egr1 expression pattern. Future studies using standardized sampling, balanced-sex cohorts, spatial transcriptomics, conditional genetic models, and direct mechanistic assays are required to clarify the spatial and temporal dynamics and causal relevance of this candidate Ccl4/Lgals3/Egr1 expression pattern.
Supporting information
Supplementary material for this article is available at https://doi.org/10.14218/NSSS.2026.00014 .
Supplementary Table. 1
Summary of detailed statistical considerations across all bioinformatic and experimental analyses
(XLSX)
Supplementary Fig. 1
Single-cell analysis of the cell-type proportions after TBI. (a) The features, counts, and percentages of mitochondrial genes in each of the analyzed samples before quality control. (b) The features, counts, and percentages after quality control. (c) Highly variable genes were identified, with the top 10 genes exhibiting the greatest variation. (d) The elbow plot inflection shows the standard deviation of each PC; the standard deviations fluctuate markedly before PC 17. The first 17 PCs were retained for cell clustering. (e) Principal component analysis of the first 17 PCs using JackStrawPlot. (f) The UMAP plot visualizes the distribution of 14 clusters of cells. mt, mitochondrial; PC, principal component; RNA, ribonucleic acid; TBI, traumatic brain injury; UMAP, Uniform Manifold Approximation and Projection.
(TIF)
Supplementary Fig. 2
Pseudotime trajectory analyses of the three major neural cell types. (a–c) Pseudotime analysis and pseudotemporal trajectories of neurons (a), astrocytes (b), and microglia (c). (d–f) Comparative cell distribution along the pseudotime trajectory between the TBI and control groups for neurons (d), astrocytes (e), and microglia (f). (g, h) Dynamic expression patterns of the top variable genes along the pseudotime trajectory for astrocytes (g) and neurons (h). (i) Targeted retrospective mapping of Ccl4, Lgals3, and Egr1 along the inferred microglial pseudotime trajectory. Egr1 expression decreased during earlier portions of the inferred trajectory, whereas Ccl4 and Lgals3 increased during later portions. This computational ordering does not establish chronological precedence or causal regulation. TBI, traumatic brain injury.
(TIF)
Supplementary Fig. 3
The cellular interaction network among seven types of cells. (a) Number of interactions. (b) Interaction weights/strength. (c) Bubble plot of interactions between different cell types.
(TIF)
Supplementary Fig. 4
Identification of module genes in microglial cells. (a) The scale-free exponent and average connectivity of various soft thresholds. (b) Dendrogram of gene clustering. (c) Correlation between modules. MG, microglia.
(TIF)
Supplementary Table. 2
Detailed module assignments and kME values for the 90 MRGs identified via hdWGCNA
(XLSX)
Supplementary Table. 3
Differentially expressed genes between TBI and Control samples
(XLSX)
Supplementary Fig. 5
Heatmap of the top differentially expressed genes. Expression profiles of the top 10 significantly upregulated and top 10 downregulated nominally differentially expressed genes between the TBI and Control groups in the GSE58485 dataset. The heatmap was generated using the ComplexHeatmap R package. Rows represent individual genes, and columns represent distinct samples. Color scaling indicates relative expression levels. TBI, traumatic brain injury.
(TIF)
Supplementary Table. 4
Key genes involved in a variety of biological functions in TBI based on GO
(XLSX)
Supplementary Fig. 6
Functional annotation of candidate genes. (a) Bubble plot indicating the enrichment pathways based on GO. (b) Enrichment pathways based on KEGG. GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes. BP, biological process; CC, cellular component; CCR, C-C chemokine receptor; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; MF, molecular function.
(TIF)
Supplementary Table. 5
Significantly enriched GSEA pathways in TBI vs. Control
(XLSX)
Supplementary Table. 6
FIMO prediction of EGR1 binding sites on the promoters of target genes
(XLSX)
Supplementary Table. 7
Database-derived candidate gene-drug associations identified via DGIdb
(XLSX)
Declarations
Funding
This work was supported by the Key Project of the Qinghai Provincial Health System (grant number 2025-wjzd-11) and the Special Project for Scientific and Technological Research of the Sichuan Provincial Administration of Traditional Chinese Medicine (grant number 25MSZX474).
Conflict of interest
The authors declare that they have no conflicts of interest.
Author contributions
Conceptualization (TY, SY), software (TY, DZ, KF, JC), formal analysis (TY, DZ, KF, JC, YY), writing - original draft (TY, YY), writing - review & editing (TY, YY), supervision (SY), project administration (SY), funding acquisition (SY), validation (DZ, KF, JC), investigation (DZ), data curation (DZ), methodology (YM), and visualization (KF, JC). All authors have approved the final version and publication of the manuscript.
Ethical statement
All animal experiments were approved by the Ethics Committee of The General Hospital of Western Theater Command (Approval No. 2023ky119-1) and were conducted in strict accordance with the NIH Guide for the Care and Use of Laboratory Animals and were reported in accordance with the ARRIVE 2.0 guidelines. All efforts were made to minimize the number of animals used and their suffering.
Data sharing statement
The transcriptomic datasets analyzed in this study are publicly available and can be found in the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). The specific datasets utilized include: Training Set: GSE58485 (GPL1261), Validation Set: GSE242025 (GPL20258), Single-Cell Dataset: GSE101901 (GPL17021). Furthermore, regulatory and functional data supporting this research were obtained from the following repositories: Transcription Factor Binding Sites: The JASPAR database (EGR1: MA0162.1) and the mouse reference genome assembly mm10 (GRCm38), Cell Markers: The CellMarker database (http://117.50.127.228/CellMarker/), Drug-Gene Interactions: The Drug Gene Interaction Database (DGIdb) (https://www.dgidb.org). The raw and processed data generated during the current study, including the results of the hdWGCNA, pseudotime analysis, and RT-qPCR experimental validation, are not publicly archived but are available from the corresponding author upon reasonable request.