Multi-Omics and Clinical Validation Identify Key Glycolysis- and Immune-Related Genes in Sepsis

Introduction

Sepsis is a severe medical syndrome characterized by a dysregulated host response to infection, leading to organ dysfunction and high mortality. Despite advancements in clinical management, its rising incidence and complex immunopathology remain major global health challenges. Traditional treatment, including the prompt administration of antibiotics, fluid resuscitation, and organ support, are often compromised by diagnostic delays, antimicrobial resistance, and the lack of targeted therapeutic interventions.¹ Therefore, a comprehensive understanding of the molecular mechanisms underlying sepsis is crucial for improving diagnostic accuracy and advancing the development of effective, personalized therapeutic interventions.

Recent studies highlight the pivotal role of metabolic reprogramming, particularly glycolysis, in modulating immune responses during sepsis. Upon activation, immune cells such as macrophages, neutrophils, dendritic cells, and T lymphocytes undergo a rapid metabolic shift from oxidative phosphorylation to aerobic glycolysis, known as the Warburg effect, to fulfill the energetic demands.² This transition enhances cytokine production, phagocytosis, and microbial clearance in the early phase of sepsis.^1,3 However, persistent or dysregulated glycolytic activity may exacerbate inflammation or lead to immunosuppression. In neutrophils, glycolytic inhibition impairs chemotaxis, reactive oxygen species (ROS) production, and neutrophil extracellular trap (NET) formation, compromising the primary host defense.⁴ In macrophages, glycolysis blockade reduces NLRP3 inflammasome activation, phagocytic capacity, and antigen presentation, thereby weakening both innate and adaptive immunity.⁵ Similarly, effector CD4⁺ and CD8⁺ T cells depend on glycolysis to support proliferation and cytokine production; disruption of this pathway results in metabolic paralysis and impaired their ability to mount effective immune responses.^6,7 Given the multifaceted and cell-type–specific roles of glycolysis in immune modulation, a deeper and high-resolution understanding of its molecular regulation in sepsis is urgently needed.

Recent advancements in high-throughput sequencing technologies have significantly enhanced our comprehension of disease mechanisms. Nonetheless, few studies have employed multi-omics approaches to dissect the key molecular mechanisms linking glycolysis to immune remodeling in sepsis. In this study, we integrate bulk transcriptomics with single-cell RNA sequencing (scRNA-seq) to provide a comprehensive characterization of the molecular regulation of glycolysis in sepsis. Causal relationships are deduced through Mendelian randomization, and key findings are experimentally validated using peripheral blood samples from a clinical sepsis cohort. Furthermore, in silico drug docking is conducted to identify potential therapeutic targets. By amalgamating multi-omics profiling, causal inference, and therapeutic prediction, this research elucidates the immunometabolic reprogramming that may drive the progression of sepsis and identifies actionable targets for precision medicine.

Materials and Methods

Data Acquired

Three publicly available transcriptomic datasets were retrieved from the Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/). Two bulk microarray datasets, GSE57065 (training set) and GSE26378 (validation set), were both generated using the GPL570 platform. GSE57065 contains whole blood transcriptomes from 82 patients with sepsis and 25 healthy controls, whereas GSE26378 includes 82 sepsis samples and 21 controls. scRNA-seq data were obtained from GSE175453 (platform: GPL18573), which includes circulating immune cells from 4 individuals with sepsis and 5 healthy donors. Genome-wide association study (GWAS) summary statistics for sepsis (ID: ieu-b-4980) were downloaded from the IEU Open GWAS Project (https://gwas.mrcieu.ac.uk/), encompassing 486,484 individuals of European ancestry (11,463 sepsis cases and 474,841 controls) and 12,243,539 single-nucleotide polymorphisms (SNPs). A curated list of 74 glycolysis marker genes was obtained from the REACTOME_GLYCOLYSIS gene set in the Molecular Signatures Database (MSigDB; https://www.gsea-msigdb.org/gsea/msigdb).

ssGSEA Scoring and WGCNA Module Analysis

Single-sample gene set enrichment analysis (ssGSEA) was performed using the GSVA package⁸ to quantify enrichment scores of 74 glycolysis marker genes. Group differences were assessed by Wilcoxon rank-sum test (P < 0.05). To identify glycolysis-associated gene modules, weighted gene co-expression network analysis (WGCNA)⁹ was applied. Outliers were excluded, and the soft-thresholding power was selected based on scale-free topology. Genes were clustered into modules (minimum size = 100; merge threshold = 0.1). Pearson correlation¹⁰ was then used to evaluate associations between module eigengenes and glycolysis scores (|r| > 0.3, P < 0.05). Genes with module membership (MM) > 0.8 and gene significance (GS) > 0.4 were defined as glycolysis-related module genes.

Differential Expression and Functional Enrichment Analyses

Differentially expressed genes (DEGs) from the GSE57065 dataset were identified using |log2FC| > 0.5 and P < 0.05.¹¹ Expression changes were visualized with volcano plots and heatmaps (ggplot2).¹² Differentially expressed glycolysis-related genes (DGRGs) were defined as the intersection of DEGs and key WGCNA module genes. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were conducted using clusterProfiler.¹³ Protein–protein interaction (PPI) networks were constructed via STRING (confidence score > 0.7).

Mendelian Randomization Analysis

To assess the causal relevance of DGRGs in sepsis, Mendelian randomization (MR) was conducted using genetic variants (SNPs) as instrumental variables (IVs), with DGRGs as exposures and sepsis as the outcome. This analysis assumed that selected SNPs were strongly associated with DGRGs, affected sepsis only via these genes, and were independent of confounders. Expression quantitative trait loci (eQTL) summary statistics were processed using the TwoSampleMR package.¹⁴ SNPs with P < 5 × 10⁻⁶ and LD pruning thresholds (R² < 0.001, distance > 100 kb) were retained. Only SNPs with F-statistics >10 (minimum of 3 SNPs) were included. Exposure and outcome datasets were harmonized to align effect alleles. Causal effects were estimated using five MR methods, with inverse variance weighting (IVW) as the primary. Associations were considered significant at IVW P < 0.05; OR > 1 and < 1 indicated risk and protective effects, respectively.

Machine Learning

PPI analysis of causal genes was performed using the STRING database (confidence score > 0.15), and the corresponding coding genes were considered candidate targets. To further refine these candidates, three machine learning algorithms were applied to the GSE57065 dataset: LASSO (glmnet),¹⁵ SVM-RFE (e1071),¹⁶ and Boruta (Boruta).¹⁷ Genes identified by all three methods were defined as signature genes.

ROC Analysis

Diagnostic performance was evaluated using receiver operating characteristic (ROC) curves generated with the pROC package.¹⁸ Genes with an area under the curve (AUC) > 0.7 were considered to have high diagnostic value.

Immune Cell Infiltration Analysis

The immune landscape of sepsis and control samples was assessed using the CIBERSORT algorithm to estimate the relative proportions of 22 immune cell types.¹⁹ Samples with P > 0.05 were excluded. Group differences in immune cell infiltration were analyzed using the Wilcoxon rank-sum test (P < 0.05). To investigate immune regulatory patterns, Spearman correlation analysis was conducted, with significant correlations defined as |r| > 0.3 and P < 0.05.

Molecular Regulatory Network Construction

The ENCORI (https://rnasysu.com/encori/) and miRnet (http://www.mirnet.ca) databases were employed to predict potential upstream regulatory microRNAs (miRNAs) of hub genes. Furthermore, the ENCODE database (http://www.genome.ucsc.edu/ENCODE/) was utilized to identify potential transcription factors (TFs) of hub genes. Cytoscape software was used to visualize the regulatory network, which includes key miRNAs, mRNAs (key genes), and TFs.

Drug Prediction and Molecular Docking

Potential therapeutic compounds targeting key genes were predicted using the DSigDB database (http://dsigdb.tanlab.org/DSigDBv1.0/). The top five candidates ranked by interaction score were selected as key drug molecules. Drug–gene interaction networks were visualized using Cytoscape. To further explore drug–target interactions, molecular docking was conducted between the highest-ranked compound and its corresponding target genes. Three-dimensional structures of protein receptors were obtained from the Protein Data Bank (PDB, https://www.rcsb.org/), while ligand structures were sourced from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/). Docking simulations were performed using the CB-Dock2 platform (https://cadd.labshare.cn/cb-dock2/php/index.php), and binding affinities were evaluated based on calculated binding energies. Lower energy values indicate more stable interactions and stronger binding affinities between ligands and receptors.

scRNA-Seq Analysis and CellChat Inference

scRNA-seq data were processed using the Seurat package.²⁰ Cells with fewer than 200 or more than 6000 detected genes, or with >5% mitochondrial transcripts, were excluded. Data were log-normalized, and the top 2000 highly variable genes were used for scaling and principal component analysis (PCA). Clustering was visualized using uniform manifold approximation and projection (UMAP). Cell types were annotated based on canonical markers from the literature²¹ and the CellMarker database (http://bio-bigdata.hrbmu.edu.cn/CellMarker); clusters lacking definitive markers were assigned using the SingleR package.²² Differences in cell-type proportions between sepsis and control groups were assessed using chi-square tests (P < 0.05), and populations with significant differences were defined as differentially abundant. Key gene expression was evaluated within these subsets. To explore intercellular communication, CellChat was used to infer signaling networks between cell types (P < 0.05), with representative ligand–receptor pairs visualized via dot plots.²³

Cell Trajectory Analysis

Cell trajectory analysis was performed using the Monocle package,²⁴ based on the top 2000 highly variable genes to reconstruct lineage trajectories across cell subpopulations. This approach enabled pseudotime inference of key cell types and dynamic profiling of key gene expression along differentiation paths.

Patient Recruitment and Sample Processing

Peripheral blood samples were collected from patients with sepsis, diagnosed according to Sepsis-3 criteria,²⁵ within 6 hours of hospital admission and prior to treatment initiation, along with age- and sex-matched healthy controls. Exclusion criteria included age <18 years, pregnancy, malignancy, HIV infection, organ transplantation, or immunosuppressive therapy. One aliquot was used for routine blood counts, and another was processed for RT-qPCR analysis. The study was approved by the Ethics Committee of Sichuan Provincial People’s Hospital (No. 2025510).

RT-qPCR Analysis

Total RNA was isolated from peripheral blood using TRIzol reagent, and reverse transcription was performed with the Takara PrimeScript RT kit (RR037B) to synthesize cDNA. For quantitative PCR (qPCR), 2 µL of cDNA was used per reaction, employing SYBR Green Premix (MCE, HY-K0524) as the detection reagent. Primer sequences (Table 1) were designed and synthesized by Beijing Tsingke Biotech Co., Ltd. Gene expression levels were normalized to GAPDH, and relative quantification was calculated using the 2^−ΔΔCt method. All reactions were performed in triplicate, and expression differences were assessed relative to the control group.

Table 1 RT–qPCR Primer Sequences for Candidate Gene Validation in Sepsis Patients

Statistical Analysis

All bioinformatic analyses were performed using R software (4.2.2). For clinical validation, differences in gene expression between groups were assessed using unpaired two-tailed t-tests. P < 0.05 was considered statistically significant.

Results

Identification and Functional Characterization of Glycolysis-Related DEGs in Sepsis

Glycolysis scores were derived from the expression profiles of 74 glycolysis marker genes. Notably, sepsis patients exhibited significantly lower glycolysis compared to healthy controls in GSE57065 dataset (Figure 1A). Furthermore, to elucidate the molecular mechanisms underlying this metabolic alteration, WGCNA was performed. Hierarchical clustering revealed no outliers (Figure 1B), and a soft-thresholding power of 5 was chosen to achieve scale-free topology (Figure 1C). Finally, a total of 26 co-expression modules were identified (Figure 1D), among which the blue module, which contains 301 genes, exhibited the strongest positive correlation with glycolysis scores (Figure 1E and F) and was therefore selected for further analysis.

Figure 1 Identification and functional characterization of glycolysis-related DEGs in sepsis. (A) Box and violin plots of glycolysis scores between control and sepsis groups. (B) Sample clustering to detect outliers prior to WGCNA. (C) Analysis of scale Independence and mean connectivity to select soft-thresholding power. (D) Gene dendrogram and module identification via WGCNA. (E) Correlation heatmap between module eigengenes and glycolysis scores. (F) Scatter plot of gene significance versus module membership in the blue module. (G) Volcano plot of DEGs between sepsis and controls. (H) Heatmap of DEG expression profiles across samples. (I) Venn diagram showing overlap between DEGs and key module genes. (J) GO enrichment analysis of intersecting genes. (K) Chord plot of the top five KEGG-enriched pathways. (L) PPI network of glycolysis-related DEGs.

Next, 4256 DEGs were identified between sepsis and controls, including 1742 upregulated and 2514 downregulated genes (Figure 1G and H). Intersecting these DEGs with the blue module from WGCNA yielded 267 DGRGs (Figure 1I). GO analysis demonstrated that these DGRGs were significantly enriched in biological processes such as ribonucleoprotein complex biogenesis, preribosome assembly, and RNA catalytic activity (Figure 1J). KEGG pathway analysis indicated significant enrichment in pathways related to nucleocytoplasmic transport and ribosome biogenesis (Figure 1K). A PPI network was also constructed using the STRING database, comprising 150 nodes and 687 edges (Figure 1L). The dense connectivity observed within the network suggests that these genes may participate in tightly coordinated biological processes.

Mendelian Randomization Identifies Causal Genes in Sepsis

In order to identify potential causal regulators of sepsis, a Mendelian randomization analysis was conducted, which identified 22 out of 267 DGRGs with significant causal associations, as detailed in Table 2. Notably, CRY1 and ZFP64 were associated with an increased susceptibility to sepsis (OR > 1), whereas CLNS1A and DDX18 demonstrated protective associations (OR < 1). The corresponding scatter plots are provided in Supplementary Figure 1A–V.

Table 2 Mendelian Randomization Identifies Causal Genes in Sepsis

Machine Learning Identifies and Validates Key Genes in Sepsis

To prioritize diagnostically relevant genes among these causal candidates, we first constructed a PPI network and identified 13 genes with strong protein-level connectivity (Figure 2A). We then applied three machine learning algorithms: LASSO, SVM-RFE, and Boruta, to refine the candidate list and successfully identified eight signature genes consistently selected by all three methods (Figure 2B–F). Of note, DDX18, EIF3L, MAK16, THUMPD1, and ZNF260 were markedly downregulated in sepsis across both datasets (GSE57065 and GSE26378) and exhibited strong diagnostic potential (AUC >0.7; Figure 2G–J).

Figure 2 Continued.

Figure 2 Machine learning identifies and validates key genes in sepsis. (A) PPI network of candidate genes. (B) Selection of the optimal λ in the LASSO model. (C) LASSO coefficients plotted against log(λ). (D) SVM-RFE cross-validation curve indicating the optimal error rate. (E) Feature importance ranked by Boruta analysis. (F) Venn diagram showing overlapping genes identified by LASSO, SVM-RFE, and Boruta. (G and H) Violin plots showing differential expression of hub genes in sepsis vs control samples (GSE57065 and GSE26378). (I and J) ROC curves assessing the diagnostic performance of hub genes in the same datasets.

Immune Landscape Alterations and Hub Gene–Immune Correlations in Sepsis

Given the immune-driven nature of sepsis, we applied the CIBERSORT algorithm to profile the immune landscape. This analysis revealed a distinct shift in cellular composition in sepsis, with increased neutrophils and reduced T cells, B cells, and resting NK cells, relative to control (Figure 3A and B). Correlation analysis showed that resting NK cells were most positively associated with CD8⁺ T cells (r = 0.724), whereas resting CD4⁺ memory T cells exhibited the strongest negative association with M0 macrophages (r = –0.64) (Figure 3C). To evaluate the immunological relevance of the identified hub genes, we assessed their expression in relation to immune cell abundance and found that all five genes were negatively correlated with neutrophils and positively associated with T and naïve B lymphocytes (Figure 3D). These findings imply a potential mechanistic role for these genes in orchestrating immune microenvironmental changes during sepsis.

Figure 3 Immune landscape alterations and hub gene–immune correlations in sepsis. (A) Stacked bar chart showing the relative abundance of 22 immune cell types in sepsis and control samples. (B) Boxplots comparing immune cell infiltration levels between the two groups. (C) Correlation heat map of 22 immune cell. (D) Correlation between five glycolysis-related genes and immune cell infiltration.

Glycolysis Genes and Their Immune Correlates in Sepsis: Clinical Validation

To substantiate the transcriptomic findings, peripheral blood samples from sepsis patients and healthy controls were analyzed. RT-qPCR results confirmed a significant downregulation of all five genes in sepsis (Figure 4A). Concurrently, routine hematological assessments indicated a distinct alteration in leukocyte composition, characterized by increased neutrophil and decreased lymphocyte proportions in sepsis patients (Figure 4B and C). Spearman correlation analysis demonstrated that the expression levels of all five genes were negatively correlated with neutrophil percentage and positively correlated with lymphocyte percentage (Figure 4D). These findings corroborate our transcriptomic data and highlight the potential involvement of these glycolysis-related genes in modulating the immune landscape in sepsis.

Figure 4 Glycolysis genes and their immune correlates in sepsis: clinical validation. (A) mRNA expression levels of five glycolysis-related hub genes in peripheral blood from healthy controls and patients with sepsis (n = 10 per group). (B and C) Percentages of neutrophils and lymphocytes in the two groups. (D) Spearman correlation heatmap illustrating associations between gene expression and immune cell proportions. Data are presented as mean ± SD.

miRNA–mRNA–TF Network Construction and Drug Prediction

The miRNA–mRNA–TF regulatory network, shown in Figure 5A, includes 80 miRNAs, 123 TFs, and 5 hub genes. Key regulators such as GABPA and hsa-miR-1-3p are linked to the hub genes, suggesting potential upstream mechanisms. A drug–gene interaction network was also created, involving three hub genes and five drugs (Figure 5B). Alsterpaullone targets MAK16, DDX18, and THUMPD1 simultaneously, while Colchicine shows the highest interaction with DDX18 and MAK16, with binding affinities of –9.4 kcal/mol and –5.4 kcal/mol, respectively. Key interacting residues were identified as SER-289 and ARG-261 for DDX18, and ARG-280 and ARG-124 for MAK16 (Figure 5C–F). These results provide potential therapeutic options for the treatment of sepsis.

Figure 5 miRNA–mRNA–TF network construction and drug prediction. (A) miRNA–mRNA–transcription factor (TF) regulatory network. (B) Drug–gene interaction network, with green triangles representing predicted drugs and red triangles denoting key target genes. (C and D) Molecular docking of colchicine with DDX18, highlighting the binding pocket and key residues (Ser289 and Arg261). (E and F) Molecular docking of colchicine with MAK16, showing the binding interface and key residues (Arg280 and Arg124).

scRNA-Seq-Based Identification of Key Cell Types and Intercellular Signaling Networks

To explore glycolysis-related mechanisms at the single-cell level, we conducted an analysis of a publicly available scRNA-seq dataset. After quality control and filtering (Supplementary Figure 2A and B), the dataset included 25,591 cells and 18,776 genes. The top 2000 highly variable genes were identified, and the top 10 were labeled (Supplementary Figure 3). PCA and permutation testing validated the selection of 30 components (Supplementary Figures 4 and 5), which were subsequently utilized for UMAP clustering. This approach led to the identification of 22 distinct clusters (Supplementary Figure 6), which were annotated using canonical markers documented in the literature and curated within the CellMarker database. These clusters were categorized into seven primary immune cell types: neutrophils, myeloid cells, T cells, natural killer T (NKT) cells, dendritic cells (DCs), monocytes, and B cells (Figure 6A and B). As illustrated in Figure 6C and D, neutrophils and monocytes exhibited increased levels in sepsis, whereas T cells, myeloid cells, and dendritic cells showed a reduction. To further delineate the expression landscape of glycolysis-related hub genes, we profiled their distribution across all immune cell subsets (Figure 6E). Based on their preferential enrichment in specific immune populations and the established roles of these cells in sepsis, T cells, B cells, and neutrophils were selected for an in-depth analysis. Among the five hub genes, DDX18, EIF3L, and THUMPD1 demonstrated significant intergroup variations in both T cells and neutrophils. In general, these genes were downregulated in cases of sepsis, aligning with findings from bulk microarray data. However, an exception was observed with DDX18, which exhibited elevated expression levels specifically in neutrophils (Figure 6F).

Figure 6 Continued.

Figure 6 Single-cell subpopulation analysis and expression patterns of hub genes. (A) UMAP visualization of annotated immune cell populations. (B) Dot plot showing the expression of canonical marker genes used for cell type identification. Dot size reflects the percentage of expressing cells; color intensity indicates average expression level. (C) Violin plots comparing the relative abundance of immune cell types in sepsis and control samples. (D) Box plots showing differences in immune cell proportions. (E) UMAP plots displaying the expression patterns of five hub genes across all identified cell types. (F) Violin plots comparing hub gene expression in T cells, B cells, and neutrophils between sepsis and control.

Subsequently, we assessed the number and strength of cell–cell interactions among different cell populations (Figure 7A and B). In sepsis, interactions were notably increased between T cells and NK T cells, as well as between neutrophils and myeloid cells. T cell–DC interactions were both more frequent and stronger, while neutrophil–DC interactions occurred more often but with reduced intensity. Interactions between T cells and neutrophils were also elevated. Ligand–receptor analysis further identified a significant involvement of the MIF–(CD74+CD44) axis, particularly in B cell–DC communication across both groups (Figure 7C), indicating its central role in sepsis-related intercellular signaling.

Figure 7 Intercellular cell-cell interactions. (A) Circle plot showing the number of predicted cell–cell interactions; edge width reflects the number of ligand–receptor pairs. (B) Circle plot illustrating interaction strength between cell types; edge width indicates overall interaction intensity. (C) Dot plot of ligand–receptor pairs across cell types. Dot color indicates interaction strength.

Trajectories of Key Cell Types and the Expression Patterns of Glycolysis-Related Key Genes

To elucidate the transcriptional dynamics of key immune cell populations, we conducted a pseudotime trajectory analysis. Initially, T cells and neutrophils were categorized into four and six subclusters, respectively (Figure 8A). The resulting trajectories illustrated dynamic transitions over pseudotime, with darker colors signifying earlier stages of differentiation (Figure 8B). Analysis of these subclusters revealed that subpopulations 1 and 2 were predominant in the terminal stages of T cell and neutrophil differentiation, respectively (Figure 8C). Subsequently, T cells and neutrophils were distributed across seven and five distinct differentiation states (Figure 8D). In the context of sepsis, T cell differentiation was modestly impaired, whereas the overall differentiation of neutrophils remained largely unaffected. Nonetheless, there was a slight increase in the proportion of neutrophil subpopulation 2 during the later stages of differentiation. Furthermore, glycolysis-associated hub genes exhibited distinct transcriptional trajectories during immune cell differentiation. In T cells, DDX18 and EIF3L showed progressive downregulation, while THUMPD1 expression increased (Figure 8E). A similar pattern of DDX18 and EIF3L downregulation was observed in neutrophils, whereas the expression of the remaining genes showed minimal change (Figure 8F). These findings highlight cell type–specific transcriptional dynamics of glycolytic hub genes during immune cell differentiation, suggesting their potential roles in shaping immune cell fate and function in sepsis.

Figure 8 Continued.

Figure 8 Pseudotime trajectories of key immune cell subsets and expression dynamics of glycolysis-related hub genes. (A) UMAP plots showing T cell and neutrophil subtypes, with each color representing a distinct cell subset. (B) Trajectory plots depicting the distribution of cells along pseudotime. (C) Pseudotime analysis colored by clusters, highlighting cluster-specific trajectories. (D) Pseudotime analysis colored by cellular states, illustrating transitions between functional states. (E and F) Expression dynamics of five key glycolysis-related genes across cell states, with expression levels changing along pseudotime.

Discussion

Recent advancements in high-throughput sequencing technologies have significantly enhanced our comprehension of disease mechanisms. Nevertheless, there is a paucity of studies utilizing multi-omics approaches to elucidate the fundamental molecular mechanisms connecting glycolysis with immune remodeling in sepsis. In this study, we observed a marked reduction in glycolytic activity in patients with sepsis. Through integrative transcriptomic analyses, Mendelian randomization, and clinical validation, five glycolysis-related hub genes were identified. These genes exhibited a strong correlation with altered immune cell infiltration, characterized by an increase in neutrophils and a decrease in lymphocytes. scRNA-seq further elucidated changes in immune cell composition and revealed cell type–specific expression patterns of these genes. Together, these findings define a glycolysis-related transcriptional program underlying immune remodeling in sepsis and highlight promising molecular targets for diagnostic and therapeutic intervention.

Previous studies have characterized the immune response in sepsis as a paradoxical state of simultaneous hyperinflammation and immunosuppression. This duality complicates clinical management, as anti-inflammatory therapies may inadvertently exacerbate immune suppression, increasing the risk of secondary infections and delaying recovery.^26,27 Emerging evidence suggests that metabolic reprogramming, especially aberrant glycolysis, plays a pivotal role in driving this immune dysfunction. Inflammatory stimuli prompt both innate and adaptive immune cells to rapidly adopt aerobic glycolysis, even under normoxic conditions, a phenomenon known as the Warburg effect, supporting ATP production, biosynthesis, and accelerated effector functions.²⁸ In macrophages and dendritic cells, the metabolic transition is driven by the PI3K–Akt–mTOR–HIF-1α signaling pathway, which upregulates glucose transporters and glycolytic enzymes.²⁹ T cell subsets exhibit distinct metabolic profiles: effector T cells predominantly rely on glycolysis to sustain proliferation and cytokine secretion, whereas regulatory T cells favor oxidative phosphorylation, with elevated glycolytic activity impairing their suppressive function.³⁰ Glycolysis also underlies trained immunity, wherein epigenetically reprogrammed innate immune cells exhibit heightened responses to subsequent stimuli; this process is facilitated by increased glycolytic flux, which promotes histone modifications and transcriptional remodeling.³¹ However, in pathological conditions such as sepsis, immune cells may exhibit dysregulated glycolytic activity—either excessive or suppressed—contributing to immune exhaustion and increased vulnerability to infection.^32,33 These observations have positioned glycolytic modulation as a promising therapeutic strategy. In macrophages, glycolysis inhibitors such as 2-deoxy-D-glucose and dioscin suppress HIF-1α–mediated gene expression to alleviate hyperinflammation.^34,35 In contrast, promoting glycolysis in lymphocytes may enhance immune competence under immunosuppressive conditions by restoring T cell activity.^7,36 In our study, glycolysis levels were markedly reduced in sepsis patients, consistent with previous reports showing pronounced glycolytic suppression during advanced stages of the disease.^4,37 Moreover, transcriptomic and clinical cohort analyses revealed higher neutrophil counts and lower lymphocyte proportions. These findings are consistent with previous reports that septic patients frequently display an elevated neutrophil-to-lymphocyte ratio (NLR), driven by neutrophilia and lymphopenia—a parameter closely associated with disease severity and poor prognosis in both retrospective analyses and meta-analyses.^38,39 Building on these findings, we performed multi-omics analyses and identified five pivotal glycolysis-related genes, three of which (THUMPD1, MAK16, and DDX18) showed significant associations with T cells and neutrophils.

EIF3L encodes a subunit of the eukaryotic translation initiation factor 3 (EIF3) complex, which is essential for translation initiation and protein synthesis. Previous studies have reported reduced EIF3L expression in T lymphocytes from patients with AIDS, where its downregulation is thought to impair T cell homeostasis and contribute to immunosuppression.⁴⁰ To date, no direct evidence links EIF3L to the regulation of glycolysis. Nonetheless, the participation of other EIF3 subunits, such as EIF3A, EIF3D, and EIF3E, in regulating glycolytic pathways via mRNA-specific translational control implies that EIF3L may participate in a similar regulatory process.^41–43 ZNF260 is a member of the zinc finger protein (ZNF) family. Although its precise role in immune regulation remains undefined, other ZNF family members have been implicated in modulating immune cell infiltration and regulating immunological responses. For example, ZHX2 expression correlates positively with the infiltration of dendritic cells, macrophages, B cells, and T helper cells in gastric cancer.⁴⁴ Similarly, ZFP64 overexpression enhances CD8⁺ T cell infiltration and promotes the secretion of immunosuppressive cytokines.⁴⁵ In addition, pharmacological or genetic upregulation of KLF2 has been shown to alleviate endothelial inflammation and reduce monocyte adhesion in response to serum from patients with COVID-19.⁴⁶ Proper ribosomal RNA (rRNA) processing is essential for ribosome biogenesis, which sustains the elevated protein synthesis required for immune cell activation. Impaired rRNA maturation can therefore limit cytokine production and compromise immune responses.⁴⁷ THUMP domain-containing protein 1 (THUMPD1) acts as an RNA adaptor that promotes the formation of N4-acetylcytidine (ac⁴C) modifications, which stabilize mRNA, optimize its secondary structure, enhance ribosome binding, and facilitate efficient and accurate protein synthesis.⁴⁸ The MAK16 gene encodes a highly conserved nucleolar protein in eukaryotes that participates in key RNA metabolic processes, including cleavage, translocation, localization, editing, and post-transcriptional regulation of mRNA.⁴⁹ Mutations in MAK16 have been linked to G1 phase cell cycle arrest, a critical stage for mRNA and protein synthesis in preparation for cell division.⁵⁰ DDX18, a member of the DEAD-box RNA helicase family, is essential for ribosome biogenesis, facilitating the processing and maturation of 18S rRNA required for ribosomal subunit assembly and function.⁵¹ Although no prior studies have linked THUMPD1, MAK16, or DDX18 to sepsis, their involvement in RNA processing, ribosome biogenesis, and cell cycle regulation implies that their dysregulation may disrupt protein synthesis during immune activation, thereby potentially impacting immune cell function.

Current sepsis management relies primarily on antibiotics and supportive care, with no approved pharmacologic agents that significantly improve outcomes. This highlights the urgent need for novel targeted therapies. Our molecular docking analysis identified colchicine as the top interactor with DDX18 and MAK16, while alsterpaullone was predicted to target DDX18, MAK16, and THUMPD1 simultaneously. Colchicine, widely used for inflammatory disorders such as gout and Behçet’s disease, exerts its effects by inhibiting neutrophil recruitment, chemotaxis, and inflammasome activation.⁵² Recent studies demonstrate its efficacy in sepsis models, where it attenuates acute lung injury, reduces oxidative stress and pyroptosis in alveolar macrophages,⁵³ and alleviates systemic inflammation and hepatic injury.⁵⁴ Colchicine also exhibits antiviral activity in vitro.^55,56 Alsterpaullone, an ATP-competitive inhibitor of cyclin-dependent kinases (CDKs) and glycogen synthase kinase-3β (GSK-3β), has demonstrated antiviral activity against HIV-1 and flaviviruses by disrupting nucleolar localization of viral core proteins.^57,58 It also suppresses Epstein–Barr virus–associated lymphoproliferative disorders in murine models.⁵⁹ These findings support the potential of colchicine and alsterpaullone as therapeutic candidates in sepsis, owing to their combined antiviral and anti-inflammatory activities.

This study has several limitations. First, bulk transcriptomic data may obscure cell type–specific signals, despite partial compensation by scRNA-seq. Second, the scRNA-seq dataset included a limited number of patients, potentially affecting generalizability. Third, although clinical samples confirmed gene expression patterns, functional validation of these genes and predicted drugs was not performed. Further mechanistic and pharmacological studies are needed to support clinical translation.

Conclusion

Sepsis is increasingly recognized as a disorder characterized by intricate metabolic-immune interactions; however, the molecular mechanisms underlying this interplay are not well elucidated. Through integrative multi-omics analysis and clinical validation, we have identified five glycolysis-associated genes that may play a role in this crosstalk. These findings offer mechanistic insights into the pathogenesis of sepsis and propose potential targets for therapeutic intervention.