The study utilized datasets obtained from the Gene Expression Omnibus (GEO) website (https://www.ncbi.nlm.nih.gov/geo/). The training set, GSE151371 (GPL20301), included sequencing data from 38 SCI samples and 10 healthy controls. The validation set, GSE47681 (GPL1261), consisted of sequencing data from mouse spinal cord tissue, comprising nine control mice, nine mice on day 3 post-SCI, and nine mice on day 7 post-SCI. A total of 10 histone lactylation modification-related genes (HLMRGs) were analyzed: HDAC1, HDAC2, HDAC3, SIRT1, SIRT2, SIRT3 (Moreno-Yruela et al., 2022), LDHA, LDHB (Yu et al., 2021); P300 (EP300) (Zhang et al., 2019), and GCN5 (KAT2A) (Wang et al., 2022).

To assess the expression of HLMRGs in SCI, a Wilcoxon test was first applied to compare the expression between SCI and control samples in the GSE151371 dataset. Additionally, a heatmap of the expression profiles of the 10 HLMRGs across all samples was generated using the Pheatmap package (v1.0.12, https://CRAN.R-project.org/package=pheatmap). Correlations between the 10 HLMRGs were analyzed using Spearman’s correlation. The RCircos package (v1.2.2) (Zhang et al., 2013) was used to map the 10 HLMRGs to their respective chromosomes.

The gene expression data was log-transformed (log2 (x+1)) to stabilize variance and meet the normality assumption. Samples and genes with missing values were filtered using the na. omit () method. Logistic regression and SVM classification were applied, with regularization to prevent overfitting. The model’s performance on imbalanced data was evaluated using the area under curve (AUC)-receiver operating characteristic (ROC) curve based on assessment metrics. The ROC curves for differentially expressed HLMRGs (DE-HLMRGs) in SCI and control samples from GSE151371 were plotted using the pROC package (v1.1.1) (Robin et al., 2011). Genes with AUC greater than 0.7 were considered as potential biomarkers for constructing a logistic regression model, and the model’s performance was evaluated using the ROC curve. The createDataPartition () method was used for stratified sampling by default, ensuring that the ratio of the training and testing sets matched the original data distribution. Finally, 10-fold cross-validation was performed using trainControl to reduce the randomness of a single data split and enhance the reliability of model stability evaluation.

Based on the identified biomarkers in the GSE151371 dataset, missing values were removed, and a logistic regression model was fitted using the lrm () function. The data distribution was then set using the datadist () function for subsequent analysis. Next, the nomogram () function from the R package rms (v1.6) (https://CRAN.R-project.org/package=rms) was used to generate a nomogram, specifying the transformation function for the predicted probabilities. The plot () function was subsequently used to visualize the nomogram to assess the clinical applicability of the model. Finally, the calibrate () function was used to perform calibration analysis, evaluating the consistency between the predicted probabilities and the actual observed probabilities. The C-statistic was used to evaluate the calibration curve, referring to the proportion of all possible positive-negative sample pairs in which the model predicted the positive samples with higher probabilities. This statistic effectively measured the model’s predictive ability and accuracy. When the C-statistic value was greater than 0.7, it indicated a stronger correlation and better performance of the model. Decision curve analysis (DCA) was performed using the rms package (v1.6) to assess whether the nomogram prediction was beneficial to patients with SCI.

To explore the biological functions and signaling pathways associated with the biomarkers, the correlation coefficients between the expression of all genes and biomarkers were calculated. Genes were then ranked based on their correlation coefficients for GSEA on each biomarker, applying the criteria | Normalized Enrichment Score (NES)| > 1, p < 0.05, and False Discovery Rate (FDR) < 0.25. Single-gene GSEA for each biomarker was conducted using the ClusterProfiler package (v3.18.1) (Yu et al., 2012) in GSE151371. The top five most significant pathways in the Kyoto Encyclopedia of Genes and Genomes (KEGG), Gene Ontology (GO) Biological Process (BP), GO Cellular Component (CC), and GO Molecular Function (MF) for each biomarker were presented.

For immune cell composition analysis, the GSE151371 samples were examined for immune infiltration using CIBERSORT, which provided the distributional proportions of various immune cell types. The proportions of immune cells between the control and SCI groups were compared using the Wilcoxon test (P < 0.05). The p-values were not corrected. The correlation of biomarkers with differential immune cell types was analyzed by Spearman correlation at P < 0.05 and |correlation (r)| > 0.3. The “HALLMARK_INFLAMMATORY_RESPONSE” was searched in the MSigDB database (https://www.gsea-msigdb.org/), and inflammation-related genes were obtained (Lu and Jia, 2022). The correlation between inflammation-related genes and differential immune cells was calculated through Spearman analysis (P < 0.05 and |correlation (r)| > 0.3). Genes were grouped, and the maximum absolute correlation was computed. The top 10 genes with the highest correlation were selected for display.

Based on the identified biomarkers, unsupervised consistent cluster analysis was performed on the 38 SCI samples from GSE151371 to determine distinct clusters. Principal component analysis (PCA) was then conducted to evaluate the differentiation ability of the clusters. Wilcoxon’s method was applied to assess differences in biomarker expression between the clusters.

Differentially expressed genes (DEGs) between the clusters were identified using the limma package (version 3.56.2) (Ritchie et al., 2015), with criteria of adj. P < 0.05 and |log₂foldchange (FC)| > 0.5. GO and KEGG analyses of these DEGs were performed using the ClusterProfiler package (P < 0.05) to explore the functional roles of the DEGs. The P-values were adjusted using the Benjamini–Hochberg (BH) correction method. In addition, the protein-protein interaction network diagram between the DEGs and HLMRGs was constructed using the Search Tool for the Retrieval of Interacting Genes (STRING) database (https://string-db.org/), showcasing the interactions between the DEGs and HLMRGs. The Confidence score was set to 0.9, the species was selected as “Homo sapiens,” and Cytoscape (v3.10.1) (Doncheva et al., 2019) was used for data visualization.

Gene set variation analysis (GSVA) was employed to calculate the scores of Hallmark pathway gene sets using the GSVA package (v1.38.2) (Hanzelmann et al., 2013). Wilcoxon’s test was applied to compare the differences in pathway activity across clusters at P < 0.05. The p-values were not corrected. Additionally, differences in immune cell infiltration between the clusters were analyzed using Wilcoxon’s test (P < 0.05), and the correlation of differential immune cells with biomarkers was examined. The P-values were adjusted using the Benjamini–Hochberg (BH) correction method. Differences in immune function pathways and HLA genes (He et al., 2018) were also analyzed using Wilcoxon’s test, and the correlation between the 10 HLMRGs and immune function pathways or HLA genes was assessed by Spearman correlation analysis.

To identify gene modules most strongly associated with the clusters, weighted gene co-expression network analysis (WGCNA) was performed using the WGCNA package (v1.72–1) (Langfelder and Horvath, 2008). Initially, all samples from different clusters were included in a sample clustering tree to remove outliers. A soft threshold (β) was determined based on the nearly scale-free topology criterion. The resulting topology matrix was clustered using gene differences (minModuleSize = 100, TOMType = “unsigned”, mergeCutHeight = 0.45, power = 16). The dynamic tree-cutting algorithm was used to partition the tree into modules, with different colors representing distinct gene modules and gray representing genes not assigned to any module. A heatmap of the relationships between module traits and clusters was plotted to evaluate the association of each module with different clusters. The most relevant module genes, with P < 0.05, were selected as key modules. Genes with |module membership (MM)| > 0.8 and |gene significance (GS)| > 0.2 were considered key module genes for constructing a protein-protein interaction (PPI) network. Protein interaction information was obtained from the Search Tool for the Retrieval of Interacting Genes (STRING) database (https://string-db.org/) to construct a PPI network. The Cytohubba plug-in (MCC, top = 10) was applied to identify the top 10 key proteins, which were recorded as key proteins for the clusters. The correlations between these 10 key proteins and biomarkers were analyzed using Spearman correlation analysis (P < 0.05 and |r| > 0.3).

Transcription factors (TFs) for the biomarkers were predicted using the ChEA3 database. The miRNA interactions for the biomarkers were predicted using miRWalk (http://mirwalk.umm.uni-heidelberg.de/), Encyclopedia of RNA Interactome (ENOCRI, http://starbase.sysu.edu.cn/index.php), and miRTarBase (https://mirtarbase.cuhk.edu.cn/) databases. The common miRNAs identified through these databases were used to construct a miRNA-TF-biomarker interaction network. Additionally, potential drugs targeting the biomarkers were predicted using the Drug-Gene Interaction (DGI) database.

To validate the expression of biomarkers, SCI samples from GSE47681 were grouped separately into SCI and control categories for comparison. This study adhered to the principles outlined in the Declaration of Helsinki and was approved by the Ethics Committee of The Second Hospital of Anhui Medical University (Approval No. YW 2023-118, Date: 2023.07.06). All patients provided informed consent. To further verify the biomarker expression, reverse transcription-quantitative polymerase chain reaction (RT-qPCR) was performed on tissues from 5 SCI samples and five control samples. Clinical information of the patients is presented in Supplementary Table S1. RT-qPCR amplification was conducted for 40 cycles, consisting of 95°C for 1 min, 95°C for 20 s, 55°C for 20 s, and 72°C for 30 s. The primers used in RT-qPCR are shown in Table 1, with GAPDH as the reference gene. The relative expression levels of the biomarkers were calculated using the 2^−ΔΔCT method.

All statistical analyses were performed using R software (v4.2.2). Differences between groups were assessed using Wilcoxon’s test, with P < 0.05 considered statistically significant.

In GSE151371, eight DE-HLMRGs exhibited significant expression differences between SCI and control samples (Figure 1A). These genes included HDAC1, HDAC2, HDAC3, SIRT1, SIRT3, LDHA, LDHB, and GCN5 (KAT2A). LDHA was upregulated in SCI samples compared to controls, while the remaining genes were downregulated. The heatmap of the expression profiles of the 10 HLMRGs in both groups is shown in Figure 1B. Correlation analysis revealed significant positive correlations among most of the HLMRGs. The strongest positive correlation was observed between HDAC2 and SIRT1 (r = 0.811, P < 0.01), while the strongest negative correlation was found between LDHA and KAT2A (r = −0.489, P < 0.01) (Figure 1C). According to the chromosome localization map, HDAC1 was located on chromosome 1, HDAC3 on chromosome 5, and HDAC2 on chromosome 6, among others (Figure 1D).

The ROC curves for the eight DE-HLMRGs (biomarkers) showed AUC values greater than 0.7, indicating their potential to differentiate between SCI and control samples (Figure 2A). Furthermore, the ROC curve for the logistic regression model yielded an AUC of 0.979, suggesting excellent performance of the model (Figure 2B). The results from 10-fold cross-validation showed that the model exhibited excellent classification performance (AUC = 0.989) (Supplementary Figure S1).

A nomogram model was constructed based on the biomarkers from GSE151371 to predict SCI (Figure 3A). The calibration curve had a C-statistic of 0.987, indicating strong concordance between the apparent and predicted values, and suggesting the model’s high accuracy (Figure 3B). DCA demonstrated that the nomogram model offered a higher overall benefit to patients with SCI compared to individual biomarkers (Figure 3C).

For functional analysis, HDAC1 was enriched in 37 KEGG pathways (Figure 4A), 917 GO-BP (Figure 4B), 209 GO-CC (Figure 4C), and 218 GO-MF (Figure 4D) by GSEA, including pathways such as ‘ubiquitin-mediated proteolysis’, ‘oxidative phosphorylation’, and ‘valine, leucine, and isoleucine degradation’. The enrichment results for the remaining seven biomarkers (HDAC2, HDAC3, SIRT1, SIRT3, LDHA, LDHB, and GCN5 [KAT2A]) are shown in Supplementary Figures S2-8.

The immune cell composition in each sample is illustrated in Figure 5A. Between the SCI and control groups, 15 immune cell types showed significant differences, including memory B cells, eosinophils, and M0 macrophages (Figure 5B). Spearman correlation analysis revealed the differential immune cells most significantly associated with the biomarkers, including HDAC1-memory B cells (r = 0.440, P < 0.01), HDAC2-memory B cells (r = 0.510, P < 0.01), HDAC3-neutrophils (r = −0.310, P < 0.05), KAT2A-naive CD4 T cells/macrophages (both r = 0.490, P < 0.01), LDHA-plasma cells (r = 0.470, P < 0.01), LDHB-resting memory CD4 T cells (r = 0.750, P < 0.01), SIRT1-resting memory CD4 T cells (r = 0.640, P < 0.01), and SIRT3-neutrophils (r = −0.390, P < 0.01) (Figure 5C). Monocytes were significantly positively correlated with LPAR1 and MYC, while CD4 memory activated T cells were significantly negatively correlated with LPAR1 and MYC. The inflammatory factor IL7R was negatively correlated with T cells gamma delta and significantly positively correlated with activated dendritic cells (Supplementary Figure S9).

The 38 SCI samples from GSE151371 were categorized into two clusters based on biomarkers (Figure 6A). PCA indicated distinct differentiation between the two clusters (Figure 6B). Cluster two exhibited higher expression levels of all biomarkers, with significant differences observed only for HDAC1, HDAC3, GCN5 (KAT2A), and SIRT3 (P < 0.05) between the two groups (Figure 6C). The inter-cluster expression of the eight biomarkers is depicted in the heatmap (Figure 6D).

A total of 700 DEGs were identified between clusters 1 and 2, including 287 upregulated and 413 downregulated genes (Figure 7A). These DEGs were enriched in 140 GO terms and four KEGG pathways, including ‘Base excision repair’, ‘phosphatidylinositol signaling system’, and ‘inositol phosphate metabolism’ (Figure 7B). In addition, the PPI network diagram between the 700 DEGs and the eight HLMRGs is shown in Figure 7C.

Seven pathways showed significant differences between the clusters, including ‘apical junction’, ‘apical surface’, and ‘estrogen response early’ in the Hallmark set (Figure 8A). Four immune cell types, including activated dendritic cells, monocytes, activated CD4 memory T cells, and gamma delta T cells, displayed differential expression between the clusters (Figure 8B). Monocytes were positively correlated with HDAC3 and KAT2A, while activated dendritic cells showed a significant negative correlation with HDAC3 and SIRT3 (Figure 8C). Differences in cytolytic activity and parainflammation were also observed between the clusters (Figure 8D). Expression differences of HLA genes, including HLA-A, HLA-B, HLA-DMA, HLA-DMB, HLA-DPB1, HLA-DQB2, HLA-DRA, and HLA-E, were noted across the clusters (Figure 8E). Regarding immune function, strong positive correlations were observed between HDAC2 and APC co-stimulation (r = 0.690, P < 0.01), as well as between LDHB and HLA (r = 0.690, P < 0.01) (Figure 8F). Conversely, a strong negative correlation was found between LDHA and Type I IFN Response (r = −0.730, P < 0.01). The most significant positive correlation for HLA was observed between LDHB and HLA-DQA1 (r = 0.730, P < 0.01) (Figure 8G), while the strongest negative correlation was between LDHA and HLA-DMB (r = −0.430, P < 0.01).

For WGCNA, no outlier samples were identified within the clusters (Figure 9A). The soft threshold (β) was set to 16, which closely approximated a scale-free network (Figure 9B). Subsequently, six gene modules were identified through dynamic tree cutting (Figure 9C). Correlation analysis revealed that the blue module (r = −0.600, P < 0.05) was the most strongly associated with the clusters and was considered a key module (Figure 9D). A total of 1,615 module genes were selected for further analysis (Figure 9E). Based on the PPI network, 10 key proteins were identified: RPS9, RPS28, RPS2, RPS15, RPL13, FAU, RPL18, RPL36, RPL28, and RPL8 (Figure 9F). Correlation analysis showed that KAT2A (GCN5) was most strongly positively correlated with RPL28 (r = 0.797, P < 0.01), while LDHA showed the strongest negative correlation with RPL28 (r = −0.478, P < 0.01) (Figure 9G).

A total of 41 TFs associated with the biomarkers were predicted using the ChEA3 database, from which the top 10 TFs were selected for network construction: ZNF883, ZNF614, ZNF644, ZNF280C, TP53, ZNF140, RBPJ, ZNF692, E2F4, and USF2. Additionally, 57 common miRNAs were predicted using miRWalk, ENOCRI, and miRTarBase databases. A miRNA-TF-biomarker interaction network was constructed based on the eight biomarkers, the top 10 TFs, and the 57 common miRNAs, revealing complex regulatory relationships (Figure 10A). Using the DGIdb database, 47 potential drugs were predicted for HDAC1, 39 for HDAC2, 38 for HDAC3, 24 for SIRT1, three for SIRT3, one for LDHA, and none for LDHB. GCN5 (KAT2A) was predicted to interact with 126 potential drugs. A drug-biomarker interaction network comprising 279 edges and 208 nodes was created, including interactions such as niacinamide-SIRT3, suramin-SIRT3, and lapachone-SIRT1 (Figure 10B).

In the GSE151371 dataset (Figure 1A), only LDHA was highly expressed in the SCI group. In the GSE47681 dataset, HDAC1, GCN5 (KAT2A), LDHA, and SIRT1 were highly expressed in SCI samples, whereas HDAC2, HDAC3, LDHB, and SIRT3 were downregulated in SCI samples (Figure 11A). RT-qPCR validation revealed that HDAC1 (not significant), HDAC2, GCN5 (KAT2A), LDHA, and SIRT1 (not significant) were highly expressed, while HDAC3, LDHB (not significant), and SIRT3 were downregulated in SCI samples (Figure 11B). Although HDAC1, SIRT1, and LDHB did not reach statistical significance, their expression trends were consistent with the results of the computational analysis in the dataset and require further validation. Among the eight biomarkers, LDHA exhibited high expression in SCI, while HDAC3, LDHB, and SIRT3 were downregulated.