Identification and Mechanistic Studies of Key Genes in Thalamic Hemorrhage Pain by Multi-omics

Abstract

Background:

Thalamic hemorrhage pain (THP), a subtype of central post-stroke pain (CPSP), commonly develops following ischemic or hemorrhagic injury to the thalamus. Current therapeutic options remain inadequate due to the absence of well-defined molecular targets. This study aimed to elucidate critical genes implicated in THP pathogenesis through an integrated multi-omics approach.

Methods:

A mouse model of THP was established and mice were divided into THP and control groups. Comprehensive multi-omics profiling involving transcriptomics, proteomics, metabolomics, ribosome profiling (Ribo-seq), and single-cell RNA sequencing (scRNA-seq) was conducted. Differentially expressed genes (DEGs), differentially expressed proteins (DEPs), ribosome footprint-associated DEGs (RF-DEGs), and differentially expressed metabolites (DEMs) were identified via comparative expression analyses. Hub genes were extracted from the DEGs and subsequently intersected with scRNA-seq DEGs, DEPs, and RF-DEGs to define key gene candidates. These genes underwent gene set enrichment analysis (GSEA), disease association mapping, and drug prediction. Expression levels of key genes were used to delineate critical cell populations, followed by analyses of intercellular communication and pseudotemporal differentiation trajectories. Orthogonal partial least squares discriminant analysis was used to validate the model.

Results:

The THP mouse model was successfully validated. Multi-omics analyses yielded distinct profiles of DEGs, single-cell DEGs, DEPs, RF-DEGs, and DEMs, which were functionally annotated through enrichment strategies. Notably, 12 hub genes were prioritized, of which eight key genes (ferritin light chain 1 (Ftl1), tropomyosin 4 (Tpm4), C-C motif chemokine ligand 3 (Ccl3), C-C motif chemokine ligand 4 (Ccl4), C-C motif chemokine receptor 2 (Ccr2), interleukin 33 (Il33), C-X-C motif chemokine ligand 2 (Cxcl2), and Lymphocyte antigen 6 complex, locus C2 (Ly6c2) were identified. These genes were predominantly associated with oxidative phosphorylation and ribosomal pathways. Further analyses revealed strong associations with necrotic and inflammatory processes, and compounds such as alprostadil and anisomycin were identified as potential therapeutic agents. Single-cell analyses highlighted six pivotal cell types, including endothelial cells and macrophages. Intercellular communication networks and lineage progression patterns of these cells were systematically characterized, alongside spatial and temporal expression profiles of key genes.

Conclusions:

This study established a validated THP mouse model and employed a multi-omics integration strategy to identify eight key genes and associated molecular pathways. These findings provide novel mechanistic insights into THP pathogenesis and highlight promising targets for therapeutic intervention.

Graphical abstract

Keywords

thalamic hemorrhage / pain / genomics / proteomics / metabolomics / RNA sequencing / genes

1. Introduction

Stroke represents the primary cause of disability and mortality in the aging population [1], profoundly impairing quality of life and imposing significant financial burden on caregivers. Central post-stroke pain (CPSP), classified as a central neuropathic pain syndrome, arises from lesions in the somatosensory system following stroke. The prevalence of CPSP following hemorrhagic stroke has been reported to range between 8% and 46%. The thalamus, a vital relay center for nociceptive transmission, is frequently implicated in the pathogenesis of CPSP, with thalamic hemorrhage identified as a major etiological factor. When hemorrhagic or ischemic events involve the thalamus, the incidence of CPSP may escalate to 80% [2, 3, 4]. Lesions associated with CPSP frequently localize to the ventral posterolateral (VPL) and ventral posteromedial (VPM) nuclei of the thalamus—regions notably susceptible to CPSP development [5]. Current pharmacological interventions are largely inadequate, often yielding limited therapeutic benefit and frequent adverse effects [6, 7, 8]. Despite the clinical burden, the mechanistic role of CPSP remains insufficiently elucidated, thereby constraining the advancement of effective and targeted therapeutic strategies.

Recent preclinical investigations into the molecular pathophysiology of thalamic hemorrhage (TH) have revealed novel biomarkers with potential prognostic and therapeutic relevance [9, 10, 11]. Cerebral hemorrhage has shown to induce widespread transcriptional alterations, highlighting dysregulated genes as viable targets for pharmacological intervention in CPSP [12, 13]. Technological advancements in high-throughput sequencing have enabled the generation of diverse omics datasets, including transcriptomics, metabolomics, ribosome profiling, and proteomics, facilitating multidimensional understanding of disease processes that surpasses the limitations of single-layer analyses [14, 15]. Integrative multi-omics approaches allow for the identification of biologically meaningful biomarkers by revealing cross-talk among various molecular pathways. In parallel, single-cell RNA sequencing (scRNA-seq) significantly enhances data resolution by enabling gene expression profiling at the individual cell level. This technique captures cellular heterogeneity and delineates dynamic transitions in cellular states and transcriptional programs during disease progression [16, 17, 18]. Despite these advances, current CPSP research has predominantly concentrated on isolated molecules or discrete signaling pathways, lacking a comprehensive systems-level perspective. A biologically integrated understanding of CPSP pathogenesis necessitates holistic analyses that unify data across multiple molecular layers.

To address this gap, the present study established a thalamic hemorrhage pain (THP) mouse model and conducted comprehensive multi-omics analyses, including transcriptomics, proteomics, metabolomics, ribosome profiling, and single-cell transcriptomics, on thalamic tissue. The resulting datasets were systematically processed and subjected to bioinformatics tools to identify critical genes. By integrating transcriptomic, single-cell, and proteomic data, key molecular regulators were identified, providing mechanistic insights into THP pathogenesis and highlighting potential diagnostic and therapeutic targets. This integrative approach enabled the construction of a dynamic regulatory network governing THP, providing a robust framework for future therapeutic discovery.

2. Materials and Methods

2.1 Animal Modeling and Sample Collection

A total of 120 adult C57BL/6 mice (7–8 weeks old, 25–30 g) were acquired from Beijing Weitong Lihua Experimental Animal Technical Co., Ltd. (Beijing, China). The TH model was established by Wuhan Servicebio Biotechnology Co., Ltd. (Wuhan, Hubei, China) and approved by Beijing Weitong Lihua Experimental Animal Technical Co., Ltd. (Approval No. 2022045). All subsequent procedures, including behavioral assessments, tissue collection, and downstream analyses, were conducted under ethical approval of the Animal Care and Use Committee of the Medical College of Yangzhou University (Approval No. 202303874). All animal experiments were performed in strict accordance with the National Institutes of Health Guide for the Care and Use of Laboratory Animals. Mice were randomly assigned to two experimental groups: a model group (designated as the THP group) and a control group. Validation of the THP mouse model was confirmed by both hematoxylin and eosin staining and behavioral assessments as previously described [19]. Thalamic tissues were collected from both groups for comprehensive multi-omics sequencing. Sample allocation for each sequencing modality was implemented as follows: transcriptomics (n = 3 mice per group), proteomics (n = 3 mice per group), metabolomics (n = 9 mice per group), single-cell transcriptomics (n = 2 mice per group), and ribosome profiling (n = 3 mice per group). Mice were anesthetized with vaporized isoflurane (RWD Life Science, Shenzhen, Guangdong, China; Cat. No. R510-22-10), administered at 2.5% for induction and maintained at 2.0% during the procedure. Animals were then positioned in a stereotaxic apparatus (RWD, 68025, Shenzhen, Guangdong, China). Under stereotactic guidance, Collagenase IV (Coll IV; #C2139-1G, Sigma-Aldrich, St. Louis, MO, USA; 0.01 U/10 nL in sterile saline) was microinjected into the right ventral posteromedial (VPM) and ventral posterolateral (VPL) thalamic nuclei. Stereotaxic coordinates relative to bregma were: anterior-posterior (AP): –0.82 mm to –2.30 mm; medial-lateral (ML): +1.30 mm to –1.95 mm; and dorsal-ventral (DV): –3.01 mm to –4.25 mm from the skull surface, based on the Paxinos and Franklin mouse brain atlas. In the sham-operated group, an equal volume of sterile saline was injected at the same coordinates. Following injection, the glass micropipette was held in place for 10 min to ensure adequate dispersion of the Coll IV solution before being slowly withdrawn. The surgical site was then irrigated with sterile saline and disinfected with iodophor, followed by closure with wound clips.

Following decapitation, the scalp was incised to expose the skull, and all dissection procedures were performed on ice to preserve tissue integrity. The skull was carefully opened with fine scissors (Lige Science, Guangzhou, Guangdong, China; Cat No. LG01-107-5), the brain was excised, and the thalamus was isolated. Tissue samples were promptly transferred into pre-labeled EP tubes (Thermo Fisher Scientific, Waltham, MA, USA; Cat No. 508-GRD-Q), snap-frozen in liquid nitrogen (Huate Gas, Foshan, Guangdong, China; Cat No. 7727-37-9), and subsequently stored at –80 °C for long-term preservation. High-throughput sequencing and omics analyses were outsourced to certified third-party providers. Single-nucleus RNA sequencing (snRNA-seq), quantitative proteomics, and untargeted metabolomics were conducted by Hangzhou Lianchuan Biotechnology Co., Ltd. (Hangzhou, Zhejiang, China), while bulk RNA sequencing (RNA-seq) and ribosome profiling (Ribo-seq) were performed by Guangzhou Epigenetic Biotechnology Co., Ltd. (Guangzhou, Guangdong, China). All raw and processed data have been deposited in publicly accessible repositories. Bulk RNA-seq data are available in the NCBI Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) under accession number of GSE275389; proteomic data have been deposited in the ProteomeXchange Consortium via the PRIDE (https://www.ebi.ac.uk/pride) partner repository under identifier of PXD058910; metabolomic data are available in the MetaboLights database (https://www.ebi.ac.uk/metabolights/) under accession number of MTBLS11941; and snRNA-seq data have been deposited in GEO under accession number of GSE227033.

2.2 Single-cell Data Analysis

Single-cell transcriptomic analysis was performed using the Seurat v5.0.1 package (https://satijalab.org/seurat) [20], following stringent quality control (QC) procedures. Key cellular metrics, including the number of detected genes per cell (nFeature_RNA), total RNA count per cell (nCount_RNA), and mitochondrial gene expression percentage (percent.mt), were computed and visualized. Genes expressed in fewer than three cells and cells with fewer than 200 detected genes were excluded from downstream analysis. Cells meeting the criteria of nFeature_RNA \le95th percentile, nCount_RNA between the 5th and 95th percentiles, and percent.mt 5% were retained for downstream analysis. Data normalization was executed using the “NormalizeData” function, and the top 2000 highly variable genes were identified via the “FindVariableFeatures” function and exported for subsequent analyses. Dimensionality reduction was carried out using principal component analysis (PCA) through the “RunPCA” function, and significant principal components (PCs) were determined based on “JackStrawPlot” and “ElbowPlot” outputs (p 0.05). The selected PCs were utilized to perform unsupervised clustering using the FindNeighbors and FindClusters functions. Cluster resolution was further refined using t-distributed stochastic neighbor embedding (t-SNE), which was implemented via the RunTSNE function, with a resolution parameter set to 0.1 to delineate distinct cellular subpopulations.

2.3 Identification of Single-cell DEGs

Cell clusters identified via t-SNE were annotated into specific cell types using canonical marker genes obtained from the Tabula Muris Senis database (https://tabula-muris.ds.czbiohub.org/). The relative abundance of each cell type in the THP and control groups was quantified and visualized, and cell types exhibiting statistically significant differences (p 0.05) were designated as differentially represented. Differential gene expression analysis was performed using the “FindAllMarkers” function in Seurat v5.0.1 to define marker genes uniquely associated with each distinct cell subpopulation. The functional parameters were configured as follows: only.pos = TRUE to retain only upregulated markers in each subpopulation, min.pct = 0.25 to ensure markers were expressed in at least 25% of the cells, and logfc.threshold = 0.25 to include only genes with a log-transformed fold change \ge0.25. Marker gene expression profiles across differential cell types were further analyzed and visualized using scRNAtoolVis v0.1.0 (https://github.com/junjunlab/scRNAtoolVis), employing the Wilcoxon rank-sum test and specified filtering thresholds. Differential cell types were selected, and marker genes were identified using the FindAllMarkers function (parameters: only.pos = TRUE, min.pct = 0.1, logfc.threshold = 0). The direction of gene expression was then corrected to accurately reflect the subgroup-specificity. Subsequently, the gene symbols (SYMBOL) were converted to ENTREZ IDs using the bitr function. Gene Ontology (GO) biological process (BP) enrichment analysis was performed using the compareCluster function from the clusterProfiler v4.10.0 package (https://bioc.r-universe.dev/clusterProfiler). The top 5 significant pathways for each subgroup were selected based on p-values, and the ENTREZ IDs from the enrichment results were reverse-mapped to gene symbols. The enrichment results were visualized using the ggforce v0.4.2 package (https://CRAN.R-project.org/package=ggforce), and heatmaps and volcano plots were generated to display the expression patterns of the marker genes, providing a comprehensive analysis of the functional characteristics of each cell subgroup. To identify differentially expressed genes (DEGs) in specific cell types (cell-DEGs), the “FindMarkers” function in Seurat v5.0.1 was employed, with significance thresholds of p 0.05 and |log₂fold change (FC)| >1. The five most upregulated and five most downregulated cell-DEGs, which were ranked by absolute log₂FC, were labeled and exported for downstream analyses. Subsequent GO and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed to elucidate the biological functions and signaling pathways associated with these cell-DEGs, adhering to the methodology described in Section 2.1 (p 0.05). The top five enriched GO terms from each category, involving BP, cellular component (CC), and molecular function (MF), as well as the most significant KEGG pathways, were prioritized based on the number of enriched genes.

2.4 Transcriptome Analysis and Identification of Hub Genes

Transcriptomic count data were quantified using FeatureCounts v2.0.3 (https://subread.sourceforge.net/featureCounts.html) [21] and normalized to fragments per kilobase of transcript per million mapped reads (FPKM), enabling assessment of gene expression distributions in individual samples and facilitating cross-sample expression comparisons. PCA was subsequently applied for dimensionality reduction, and sample distribution across PCs was visualized via a PCA plot. Differential expression analysis between the THP and control groups was conducted using DESeq2 v1.36.0 (https://bioconductor.org/packages/release/bioc/html/DESeq2.html) [22], with significance thresholds set at p 0.05 and |log₂FC| >1. The resulting DEGs were visualized using a volcano plot generated by the ggplot2 v3.4.4 package (https://ggplot2.tidyverse.org/) [23], and the top 10 upregulated and top 10 downregulated DEGs, which were ranked by absolute log₂FC, were highlighted. Expression patterns of these DEGs were further visualized in a heatmap using the Pheatmap v1.0.12 package (https://cran.r-project.org/web/packages/pheatmap/index.html) [24]. DEGs identified through transcriptome profiling were collectively termed THP-DEGs. To elucidate the functional roles and signaling pathways associated with THP-DEGs (p 0.05 and |log₂FC| >1), GO and KEGG pathway enrichment analyses were performed using clusterProfiler v4.10.0 [23] and org.Mm.eg.db v3.18.0 packages (https://bioconductor.org/packages/release/data/annotation/html/org.Mm.eg.db.html), with p 0.05 and |log₂FC| >1 as significance criteria. GO enrichment encompassed three categories—BP, CC, and MF—and the top five terms in each category, along with enriched KEGG pathways, were ranked by gene count and visualized via ggplot2. Protein-protein interaction (PPI) analysis of THP-DEGs was conducted by uploading gene lists to the Search Tool for the Retrieval of Interaction Gene/Proteins (STRING) database (http://string-db.org/) using an interaction score threshold >0.4, excluding unconnected nodes. The resulting PPI network was visualized using Cytoscape v3.10.2 (https://cytoscape.org/) [25]. Key functional modules in the network were identified using the MCODE plug-in, while hub genes were further prioritized using five algorithms, involving Maximal Clique Centrality (MCC), Density of Maximum Neighborhood Component (DMNC), Maximum Neighborhood Component (MNC), Edge Percolated Component (EPC), and Degree, through the CytoHubba plug-in in Cytoscape. The top 30 genes identified by each algorithm were intersected to define the final set of hub genes. Their interactions were subsequently reconstructed to generate a dedicated PPI network representing core regulatory components.

2.5 Proteomic Analysis

Proteomic data were normalized using Spectronaut 15.0 software (https://biognosys.com/software/spectronaut/) [26], enabling stringent QC through spectral matching and scoring of detected peptides against theoretical spectra. Score-based probability assessments ensured accuracy, and peptides meeting the criteria of high confidence scores, false discovery rate (FDR) 0.01, and p 0.05 were retained as valid identifications. Protein annotation was conducted using the UniProt database (http://www.uniprot.org/). To evaluate sample distribution and inter-sample variability, PCA was performed with the “psych” v2.4.6 package (https://cran.r-project.org/web/packages/psych/index.html) [27], while density plots and boxplots, generated via the “ggplot2” 3.4.4 package, were employed to visualize sample dispersion, protein expression levels, and distributional characteristics. Differentially expressed proteins (DEPs) between the THP and control groups were identified using the “DEqMS” package (https://www.bioconductor.org/packages/release/bioc/html/DEqMS.html). Transcription is governed by a complex regulatory network involving transcription factors and chromatin structure, with transcript-level fluctuations subjected to amplification or attenuation during translation. A stringent threshold of |log₂FC| >1 and p 0.05 was applied to enhance the reliability of identifying transcriptional alterations and to minimize spurious findings. For proteomic analysis, a cutoff of |log₂FC| >0.26 was adopted, reflecting the typically subtler nature of protein-level changes due to regulatory influences at the post-transcriptional, translational, and degradation stages. Even modest shifts in protein expression may yield significant biological effects; thus, applying a |log₂FC| >0.25 threshold facilitates broader detection of DEPs with potential functional relevance. Volcano plots and heatmaps were generated using the ggplot2 3.4.4 and “Pheatmap” 1.0.12 packages, respectively, to illustrate the number and expression profiles of THP-DEPs. Cluster resolution in Seurat was defined via the FindClusters function (scRNA.norm.pca.clu, resolution = 0.1). Differential expression analysis in DESeq2 employed a generalized linear model (GLM). For visualization, volcano plots were constructed with ggplot. Functional annotation and pathway enrichment analyses of THP-DEPs were performed using GO and KEGG methodologies, as outlined in Section 2.1, with significance determined at p 0.05. The top 20 GO terms in each subdomain, including BP, CC, and MF, along with KEGG pathways, were ranked by ascending p-values and visualized as bubble plots created with using “ggplot2” 3.4.4. PPI networks were constructed based on the STRING database (http://string-db.org), restricted to mouse (taxid:10090), with an interaction score threshold >0.4 to retain medium-confidence interactions. During network refinement, isolated proteins and those with single-edge connectivity were excluded to improve network integrity. Interaction strength was evaluated using the combined STRING score (0–1), in which values >0.4 were considered biologically meaningful, and higher scores indicated stronger evidence. Network centrality was quantified by degree of connectivity, with nodes exhibiting higher degrees identified as central hubs. Final visualization was conducted using Cytoscape v3.10.2, where node size reflected connectivity and edge thickness denoted interaction strength.

2.6 Metabonomic Analysis

PCA for metabonomic profiling was performed using the “Gmodels” 2.19.1 package (https://github.com/r-gregmisc/gmodels). Partial least squares discriminant analysis (PLS-DA) and orthogonal PLS-DA (OPLS-DA) were implemented via the “mixOmics” 3.19 (https://mixomics.org/install/) [28] and “ropls” v1.36.0 (https://bioconductor.org/packages/release/bioc/html/ropls.html) [29] packages, with the latter serving to validate model robustness. Visualization was carried out using the “ggplot2” 3.4.4 package. To mitigate overfitting, model reliability was evaluated through seven-fold cross-validation and 200 iterations of response permutation testing (RPT). Differential metabolites (THP-DEMs) between the THP and control groups were identified through combined multidimensional and unidimensional statistical approaches, with significance defined as p 0.05 and variable importance in projection (VIP) >1. Expression levels and quantitative distribution of THP-DEMs were visualized using volcano plots and heatmaps generated via “ggplot2” 3.4.4 and “Pheatmap” 1.0.12. Functional interpretation was achieved through KEGG enrichment analysis using the “MetaboSignal” 3.19 package (https://www.bioconductor.org/packages/release/bioc/html/MetaboSignal.html) [30], applying a significance threshold of p 0.05. Enriched pathways were ranked by ascending p-values to highlight the most relevant metabolic processes.

2.7 Ribosome Profiling Sequencing (Ribo-seq)

For Ribo-seq data preprocessing, QC was applied to exclude reads containing adapter sequences, over 10% ambiguous bases (N), or a quality score below 20 (Q20) for \ge50% of the read length. Polyadenine (poly-A) reads were also discarded. Non-coding RNAs including ribosomal RNAs (rRNAs), transfer RNAs (tRNAs), small nucleolar RNAs (snoRNAs), small nuclear RNAs (snRNAs), and microRNAs (miRNAs) were filtered out to improve specificity. The remaining reads were aligned to the mouse reference genome to generate ribosome footprints (RFs), which were subsequently normalized to FPKM for downstream analyses. Differentially expressed RFs (RF-DEGs) between the THP and control groups were identified using the “DESeq2” 1.36.0 package, with statistical cutoffs set at p 0.05 and |log₂FC| >1. A volcano plot constructed with “ggplot2” 3.4.4 highlighted the top 10 upregulated and downregulated RF-DEGs based on descending |log₂FC| values. Heatmaps generated via the “Pheatmap” 1.0.12 package further depicted expression patterns and clustering across samples. Functional enrichment analysis of RF-DEGs was performed using the “clusterProfiler” 4.10.0 and “org.Mm.eg.db” 3.18.0 packages, with enriched GO terms and KEGG pathways identified at p 0.05. The top five GO terms for BP, CC, and MF, ranked by gene count, were presented alongside KEGG pathway results to elucidate the underlying functional landscape.

2.8 Identification of Key Genes and Functional Analyses of Key Genes

Key genes were identified by integrating hub gene sets derived from THP-DEPs, cell-DEGs, and RF-DEGs. Specifically, “protein_hub” and “ribosome_hub” genes were defined by the intersection of hub genes with THP-DEPs and RF-DEGs, respectively. Overlaps between these hub gene subsets and cell-DEGs were further analyzed using the “UpSetR” 1.4.0 package (https://doi.org/10.1093/bioinformatics/btx364), enabling the final selection of key genes. To elucidate the biological functions and pathways associated with these key genes, gene set enrichment analysis (GSEA) was performed utilizing the mouse KEGG gene sets from the Molecular Signatures Database (MSigDB, http://software.broadinstitute.org/gsea/msigdb/index.jsp). Spearman correlation coefficients were computed between each key gene and all transcriptome-wide genes, followed by ranking based on descending correlation values. GSEA was then conducted using the “clusterProfiler” 4.10.2 package, applying thresholds of adj.p 0.05 and |normalized enrichment score (NES)| >1. The top five enriched pathways for each key gene, ranked by descending |NES|, were reported. To further characterize the biological functions and interaction profiles of the key genes, a gene-gene interaction (GGI) network was constructed using the GeneMANIA platform (https://genemania.org/). Functional similarity among key genes was assessed through “Friends” analysis using the “GOSemSim” 2.24.0 package (https://www.bioconductor.org/packages/release/bioc/html/GOSemSim.html) [31], with statistical significance defined as p 0.05. Elevated Friends scores indicated stronger functional relatedness. Chromosomal localization of key genes was determined and visualized via the “RCircos” 1.2.2 package (https://cran.r-project.org/web/packages/RCircos/index.html) [32].

2.9 Disease Correlation Analysis and Drug Prediction

To investigate potential disease associations, key genes were analyzed using the Comparative Toxicogenomics Database (CTD, https://ctdbase.org/search/), identifying the top ten genes with the highest disease relevance. Drug-target associations were predicted using the DSigDB database (https://dsigdb.tanlab.org/DSigDBv1.0/), revealing compounds with potential therapeutic relevance for the identified key genes. The resulting disease-gene and drug-gene interaction networks were visualized using Cytoscape 3.10.2.

2.10 Identification of Key Cells, Cellular Communication and Pseudotime Analysis

To identify key cell populations associated with the expression of pivotal genes at the single-cell level, the AddModuleScore function from the Seurat package (v5.0.1) was utilized, using a predefined set of key genes. Differential expression analysis between THP and control groups was conducted in immune-related cell subsets using the Wilcoxon test (p 0.05), enabling the identification of cell types with significant alterations in gene expression. To further elucidate the roles of these key cells, cell–cell communication and pseudotime trajectory analyses were performed. Intercellular signaling networks were inferred using the CellChat package (v1.6.1) (https://github.com/jinworks/CellChat) [33], with ligand–receptor interactions mapped based on annotations from the CellChatDB database (https://deepwiki.com/jinworks/CellChat/2.2-cellchatdb/). Temporal dynamics and differentiation pathways of the identified key cells were explored via pseudotime analysis using Monocle2 (v2.26.0) (https://bioconductor.org/packages/release/bioc/html/monocle.html) [34]. The Branched Expression Analysis Modeling (BEAM) method was applied to capture gene expression trends along developmental branches, and results were visualized using the “plot pseudotime heatmap” function.

2.11 RNA Extraction and RT-qPCR

Total RNA was isolated from the thalamic tissue of mice three days post-Coll IV injection using the RNA Easy Fast Tissue Kit (DP451, TIANGEN, Beijing, China), followed by reverse transcription into cDNA with the FastQuant RT Kit containing gDNase (KR116, TIANGEN), in accordance with the manufacturer’s protocol. Primers targeting eight key gene transcripts were designed based on NCBI database sequences (https://www.ncbi.nlm.nih.gov/) using Primer3 (https://primer3.org/) (Table 1). Quantitative real-time PCR was performed using ChamQ SYBR qPCR Master Mix (Q311-02, Vazyme, Nanjing, Jiangsu, China) on the Step One Plus Real-Time PCR System (Applied Biosystems, Foster City, CA, USA). Relative transcript levels were quantified using the 2^{-Δ⁢Δ⁢CT} method [35], with ACTB serving as the internal control. Normalized gene expression was expressed as a fold change relative to the mean expression observed in the sham group.

2.12 Statistical Analysis

Bioinformatic analyses were performed using R 4.2.2 software (https://cran.r-project.org/bin/windows/base/old/). Group comparisons were performed using the Wilcoxon test, and statistical significance was defined as p 0.05.

3. Results

3.1 Selection of Differential Cells

Following QC, 28,750 cells and 32,935 genes were retained for single-cell analysis (Supplementary Fig. 1A), and the top 2000 most variable genes were identified (Supplementary Fig. 1B). Based on the PCA scree plot and fragmentation profile, the first 20 PCs were selected for downstream clustering (p 0.05) (Supplementary Fig. 1C,D). Notably, 14 distinct cell clusters were delineated using t-SNE (Fig. 1A) and annotated into 10 cell types according to canonical marker genes: endothelial cells (ECs), oligodendrocytes, astrocytes, brain pericytes, Bergmann glia cells, macrophages, oligodendrocyte precursor cells (OPCs), neurons, microglial cells, and natural killer cells (Fig. 1B). Among them, ECs, oligodendrocytes, astrocytes, brain pericytes, Bergmann glia cells, macrophages, OPCs, neurons, and natural killer cells displayed altered distributions between the THP and control groups and were designated as differential cell types (Fig. 1C). A substantial number of marker genes demonstrated significant expression changes across these cell populations, delineating distinct transcriptional signatures for each subpopulation (Supplementary Fig. 2A,B). The results of the functional enrichment analysis indicated that endothelial cells were primarily enriched in blood coagulation and fibrin clot formation, while oligodendrocyte precursor cells were primarily enriched in protein maturation (Supplementary Fig. 2C). Differential expression analysis yielded 14,643 cell-specific DEGs, comprising 5739 upregulated and 8904 downregulated genes (Fig. 1D). GO enrichment analysis highlighted biological processes, such as synapse organization and axonogenesis, cellular components including the synaptic membrane and neuron-to-neuron synapse, and molecular functions involving passive transmembrane transporter and channel activities (Fig. 1E). KEGG pathway analysis further revealed enrichment in neurodegenerative disease pathways, including those related to Alzheimer’s disease and human papillomavirus infection (Fig. 1F).

3.2 Identification of Hub Genes Using Transcriptome Analysis

Transcriptome profiling identified 12 hub genes with potential relevance to THP-associated pathology. QC metrics for the transcriptomic dataset are detailed in Supplementary Fig. 3. Differential expression analysis detected 376 THP-DEGs, including 211 upregulated and 165 downregulated transcripts (Fig. 2A). The top 20 THP-DEGs, ranked by absolute log₂FC, were further analyzed for expression patterns and distribution (Fig. 2B). GO enrichment analysis yielded 1569 significant terms (Supplementary Table 1), emphasizing biological processes, such as the positive regulation of responses to external stimuli and secretion, cellular components including secretory granules and membrane microdomains, and molecular functions, involving G protein-coupled receptor binding and passive transmembrane transporter activity (Fig. 2C). KEGG pathway analysis indicated significant enrichment in pathways related to neuroactive ligand-receptor interactions, cytokine-cytokine receptor interactions, calcium signaling, and additional intercellular signaling cascades (Fig. 2D), underscoring their potential role in cellular communication and signal mediation. A PPI network was constructed after removing isolated and singly connected nodes, revealing Il6 as the central node with the highest connectivity (Fig. 2E). Integration of five distinct analytical algorithms led to the consistent identification of 12 hub genes: Ccl3, Ccl4, Ccr2, Il33, Cxcl2, interleukin 1 receptor antagonist (Il1rn), selectin P (Selp), C-C motif chemokine ligand 17 (Ccl17), selectin L (Sell), Ly6c2, selectin E (Sele), and indoleamine 2,3-dioxygenase 1 (Ido1) (Fig. 2F). In the PPI network of hub genes, Ccr2, Ccl4, and Ccl3 displayed prominent interaction patterns, suggesting a core role in the regulatory network (Fig. 2G).

3.3 Identification of THP-DEPs

THP-DEPs were identified via differential expression analysis. QC of the proteomic data confirmed uniform sample distribution and appropriate clustering patterns, indicating reliable data consistency across samples (Supplementary Fig. 4). A total of 106 THP-DEPs were detected, including 40 upregulated and 66 downregulated proteins (Fig. 3A,B). GO analysis yielded 217 enriched terms (Supplementary Table 2), with the most significant categories involving cytoplasm and intracellular anatomical structures (CC), protein binding (e.g., glycosaminoglycan binding), regulation of protein activity (e.g., protein phosphatase inhibitor activity) (MF), cellular component morphogenesis, and postsynapse organization (BP) (Fig. 3C). These results indicate that THP-DEPs may exert their effects predominantly through protein binding and modulation of enzymatic activity. KEGG pathway enrichment revealed significant involvement of differentially expressed proteins (THP-DEPs) in pathways related to mineral absorption, spliceosome function, cardiac muscle contraction, and other associated biological processes (Fig. 3D). PPI network analysis identified 50 highly interactive THP-DEPs, with Rpl4 and Hnrnpc demonstrating the highest interaction strengths (Fig. 3E).

3.4 Metabonomic Analysis

Prior to metabonomic profiling, PCA, PLS-DA, and RPT validation confirmed the robustness and reliability of the dataset (Supplementary Fig. 5). A total of 2387 THP-DEMs were identified, including 1022 upregulated metabolites (e.g., 6-formylumbelliferone, DL-indole-3-lactic acid, benzyl chloride) and 1365 downregulated metabolites (e.g., isopentenyladenosine, uracil, cystine) (Fig. 4A,B). Notably, 7 KEGG pathways were enriched, encompassing eicosanoids (drug development), arachidonic acid metabolism, linoleic acid metabolism (metabolism), intestinal immune network for IgA production, serotonergic synapse, and ovarian steroidogenesis (organismal systems) (Fig. 4C). These data indicate that differentially expressed metabolites are primarily associated with signaling cascades related to organismal systems.

3.5 Ribo-seq Results

Following filtration and normalization of the Ribo-seq dataset, high-quality ribosome profiling results were obtained, as presented in Supplementary Fig. 6. This process yielded 436 RF-DEGs, comprising 355 upregulated and 81 downregulated genes (Fig. 5A,B). GO enrichment analysis revealed 2163 terms (Supplementary Table 3), with RF-DEGs predominantly associated with leukocyte migration and positive regulation of response to external stimuli (BP), collagen-containing extracellular matrix and membrane microdomains (CC), and cell adhesion molecule binding and glycosaminoglycan binding (MF) (Fig. 5C). Notably, the GO terms “positive regulation of response to external stimulus” and “membrane microdomain” were commonly enriched between RF-DEGs and THP-DEGs (Figs. 2C,5C), suggesting a potential convergence of regulatory mechanisms. KEGG pathway enrichment analysis identified 66 enriched pathways, implicating RF-DEGs in lipid metabolism, atherosclerosis, cytokine–cytokine receptor interactions, and human papillomavirus infection (Fig. 5D). Integration of genes from the “ribosome_hub” and “protein_hub” categories with cell-type-specific DEGs led to the identification of 8 key genes: Ftl1, Tpm4, Ccl3, Ccl4, Ccr2, Il33, Cxcl2, and Ly6c2 (Fig. 5E).

3.6 Functional Analysis, Interactions, and Localization of Key Genes

GSEA based on the 8 key genes revealed distinct functional associations. Ccl3, Ccl4, and Cxcl2 were enriched in pathways related to ribosome function, Parkinson’s disease, oxidative phosphorylation, and Huntington’s disease. Ccr2, Il33, Ly6c2, and Tpm4 were linked to cytokine–cytokine receptor interactions, oxidative phosphorylation, steroid biosynthesis, and CAMs (Fig. 6A). The related signaling pathways of key genes are illustrated in Supplementary Fig. 7. Ftl1 exhibited significant enrichment in long-term potentiation, leukocyte transendothelial migration, and complement and coagulation cascades. The recurrent enrichment of oxidative phosphorylation and ribosome-related pathways across multiple genes underscores their roles in translational control and cellular energy homeostasis. The GGI network identified 20 genes exhibiting strong connectivity with the 8 key genes, including Tpm3, Cd1, and Cxcl5 (Fig. 6B). Functional similarity analysis revealed that Ccl4 had the highest similarity score, indicating substantial overlap between its GO annotation terms and those of other genes. This finding suggests that Ccl4 has the highest similarity with other key genes and may play a central regulatory role in THP development (Fig. 6C). Chromosomal mapping localized Cxcl2 to chromosome 5; Ftl1, Tpm4, and Ccr2 to chromosomes 7, 8, and 9, respectively; Ccl3 and Ccl4 to chromosome 11; and Ly6c2 and Il33 to chromosomes 15 and 19, respectively (Fig. 6D).

3.7 Key Gene-disease Network and Drug Prediction Network

The key gene–disease network delineated complex relationships between the eight key genes and various pathological conditions, notably chemical- and drug-induced liver injury, necrosis, and inflammation, underscoring their potential role in THP pathophysiology (Fig. 7A). Drug prediction analysis revealed no pharmacological agents targeting Ftl1 and Ly6c2. However, four compounds (alprostadil, anisomycin, 8-azaguanine, and trichostatin) were predicted to concurrently target four of the remaining key genes, indicating potential synergistic therapeutic value for THP intervention (Fig. 7B).

3.8 Identification of Six Key Cell Types, Along With Analysis of Cellular Communication and Pseudotime for These Key Cells

Single-cell sequencing data encompassing 10 immune-related cell types were further analyzed to identify cell populations associated with THP. Differential expression profiling between THP and control groups identified six key cell types with significant alterations (p 0.05): ECs, astrocytes, brain pericytes, macrophages, OPCs, and neurons (Fig. 8). Intercellular communication analysis revealed a highly interactive network among these key cells. Astrocytes, brain pericytes, and macrophages displayed enhanced communication in the THP group, while ECs maintained consistently strong interaction networks across both conditions (Fig. 9A, Supplementary Tables 4,5). The changes in ligands and receptors between the two groups are shown in Fig. 9B. Individual communication networks were further reconstructed for each key cell type (Supplementary Fig. 8). Ligand–receptor pair analysis identified Nrg3–Erbb4 as a conserved signaling axis mediating intercellular signaling across both experimental groups (Fig. 9C). Pseudotime trajectory analysis provided insights into cellular differentiation and transcriptional dynamics of key genes:

ECs followed a five-stage differentiation path, originating at stage 1 and terminating at stage 4 (Fig. 10A). Ftl1, Il33, and Tpm4 exhibited gradual upregulation, while the other five key genes remained minimally expressed.

Macrophages exhibited a nine-stage trajectory with cells predominantly derived from the THP group. Stage 1 represented early differentiation, while stage 7 marked terminal maturation (Fig. 10B). Ccr2, Cxcl2, Ftl1, Ly6c2, and Tpm4 peaked at early stages; Ccl4 peaked at intermediate stages; and Il33 exhibited progressive upregulation toward maturation.

Pseudotime analysis of the remaining four key cell types (astrocytes, brain pericytes, OPCs, and neurons) revealed comparably complex differentiation dynamics. Notably, Ftl1 and Il33 maintained relatively high expression levels, whereas the remaining six key genes exhibited persistently low transcriptional activity (Supplementary Fig. 9). These findings highlight the context-dependent, cell-specific regulatory roles of key genes in modulating differentiation and intercellular signaling during THP progression.

3.9 Expression Levels of the Key Genes

To investigate the involvement of candidate genes in CPSP, RT-qPCR was performed to validate the expression levels of eight key genes identified through multi-omics integration. Among them, four genes exhibited statistically significant differential expression between the sham and CPSP groups. Specifically, Ftl1 and Ly6c2 were significantly downregulated, while Ccl4 and Ccr2 were significantly upregulated in the thalamic tissue of CPSP mice (Fig. 11A–H). In the transcriptome sequencing data, Ccl3, Ccl4, Ccr2, Cxcl2, Il33, and Ly6c2 were significantly upregulated in the THP group (Fig. 12). The expression trends of Ccl4 and Ccr2 were consistent with the RT-qPCR results.

4. Discussion

In this study, a THP mouse model was established, and thalamic tissue was analyzed using an integrative multi-omics strategy, involving bulk RNA sequencing, single-cell transcriptomics, proteomics, metabolomics, and ribosome profiling. Transcriptomic analysis identified 376 DEGs, which promoted the construction of a PPI network and the selection of 12 hub genes. Parallel omics analyses revealed DEPs, metabolites, and ribosome-associated genes (RF-DEGs), with subsequent enrichment analyses (GO, KEGG) revealing their involvement in immune regulation, oxidative phosphorylation, and ferroptosis. Through integrative analysis, eight key genes were prioritized: Ftl1, Tpm4, Ccl3, Ccl4, Ccr2, Il33, Cxcl2, and Ly6c2. Expression patterns across single-cell datasets highlighted six critical cell populations contributing to THP pathophysiology: ECs, astrocytes, brain pericytes, macrophages, OPCs, and neurons. Intercellular communication networks among these cell types exhibited THP-specific alterations. Pseudotime analysis further revealed distinct differentiation trajectories and dynamic expression patterns of the key genes.

Ftl1, encoding the ferritin light chain, is a central regulator of ferroptosis. Previous study has shown that Ftl1 modulates microglial homeostasis and activation-related genes in ischemic stroke [36]. Mutations in Ftl1 are causally linked to neurodegeneration with brain iron accumulation (NBIA), and its upregulation has been associated with enhanced ferroptotic activity [37]. Notably, DEPs in the THP model were significantly enriched in the ferroptosis pathway (Fig. 3D), suggesting that Ftl1 may serve as a therapeutic target for modulating ferroptosis in THP. To date, no prior studies have reported a direct association between Ftl1 and THP.

Tropomyosin alpha-4 chain (TPM4) is expressed in human aortic endothelial cells and cultured human brain microvascular ECs [38, 39]. It has been identified as a prognostic biomarker in several malignancies. Furthermore, missense mutations in TPM4 disrupt platelet function, impair promyosin cytoskeletal architecture, and are associated with hemorrhagic and thrombotic disorders [38]. These findings suggest a mechanistic link between TPM4 dysfunction and THP, potentially through compromised platelet function and dysregulated immune–endothelial interactions. By influencing inflammatory responses and cellular integrity, TPM4 may contribute to pain pathogenesis. Therefore, therapeutic modulation of TPM4 could enhance platelet stability, attenuate immune-mediated inflammation, and ameliorate symptoms associated with thalamic hemorrhagic pain.

Interleukin-33 (IL-33), a member of the IL-1 cytokine family, is implicated in the regulation of inflammation, host defense, and autoimmunity [40]. It is abundantly expressed in secondary lymphoid organs, such as lymph nodes and appendix, and is distributed throughout the vascular system [41]. IL-33 exerts its immunomodulatory effects through the ST2 receptor (IL-33R), and accumulating evidence supports its neuroprotective role in ischemic stroke models [42, 43]. Notably, IL-33 upregulation has shown to enhance blood–brain barrier integrity and promote functional recovery following stroke [43]. These findings suggest IL-33 as a promising neuroprotective factor and potential therapeutic candidate in THP.

Ccl3, Ccl4, Ccr2, and Cxcl2 are chemokine-related genes with distinct immunological and inflammatory functions. CCL3 plays a pivotal role in the recruitment of diverse immune cell subsets to intra-tumoral microenvironments [44]. In contrast, CCL4 is associated with chronic inflammation, particularly during cellular senescence and vascular dysfunction, and is notably upregulated in endothelial cells. It contributes to senescence and functional decline by promoting reactive oxygen species (ROS) production and inflammatory signaling, potentially exacerbating the pathology of non-healing wounds [45, 46]. Similarly, CCR2 and CXCL2 are implicated in senescence-related pathways and play remarkable roles in the development and progression of cardiovascular diseases. A significant upregulation of CCR2 following subarachnoid hemorrhage is closely associated with neuroinflammation and neuronal apoptosis, while inhibition of CCR2 can effectively alleviate brain edema, decrease blood-brain barrier permeability, and mitigate neuroinflammation, thus improving clinical outcomes [47]. In acute myocardial infarction, CCR2-expressing macrophages possess pro-inflammatory characteristics that promote inflammation and tissue damage. In contrast, M2-polarized macrophages, through the secretion of small extracellular vesicles, regulate CCR2-expressing macrophage function to suppress cardiac inflammation and promote repair, demonstrating the potential of extracellular vesicle-based therapies [48]. Additionally, CXCL2 plays a critical role in age-related liver injury by recruiting neutrophils via the CXCL2-CXCR2 axis, exacerbating liver inflammation. Targeting this pathway may represent a novel therapeutic approach for aging-associated liver injury [49]. These mechanisms highlight the involvement of CCR2 and its associated signaling pathways in both acute and chronic inflammation, suggesting that they may provide novel therapeutic directions for thalamic hemorrhagic pain. The upregulation of CCR2 expression in macrophages may be associated with pain onset, and inhibition of CCR2 may alleviate inflammation-induced pain perception. Therefore, therapeutic strategies targeting CCR2, such as CCR2 antagonists, may provide new treatment approaches for the clinical management of thalamic hemorrhagic pain.

Ly6c2 is a glycosylphosphatidylinositol-anchored cell surface glycoprotein that is widely expressed on immune cells, such as monocytes, macrophages, and dendritic cells [50]. Recent study has highlighted the crucial role of Ly6c2 in various immune-related diseases, particularly in liver diseases, liver fibrosis, and tumor immunity. During liver fibrosis repair, Ly6c2 contributes to fibrosis reversal by modulating the expression level of Tlr4 in CD11b cells and by maintaining a low expression level of Ly6c2 [51]. Additionally, Ly6c2-positive monocytes can differentiate into macrophages and dendritic cells, promoting the formation of memory CD8 T cells, thereby enhancing anti-tumor immune responses [52]. These findings provide novel therapeutic insights for immunotherapy in liver diseases and liver cancer treatment. In study of obesity-induced stroke, Ly6c2 has been significantly upregulated in brain tissues, suggesting its potential role in the immune response during stroke, particularly in regulating immune cells [53]. For the treatment of THP, Ly6c2 may alleviate pain by modulating the immune response in brain tissue. Given the close association between THP and immune inflammation, Ly6c2 may further regulate neuroinflammatory responses triggered by thalamic hemorrhage through influencing the function of immune cells, such as monocytes and macrophages. Therefore, modulating Ly6c2 expression level may provide novel therapeutic strategies for the treatment of THP. By regulating the expression or function of Ly6c2, it may be feasible to reduce immune responses in the brain following stroke, thereby mitigating neuroinflammation and improving pain management and recovery for patients.

Using GSEA, we identified eight key genes potentially involved in the cytokine-cytokine receptor interaction (i.e., a biological process that is closely associated with neuroinflammation and neuropathic pain). Following neural injuries, such as spinal cord injury, cytokine-cytokine receptor interactions activate downstream signaling pathways, promoting inflammatory responses and modulating pain perception. These interactions, particularly in microglia, have shown to play a critical role in neuroinflammation [54]. Given the central role of the thalamus in pain processing, immune dysregulation in this structure may be closely linked to these interactions, further exacerbating pain. The chemokine ligand 5/C-C chemokine receptor type 5 signaling pathway is implicated in neuropathic pain [55]. Although C-C chemokine receptor type 5 and CCR2 bind to different chemokines, they mainly exhibit synergistic effects during immune responses [56]. Therefore, CCR2 may contribute to inflammatory responses in the thalamus by participating in cytokine-cytokine receptor interactions. Additionally, other key genes identified in this study may similarly modulate thalamic immunoregulatory responses, influencing the development and treatment of THP. Further investigation into the precise mechanisms of these genes may provide novel insights into their roles in THP and lay a stronger theoretical basis for clinical interventions.

In this study, scRNA-seq was employed to characterize the thalamic cellular landscape in mice with THP compared to controls, revealing microenvironmental and developmental remodeling. Six principal cell types, including endothelial cells, astrocytes, brain pericytes, macrophages, OPCs, and neurons, were identified. A comprehensive cell–cell interaction network was constructed based on ligand–receptor pair analysis, illustrating the coordinated functional interaction of these cells in THP. Astrocytes, for instance, play multifaceted roles in post-stroke brain repair [57, 58]. Stroke-induced alterations in ECs can lead to the release of extracellular microvesicles that cross the blood–brain barrier and enter the brain parenchyma [59]. These microvesicles may deliver proneural transcription factors, such as Ascl1, reprogramming astrocytes into neural progenitor-like cells [60], thereby implicating a vascular-glial axis in modulating neuroregeneration. Macrophages, commonly classified into pro-inflammatory (M1) and anti-inflammatory (M2) subtypes, are essential in modulating both the initiation and resolution phases of inflammation. Following CPSP-induced brain injury, thalamic recruitment of macrophages has shown to alleviate nociceptive behaviors and neuronal dysfunction in the thalamic–cingulate pathway by attenuating pro-inflammatory activity and cytokine release [61]. Similarly, macrophages play essential roles in resolving neuroinflammation in ischemic stroke models [62, 63, 64].

Several limitations of this study should be acknowledged. The relatively small sample size might constrain the statistical power and generalizability of the results, potentially affecting the precision of key gene and cell-type identification. Additionally, inter-model variability and species-specific differences between animal models and humans may influence the interpretation and translatability of the findings. Future cross-species comparative analyses are warranted to elucidate the extent of conservation and divergence in cellular functions, thereby providing deeper insights into species-specific biological mechanisms. Furthermore, while several DEGs were identified, their mechanistic roles in THP pathogenesis remain to be elucidated and will be the focus of future investigations.

5. Conclusions

In conclusion, this study employed a multi-omics integrative analysis framework to identify and validate key genetic and cellular markers associated with THP. The combined transcriptomic and single-cell analyses elucidated the molecular and cellular architecture of THP, revealed putative gene–cell interaction networks involved in disease progression, and identified candidate targets for therapeutic intervention. These findings advance the understanding of THP pathogenesis and lay the basis for the development of novel diagnostic and therapeutic strategies.

Availability of Data and Materials

The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

References

[1] Virani SS, Alonso A, Benjamin EJ, Bittencourt MS, Callaway CW, Carson AP, et al. Heart Disease and Stroke Statistics-2020 Update: A Report From the American Heart Association. Circulation. 2020; 141: e139–e596. https://doi.org/10.1161/CIR.0000000000000757.

[2] Klit H, Finnerup NB, Jensen TS. Central post-stroke pain: clinical characteristics, pathophysiology, and management. The Lancet. Neurology. 2009; 8: 857–868. https://doi.org/10.1016/S1474-4422(09)70176-0.

[3] Kumar S, Selim MH, Caplan LR. Medical complications after stroke. The Lancet. Neurology. 2010; 9: 105–118. https://doi.org/10.1016/S1474-4422(09)70266-2.

[4] Delpont B, Blanc C, Osseby GV, Hervieu-Bègue M, Giroud M, Béjot Y. Pain after stroke: A review. Revue Neurologique. 2018; 174: 671–674. https://doi.org/10.1016/j.neurol.2017.11.011.

[5] Kumar G, Soni CR. Central post-stroke pain: current evidence. Journal of the Neurological Sciences. 2009; 284: 10–17. https://doi.org/10.1016/j.jns.2009.04.030.

[6] Gurba KN, Chaudhry R, Haroutounian S. Central Neuropathic Pain Syndromes: Current and Emerging Pharmacological Strategies. CNS Drugs. 2022; 36: 483–516. https://doi.org/10.1007/s40263-022-00914-4.

[7] Mulla SM, Wang L, Khokhar R, Izhar Z, Agarwal A, Couban R, et al. Management of Central Poststroke Pain: Systematic Review of Randomized Controlled Trials. Stroke. 2015; 46: 2853–2860. https://doi.org/10.1161/STROKEAHA.115.010259.

[8] Choi HR, Aktas A, Bottros MM. Pharmacotherapy to Manage Central Post-Stroke Pain. CNS Drugs. 2021; 35: 151–160. https://doi.org/10.1007/s40263-021-00791-3.

[9] Fu G, Du S, Huang T, Cao M, Feng X, Wu S, et al. FTO (Fat-Mass and Obesity-Associated Protein) Participates in Hemorrhage-Induced Thalamic Pain by Stabilizing Toll-Like Receptor 4 Expression in Thalamic Neurons. Stroke. 2021; 52: 2393–2403. https://doi.org/10.1161/STROKEAHA.121.034173.

[10] Infantino R, Schiano C, Luongo L, Paino S, Mansueto G, Boccella S, et al. MED1/BDNF/TrkB pathway is involved in thalamic hemorrhage-induced pain and depression by regulating microglia. Neurobiology of Disease. 2022; 164: 105611. https://doi.org/10.1016/j.nbd.2022.105611.

[11] Shi ZM, Jing JJ, Xue ZJ, Chen WJ, Tang YB, Chen DJ, et al. Stellate ganglion block ameliorated central post-stroke pain with comorbid anxiety and depression through inhibiting HIF-1α/NLRP3 signaling following thalamic hemorrhagic stroke. Journal of Neuroinflammation. 2023; 20: 82. https://doi.org/10.1186/s12974-023-02765-2.

[12] Walsh KB, Zhang X, Zhu X, Wohleb E, Woo D, Lu L, et al. Intracerebral Hemorrhage Induces Inflammatory Gene Expression in Peripheral Blood: Global Transcriptional Profiling in Intracerebral Hemorrhage Patients. DNA and Cell Biology. 2019; 38: 660–669. https://doi.org/10.1089/dna.2018.4550.

[13] Walsh KB, Zimmerman KD, Zhang X, Demel SL, Luo Y, Langefeld CD, et al. miR-181a Mediates Inflammatory Gene Expression After Intracerebral Hemorrhage: An Integrated Analysis of miRNA-seq and mRNA-seq in a Swine ICH Model. Journal of Molecular Neuroscience: MN. 2021; 71: 1802–1814. https://doi.org/10.1007/s12031-021-01815-9.

[14] Ma S, Zhang B, LaFave LM, Earl AS, Chiang Z, Hu Y, et al. Chromatin Potential Identified by Shared Single-Cell Profiling of RNA and Chromatin. Cell. 2020; 183: 1103–1116.e20. https://doi.org/10.1016/j.cell.2020.09.056.

[15] Hao Y, Hao S, Andersen-Nissen E, Mauck WM, 3rd, Zheng S, Butler A, et al. Integrated analysis of multimodal single-cell data. Cell. 2021; 184: 3573–3587.e29. https://doi.org/10.1016/j.cell.2021.04.048.

[16] Zhang J, Li Z, Chen Z, Shi W, Xu Y, Huang Z, et al. Comprehensive analysis of macrophage-related genes in prostate cancer by integrated analysis of single-cell and bulk RNA sequencing. Aging. 2024; 16: 6809–6838. https://doi.org/10.18632/aging.205727.

[17] Nip KM, Chiu R, Yang C, Chu J, Mohamadi H, Warren RL, et al. RNA-Bloom enables reference-free and reference-guided sequence assembly for single-cell transcriptomes. Genome Research. 2020; 30: 1191–1200. https://doi.org/10.1101/gr.260174.119.

[18] Arzalluz-Luque Á Conesa A. Single-cell RNAseq for the study of isoforms-how is that possible? Genome Biology. 2018; 19: 110. https://doi.org/10.1186/s13059-018-1496-z.

[19] Huang T, Xiao Y, Zhang Y, Ge Y, Gao J. Combination of single-nucleus and bulk RNA-seq reveals the molecular mechanism of thalamus haemorrhage-induced central poststroke pain. Frontiers in Immunology. 2023; 14: 1174008. https://doi.org/10.3389/fimmu.2023.1174008.

[20] Satija R, Farrell JA, Gennert D, Schier AF, Regev A. Spatial reconstruction of single-cell gene expression data. Nature Biotechnology. 2015; 33: 495–502. https://doi.org/10.1038/nbt.3192.

[21] Liao Y, Smyth GK, Shi W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics (Oxford, England). 2014; 30: 923–930. https://doi.org/10.1093/bioinformatics/btt656.

[22] Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology. 2014; 15: 550. https://doi.org/10.1186/s13059-014-0550-8.

[23] Wu T, Hu E, Xu S, Chen M, Guo P, Dai Z, et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innovation (Cambridge (Mass.)). 2021; 2: 100141. https://doi.org/10.1016/j.xinn.2021.100141.

[24] Zhang X, Chao P, Zhang L, Xu L, Cui X, Wang S, et al. Single-cell RNA and transcriptome sequencing profiles identify immune-associated key genes in the development of diabetic kidney disease. Frontiers in Immunology. 2023; 14: 1030198. https://doi.org/10.3389/fimmu.2023.1030198.

[25] Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Research. 2003; 13: 2498–2504. https://doi.org/10.1101/gr.1239303.

[26] Baker CP, Bruderer R, Abbott J, Arthur JSC, Brenes AJ. Optimizing Spectronaut Search Parameters to Improve Data Quality with Minimal Proteome Coverage Reductions in DIA Analyses of Heterogeneous Samples. Journal of Proteome Research. 2024; 23: 1926–1936. https://doi.org/10.1021/acs.jproteome.3c00671.

[27] Kasyanov ED, Yakovleva YV, Mudrakova TA, Kasyanova AA, Mazo GE. Comorbidity patterns and structure of depressive episodes in patients with bipolar disorder and major depressive disorder. Zhurnal Nevrologii i Psikhiatrii Imeni S.S. Korsakova. 2023; 123: 108–114. https://doi.org/10.17116/jnevro2023123112108.

[28] Rohart F, Gautier B, Singh A, Lê Cao KA. mixOmics: An R package for ’omics feature selection and multiple data integration. PLoS Computational Biology. 2017; 13: e1005752. https://doi.org/10.1371/journal.pcbi.1005752.

[29] Wang Y, Huang J. Untargeted metabolomic analysis of metabolites related to body dysmorphic disorder (BDD). Functional & Integrative Genomics. 2023; 23: 70. https://doi.org/10.1007/s10142-023-00995-4.

[30] Rodriguez-Martinez A, Ayala R, Posma JM, Neves AL, Gauguier D, Nicholson JK, et al. MetaboSignal: a network-based approach for topological analysis of metabotype regulation via metabolic and signaling pathways. Bioinformatics (Oxford, England). 2017; 33: 773–775. https://doi.org/10.1093/bioinformatics/btw697.

[31] Yu G. Gene Ontology Semantic Similarity Analysis Using GOSemSim. Methods in Molecular Biology (Clifton, N.J.). 2020; 2117: 207–215. https://doi.org/10.1007/978-1-0716-0301-7_11.

[32] An J, Lai J, Sajjanhar A, Batra J, Wang C, Nelson CC. J-Circos: an interactive Circos plotter. Bioinformatics (Oxford, England). 2015; 31: 1463–1465. https://doi.org/10.1093/bioinformatics/btu842.

[33] Jin S, Guerrero-Juarez CF, Zhang L, Chang I, Ramos R, Kuan CH, et al. Inference and analysis of cell-cell communication using CellChat. Nature Communications. 2021; 12: 1088. https://doi.org/10.1038/s41467-021-21246-9.

[34] Qiu X, Mao Q, Tang Y, Wang L, Chawla R, Pliner HA, et al. Reversed graph embedding resolves complex single-cell trajectories. Nature Methods. 2017; 14: 979–982. https://doi.org/10.1038/nmeth.4402.

[35] Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) Method. Methods (San Diego, Calif.). 2001; 25: 402–408. https://doi.org/10.1006/meth.2001.1262.

[36] Shi CL, Han XL, Chen JC, Pan QF, Gao YC, Guo PY, et al. Single-nucleus transcriptome unveils the role of ferroptosis in ischemic stroke. Heliyon. 2024; 10: e32727. https://doi.org/10.1016/j.heliyon.2024.e32727.

[37] Hinarejos I, Machuca-Arellano C, Sancho P, Espinós C. Mitochondrial Dysfunction, Oxidative Stress and Neuroinflammation in Neurodegeneration with Brain Iron Accumulation (NBIA). Antioxidants (Basel, Switzerland). 2020; 9: 1020. https://doi.org/10.3390/antiox9101020.

[38] Kimura A, Sakurai T, Yamada M, Koumura A, Hayashi Y, Tanaka Y, et al. Anti-endothelial cell antibodies in patients with cerebral small vessel disease. Current Neurovascular Research. 2012; 9: 296–301. https://doi.org/10.2174/156720212803530726.

[39] Ziegler ME, Souda P, Jin YP, Whitelegge JP, Reed EF. Characterization of the endothelial cell cytoskeleton following HLA class I ligation. PloS One. 2012; 7: e29472. https://doi.org/10.1371/journal.pone.0029472.

[40] Dinarello CA. Biologic basis for interleukin-1 in disease. Blood. 1996; 87: 2095–2147.

[41] Moussion C, Ortega N, Girard JP. The IL-1-like cytokine IL-33 is constitutively expressed in the nucleus of endothelial cells and epithelial cells in vivo: a novel ‘alarmin’? PloS One. 2008; 3: e3331. https://doi.org/10.1371/journal.pone.0003331.

[42] Xie D, Liu H, Xu F, Su W, Ye Q, Yu F, et al. IL33 (Interleukin 33)/ST2 (Interleukin 1 Receptor-Like 1) Axis Drives Protective Microglial Responses and Promotes White Matter Integrity After Stroke. Stroke. 2021; 52: 2150–2161. https://doi.org/10.1161/STROKEAHA.120.032444.

[43] Guo S, Qian C, Li W, Zeng Z, Cai J, Luo Y. Modulation of Neuroinflammation: Advances in Roles and Mechanisms of the IL-33/ST2 Axis Involved in Ischemic Stroke. Neuroimmunomodulation. 2023; 30: 226–236. https://doi.org/10.1159/000533984.

[44] Schaller TH, Batich KA, Suryadevara CM, Desai R, Sampson JH. Chemokines as adjuvants for immunotherapy: implications for immune activation with CCL3. Expert Review of Clinical Immunology. 2017; 13: 1049–1060. https://doi.org/10.1080/1744666X.2017.1384313.

[45] Chang TT, Lin LY, Chen C, Chen JW. CCL4 contributes to aging related angiogenic insufficiency through activating oxidative stress and endothelial inflammation. Angiogenesis. 2024; 27: 475–499. https://doi.org/10.1007/s10456-024-09922-y.

[46] Chang TT, Chen C, Lin LY, Chen JW. CCL4 Deletion Accelerates Wound Healing by Improving Endothelial Cell Functions in Diabetes Mellitus. Biomedicines. 2022; 10: 1963. https://doi.org/10.3390/biomedicines10081963.

[47] Tian Q, Guo Y, Feng S, Liu C, He P, Wang J, et al. Inhibition of CCR2 attenuates neuroinflammation and neuronal apoptosis after subarachnoid hemorrhage through the PI3K/Akt pathway. Journal of Neuroinflammation. 2022; 19: 312. https://doi.org/10.1186/s12974-022-02676-8.

[48] Li L, Cao J, Li S, Cui T, Ni J, Zhang H, et al. M2 Macrophage-Derived sEV Regulate Pro-Inflammatory CCR2⁺ Macrophage Subpopulations to Favor Post-AMI Cardiac Repair. Advanced Science (Weinheim, Baden-Wurttemberg, Germany). 2023; 10: e2202964. https://doi.org/10.1002/advs.202202964.

[49] Liu Y, Xiao J, Cai J, Li R, Sui X, Zhang J, et al. Single-cell immune profiling of mouse liver aging reveals Cxcl2+ macrophages recruit neutrophils to aggravate liver injury. Hepatology (Baltimore, Md.). 2024; 79: 589–605. https://doi.org/10.1097/HEP.0000000000000590.

[50] Morimoto J, Matsumoto M, Miyazawa R, Oya T, Tsuneyama K, Matsumoto M. No Major Impact of Two Homologous Proteins Ly6C1 and Ly6C2 on Immune Homeostasis. ImmunoHorizons. 2022; 6: 202–210. https://doi.org/10.4049/immunohorizons.2100114.

[51] Takimoto Y, Chu PS, Nakamoto N, Hagihara Y, Mikami Y, Miyamoto K, et al. Myeloid TLR4 signaling promotes post-injury withdrawal resolution of murine liver fibrosis. iScience. 2023; 26: 106220. https://doi.org/10.1016/j.isci.2023.106220.

[52] Liu J, Fang C, Jin X, Tian G, Sun Z, Hong L, et al. Nanosecond pulsed electric field ablation-induced modulation of sphingolipid metabolism is associated with Ly6c2⁺ mononuclear phagocyte differentiation in liver cancer. Molecular Oncology. 2023; 17: 1093–1111. https://doi.org/10.1002/1878-0261.13372.

[53] Liang J, Hu R, Wang X, Liu X, Pei L, Tian M, et al. Decoding the Transcriptional Response to Ischemic Stroke in Obese and Non-obese Mice Brain. Current Neurovascular Research. 2021; 18: 211–218. https://doi.org/10.2174/1567202618666210719150845.

[54] Sahebdel F, Zia A, Quinta HR, Morse LR, Olson JK, Battaglino RA. Transcriptomic Profiling of Primary Microglia: Effects of miR-19a-3p and miR-19b-3p on Microglia Activation. International Journal of Molecular Sciences. 2024; 25: 10601. https://doi.org/10.3390/ijms251910601.

[55] Wu S, Yang S, Li R, Ba X, Jiang C, Xiong D, et al. HSV-1 infection-induced herpetic neuralgia involves a CCL5/CCR5-mediated inflammation mechanism. Journal of Medical Virology. 2023; 95: e28718. https://doi.org/10.1002/jmv.28718.

[56] Fantuzzi L, Tagliamonte M, Gauzzi MC, Lopalco L. Dual CCR5/CCR2 targeting: opportunities for the cure of complex disorders. Cellular and Molecular Life Sciences: CMLS. 2019; 76: 4869–4886. https://doi.org/10.1007/s00018-019-03255-6.

[57] Takano T, Oberheim N, Cotrina ML, Nedergaard M. Astrocytes and ischemic injury. Stroke. 2009; 40: S8–12. https://doi.org/10.1161/STROKEAHA.108.533166.

[58] Liu Z, Chopp M. Astrocytes, therapeutic targets for neuroprotection and neurorestoration in ischemic stroke. Progress in Neurobiology. 2016; 144: 103–120. https://doi.org/10.1016/j.pneurobio.2015.09.008.

[59] Zappulli V, Friis KP, Fitzpatrick Z, Maguire CA, Breakefield XO. Extracellular vesicles and intercellular communication within the nervous system. The Journal of Clinical Investigation. 2016; 126: 1198–1207. https://doi.org/10.1172/JCI81134.

[60] Li W, Mandeville ET, Durán-Laforet V, Fukuda N, Yu Z, Zheng Y, et al. Endothelial cells regulate astrocyte to neural progenitor cell trans-differentiation in a mouse model of stroke. Nature Communications. 2022; 13: 7812. https://doi.org/10.1038/s41467-022-35498-6.

[61] Kuan YH, Shih HC, Tang SC, Jeng JS, Shyu BC. Targeting P(2)X(7) receptor for the treatment of central post-stroke pain in a rodent model. Neurobiology of Disease. 2015; 78: 134–145. https://doi.org/10.1016/j.nbd.2015.02.028.

[62] Cai W, Hu M, Li C, Wu R, Lu D, Xie C, et al. FOXP3+ macrophage represses acute ischemic stroke-induced neural inflammation. Autophagy. 2023; 19: 1144–1163. https://doi.org/10.1080/15548627.2022.2116833.

[63] Blank-Stein N, Mass E. Macrophage and monocyte subsets in response to ischemic stroke. European Journal of Immunology. 2023; 53: e2250233. https://doi.org/10.1002/eji.202250233.

[64] Cai W, Dai X, Chen J, Zhao J, Xu M, Zhang L, et al. STAT6/Arg1 promotes microglia/macrophage efferocytosis and inflammation resolution in stroke mice. JCI Insight. 2019; 4: e131355. https://doi.org/10.1172/jci.insight.131355.

Funding

National Natural Science Foundation of China(82172190)

Class A Medical Technology Key Talent of Northern Jiangsu People’s Hospital(JSGG13)

lv Yang Jin Feng program(LYJF00048)

Doctoral Startup Fund of Northern Jiangsu People’s Hospital(BSQDJ0199)

Doctoral Startup Fund of Northern Jiangsu People’s Hospital(BSQDJ0237)

Research Grant of Northern Jiangsu People’s Hospital(YJJ230056)

Yangzhou Natural Science Foundation(YZ2024159)

Yangzhou Talent Program Project(YCJH00021)

Jiangsu Province Traditional Chinese Medicine Science and Technology Development Plan - Young Talent Program(QN202423)