Evaluating the Diagnostic Value and Molecular Mechanism of Energy Metabolism-Related Gene PEA15 in Sepsis

Introduction

Sepsis, a life-threatening syndrome characterized by acute organ dysfunction due to a dysregulated response to infection, represents a significant global health challenge.¹ With a mortality rate of approximately 22.5%, sepsis accounts for 19.7% of global deaths. It is a disease marked by high incidence, mortality, and treatment costs.² Clinically, sepsis manifests as a systemic inflammatory response and multiple organ dysfunction.³ At the cellular and molecular levels, its pathogenesis is highly complex, involving a range of pathophysiological processes, including inflammatory imbalance, immune dysfunction, mitochondrial damage, coagulopathy, neuroendocrine immune network abnormalities, endoplasmic reticulum stress, and autophagy, ultimately leading to organ failure.⁴ Current treatment strategies, consisting of anti-infection therapy, hemodynamic support, and organ support, form the core of sepsis management.⁵ However, these therapies often provide suboptimal outcomes in heterogeneous patient populations. Incorporating multi-omics analyses to identify precision biomarkers—such as CD33, LY9, and S100A8/A9—for early diagnosis and targeted gene therapy is critical to overcoming this patient variability.^6,7

Energy metabolism (EM), the process of converting nutrients into cellular energy, involves metabolic fuels such as glucose, fatty acids (FAs), amino acids, ketone bodies, and lactate. These fuels are metabolized through a series of interconnected biochemical pathways to generate cellular energy, primarily in the form of adenosine triphosphate (ATP).⁸ Disruptions in EM are a central feature of sepsis and a key driver of organ dysfunction. During the early stages of sepsis, the body enters a “hypermetabolic storm”, leading to a surge in energy consumption (eg, muscle protein breakdown and enhanced gluconeogenesis), while energy utilization efficiency remains low. As mitochondrial dysfunction progresses, ATP production becomes insufficient, contributing to organ failure (eg, heart, liver, and kidneys) due to energy depletion.^9–11

Septic-induced myocardial dysfunction (SIMD) is a common complication of sepsis, characterized by metabolic disorders and myocardial energy failure. Pathologically, SIMD involves altered myocardial substrate uptake, impaired fatty acid oxidation, accumulation of intermediate metabolites, mitochondrial dysfunction, and dysregulated mitophagy/mitochondrial biogenesis, which collectively reduce ATP synthesis and trigger myocardial inflammation and cell death.¹² Notably, in burn-related sepsis, glutamine (Gln) mitigates liver injury by repairing dysfunctional mitochondrial electron transport chain (ETC) complexes, enhancing ATP synthesis, attenuating oxidative stress, and restoring energy supply and redox homeostasis.¹³ In macrophages, TGF-β (transforming growth factor-β) can promote glycolysis, suppress pro-inflammatory cytokine production, and reduce survival in sepsis models.¹⁴ In neutrophils, however, glycolysis inhibition exacerbates immunosuppression and contributes to the progression of sepsis.¹⁵ These findings highlight the complexity of EM’s role in sepsis. Despite substantial research, the identification of validated EM-related biomarkers remains an unmet need.

This study leveraged transcriptomic data from public databases to identify key EM-related genes in sepsis through bioinformatic analysis. Cellular models were established and subjected to sequencing to identify core genes. In parallel, single-cell resolution analysis was conducted to map the expression patterns of these core genes and explore their underlying regulatory mechanisms. Overall, this study provides novel insights for sepsis diagnosis, treatment, and prognostic evaluation.

Materials and Methods

Data Acquisition

Three public transcriptome datasets, GSE65682, GSE95233, and GSE167363, were obtained from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). Among these, GSE65682 and GSE95233 are bulk RNA sequencing (RNA-seq) datasets, while GSE167363 is a single-cell RNA-seq dataset. Additionally, a total of 1110 EMRGs were retrieved from the GeneCards database (https://www.genecards.org/) with a relevance score greater than 2 (supplementary file 1).

Identification of Candidate Genes and Functional Analyses

For the GSE65682 dataset, weighted gene co-expression network analysis (WGCNA) was performed using the WGCNA package (v1.70.3)¹⁶ to identify genes strongly associated with sepsis. The limma package (v3.54.0)¹⁷ was subsequently employed to identify differentially expressed genes (DEGs) between sepsis and control groups (adjusted P < 0.05, |log₂ fold change (FC)| > 1). The ggplot2 (v 3.4.1)¹⁸ and ComplexHeatmap (v 2.14.0)¹⁹ packages were used for visualization. The VennDiagram package (v1.7.1)²⁰ was then used to intersect DEGs, key module genes, and EMRGs, resulting in candidate genes related to EM in sepsis. The potential roles of these candidate genes in sepsis were then further explored, with clusterProfiler package (v 4.2.2)²¹ applied for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses (P < 0.05). For exploring the interactions of candidate genes at the protein level, Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database (https://string-db.org/) was employed, and a protein-protein interaction (PPI) network was established. The proteins with high-quality interactions were visualized by Cytoscape software (v 3.8.2)²² (confidence score > 0.4).

Recognition of EM-Related Key Genes and Diagnostic Value Analysis

In the GSE65682 dataset, machine learning algorithms, including least absolute shrinkage and selection operator (LASSO)via the glmnet package (v4.1.4)²³ and XGBoost via xgboost package (v1.7.3.1),²⁴ were applied to identify EM-related key genes in sepsis. To assess the association between these signature genes and patient survival, Kaplan-Meier (K-M) survival analysis was conducted using the survminer package (v0.4.9),²⁵ with survival differences evaluated by the Log rank test (P < 0.05), leading to the identification of hub genes. Additionally, the expression differences of hub genes between groups in the GSE65682 and GSE95233 datasets were evaluated using the Wilcoxon rank-sum test (P < 0.05) to identify key genes. To assess the diagnostic value of the EM-related key genes for sepsis, receiver operating characteristic (ROC) curves were generated using the pROC package (v1.18.0),²⁶ with an area under the curve (AUC) greater than 0.7 indicating a good diagnostic potential.

Analysis of the Regulatory Mechanism of EM-Related Key Genes

For pathway analysis, Gene Set Enrichment Analysis (GSEA) was conducted within the GSE65682 dataset to explore the signaling pathways in which the EM-related key genes were involved (P < 0.05, |NES| > 1). The reference gene set was “c2.cp.kegg.v2023.1.Hs.symbols.gmt”, which was obtained from Molecular Signatures Database (MSigDB) (https://www.gsea-msigdb.org/gsea/msigdb). Additionally, the NetworkAnalyst database (https://www.networkanalyst.ca/NetworkAnalyst/) was used to predict microRNAs (miRNAs) targeting these key genes. Corresponding long non-coding RNAs (lncRNAs) were predicted via miRNet (https://www.mirnet.ca/miRNet/home.xhtml), starBase (https://rnasysu.com/encori/), and miRcode2 (http://c1.accurascience.com/miRecords/) databases, respectively. Notably, the lncRNAs predicted in all 3 databases were selected as key lncRNAs. Finally, Cytoscape software (v 3.8.2) was applied to illustrate the key lncRNA-miRNA-mRNA regulatory network.

Construction of Cell Models and Functional Analyses

LPS-induced sepsis models were established in BEAS-2B cells (human normal lung epithelial cells) to investigate the underlying mechanisms of this EM-related biomarker at the cellular level. BEAS-2B cells were purchased from Procell Life Science & Technology Co., Ltd. (China). BEAS-2B cells were recovered from passage 2 (P2) and serially passaged 3 times at a 1:2 split ratio to reach passage 5 (P5) for subsequent plating and intervention. The cells were divided into three groups: BEAS-2B normal culture group (the blank control group), empty plasmid transfection group (BEAS-2B+ex-NC), and PEA15 overexpression vector transfection group where BEAS-2B cells were transfected with the pCDH-PEA15-3XFlag-GFP+Puro plasmid (Figure S1) (HonorGene, HG-HO003768, China) (BEAS-2B+ex-PEA15). Forty-eight hours post-transfection, the proportion of positive cells was observed and quantified via fluorescence microscopy to evaluate transfection efficiency. Subsequently, the BEAS-2B+ex-NC group (n=5) and BEAS-2B+ex-PEA15 group (n=5) were each induced with 10 µg/mL LPS (Servicebio, China) to respectively establish the sepsis model group and the PEA15 overexpression sepsis model group. And the blank control group (5 samples), which consisted of normal BEAS-2B cells. Cell viability was evaluated using the Cell Counting Kit-8 (CCK-8) assay (Tongren, Japan) according to the manufacturer’s instructions. Apoptosis was assessed using the Annexin V-PE/7-AAD Apoptosis Detection Kit (Yisheng Biotechnology Co., Ltd., China) following the manufacturer’s protocol. Cytokine levels of tumor necrosis factor-α (TNF-α), interleukin-1β (IL-1β), and interleukin-6 (IL-6) were quantified using corresponding ELISA kits (Jiangsu Meimian Industrial Co., Ltd., China) according to the manufacturer’s guidelines. The results are presented as mean ± standard deviation (SD). Differences between the two groups were analyzed using the unpaired two-tailed Student’s t-test.

Quantitative Real-Time Polymerase Chain Reaction (qRT-PCR)

RNAs were isolated from cells using TRIzol reagent (R401-01, Ambion, USA) and reverse-transcribed into cDNA with the SureScript First-strand cDNA Synthesis Kit (11141ES60, Yisheng, China). RT-qPCR was performed with 2×Universal Blue SYBR Green qPCR Master Mix (G3326-05, Servicebio, China) using primers for the biomarker and internal reference GAPDH. Following the reagent’s instructions for the reaction system and program. Primers for the biomarker and the internal reference gene (GAPDH) are listed in Table S1. The 2-^ΔΔCT method was used to calculate the gene expression level.

Western Blotting (WB)

Total protein was extracted using RIPA lysis buffer (Servicebio, China), and its concentration was determined with a BCA protein assay kit (Beyotime, China) prior to separation by SDS-PAGE electrophoresis and transfer to a PVDF membrane (Millipore, USA). The membrane was blocked with 5% bovine serum albumin (BSA) solution (Servicebio, China; 40 mL) at 37°C for 30 min, then incubated overnight at 4°C with primary antibodies (anti-PEA15, Proteintech, 21446-1-AP, 1:1000; anti-β-actin, Proteintech, 66009-1-Ig, 1:25000, USA) and subsequently with secondary antibodies (HRP-conjugated goat anti-rabbit IgG, Servicebio, GB23303, 1:3000; HRP-conjugated goat anti-mouse IgG, Servicebio, GB23301, 1:5000, China) at room temperature for 60 min, following kit instructions. Finally, the protein bands were visualized with chemiluminescence reagent (Huiying, China).

Functional Analyses of Key DEGs and Acquisition of Core Genes

RNA-seq was performed on cell samples from the biomarker overexpression sepsis model, sepsis model, and blank control groups to capture gene expression profiles (Lianchuan Biotechnology, China). DEGs resulting from biomarker overexpression were identified using the limma package (v3.54.0) (adjusted P < 0.05, |log₂FC| > 1). Low-quality bases (Q < 20) and adapter sequences were removed, and reads shorter than 20 bp were filtered out to ensure high-quality alignment. The reads were aligned to the reference genome to generate BAM files, and alignment statistics for all samples were summarized. Genome coverage distribution plots were created to assess the uniformity of sequencing depth. Using the GTF annotation file (based on the human genome reference version hg38), the number of reads mapped to each gene was counted, generating a “gene × sample” count matrix. Gene expression levels were ranked, and the maximum expression value was retained for genes with multiple expression values. Key DEGs were GO and KEGG analyses for functional exploration using the clusterProfiler package (v4.2.2) (P < 0.05). PPIs of key DEGs were evaluated using the STRING database, and the genes encoding the top 5 proteins ranked by degree in the PPI network (confidence score > 0.9) were selected as core genes.

Functional Analyses and Regulatory Mechanisms of Core Genes

To investigate the molecular mechanisms by which the EM-related biomarker influences sepsis progression, Spearman correlations between the biomarker and core genes were first assessed using the psych package (v 2.1.6) (|cor| > 0.3, P < 0.05). A gene-gene interaction (GGI) network was then constructed via the GeneMANIA database (https://genemania.org/) to predict genes associated with core gene functions and related pathways. Additionally, the GO SemSim package (v 2.24.0)²⁷ was utilized to explore functional correlations among the core genes.

To facilitate further understanding of core genes in sepsis, the RCircos package (v1.22)²⁸ was used for chromosomal localization analysis to identify their specific genomic regions. Subcellular localization information of these core genes was retrieved from the mRNALocater database (http://bio-bigdata.cn/mRNALocater/). In the GSE167363 dataset, the potential mechanisms through which the EM-related biomarker influenced sepsis progression via core genes were explored at the single-cell level. Selection criteria included nFeature RNA (the number of genes quantified per cell) within the range of 200 to 5000, nCount RNA (the total gene counts per cell) under 20,000, and percentmt (the proportion of mitochondrial gene expression per cell) limited to less than 20%. Cell clusters were annotated using literature-derived marker genes.²⁹ Core gene distribution and expression across cell types were assessed to clarify their cellular localization, and cell type functional annotation was conducted using the ReactomeGSA package (v1.16.1)³⁰ to reveal associated biological functions.

Statistical Analysis

Bioinformatics analyses were performed using R software (v 4.2.2). All experimental data are presented as means ± SD. The Wilcoxon rank-sum test, t-test, one-way ANOVA, and chi-square test were applied to evaluate intergroup differences, with statistical significance set at P < 0.05. Data visualization was performed using GraphPad Prism 5 software (v 8.0).³¹

Results

Relevant Functions of EM-Related Candidate Genes and Identification of Signature Genes

In the GSE65682 dataset, after confirming the absence of outlier samples (Figure 1A), WGCNA was performed on all samples to identify key module genes associated with sepsis. The optimal soft thresholding was determined to be 18, based on the scale-free fit index (signed R² = 0.85) and mean connectivity (approaching 0), leading to the establishment of a gene co-expression network. Genes were clustered into modules via hybrid dynamic tree cutting (minModuleSize = 30, mergeCutHeight = 0.4). Pearson correlation analyzed module eigengene (ME)-phenotype (sample grouping) associations. Genes from 14 modules were clustered, excluding those from the grey module (Figure 1B–D). Notably, after screening, 802 key module genes strongly correlated with sample grouping information were identified from the brown module (|cor| = 0.67, P < 0.0001, |GS| > 0.4, |MM| > 0.6) (Figure 1E).

Figure 1 WGCNA analysis and identification of key module genes related to sepsis. (A) Hierarchical clustering of samples in GSE65682. (B) Screening of the optimal soft threshold. Left panel: Scale-free topology fit index (y-axis) across soft-thresholding powers (x-axis); Right panel: Network connectivity. Red line denotes selected index (R2=0.85). β=8 was chosen for meeting scale-free criterion (R2>0.8) and preserving adequate connectivity. (C) Identification of co-expression modules. (D) Correlation heatmap of key modules related to traits and groups. Numbers in each cell indicate correlation and significance. (E) Scatter plot of Gene Significance (GS) vs Module Membership (MM) in the key module.

Furthermore, 1287 DEGs were identified in the sepsis group, comprising 445 upregulated and 842 downregulated genes (adjusted P < 0.05) (Figure 2A and B). Twenty-four candidate genes associated with EM in sepsis were obtained from the intersection of key module genes, DEGs, and EMRGs (Figure 2C). These candidate genes were enriched in GO terms such as “response to gamma radiation” and “regulation of apoptotic signaling pathway” (P < 0.05) (Figure 2D). Additionally, they were linked to KEGG pathways such as “NF-kappa B signaling pathway” and “apoptosis” (P < 0.05) (Figure 2E), highlighting the multiple roles of EM in sepsis progression.

Figure 2 Identification and functional analysis of candidate genes and screening of characteristic genes. (A) Volcano plot showing the distribution of differentially expressed genes (DEGs) between the sepsis and control groups. (B) Heatmap of the distribution of DEGs. (C) Venn diagram of candidate genes. (D) GO enrichment results of candidate genes. (E) KEGG enrichment results of candidate gene. (F) PPI network of the 24 candidate genes. (G) Path diagram of regression coefficients in LASSO analysis of candidate genes. (H) Cross-validation curve of LASSO analysis for candidate genes. (I) Feature gene importance ranked by XGBoost. (J) Venn diagram of signature genes. Abbreviations: BP, Biological Process; CC, Cellular Component; MF, Molecular Function.

The constructed PPI network contained 17 interaction pairs (eg, PARP1-VDAC1), involving 12 proteins, further elucidating the complex protein-level interactions among the candidate genes (Figure 2F). Screening of these candidate genes led to the identification of 7 signature genes from the intersection of 7 LASSO-signature genes (optimal lambda = 0.00271) and 15 XGBoost-signature genes, including PASK, BCL2, FTO, ATM, FXN, PEA15, and EEF2 (Figure 2G–J).

The EM-Related Key Genes PEA15 and ATM Had Outstanding Diagnostic Value for Sepsis

Based on the optimal cutoff values for signature genes (PASK = 2.539, BCL2 = 2.685, FTO = 4.654, ATM = 3.576, FXN = 3.406, PEA15 = 3.812, EEF2 = 8.514), high- and low-expression groups were generated. Significant survival differences between these groups identified ATM and PEA15 as hub genes (P < 0.05), which were associated with patient survival outcomes (Figure 3A–G). Expression analysis revealed consistent expression trends for ATM and PEA15 in both the GSE65682 and GSE95233 datasets, with significant inter-group differences (P < 0.0001), confirming them as key genes (Figure 3H and I). Notably, compared to the control group, expression of these key genes was downregulated in the sepsis group. ROC curves for these genes in both datasets exhibited AUC values exceeding 0.9, demonstrating their excellent diagnostic potential for sepsis (Figure 3J and K).

Figure 3 Identification and validation of key genes. (A–G) Kaplan-Meier survival analysis of feature genes in the high and low expression groups in GSE65682. (A) PASK. (B) BCL2. (C) FTO. (D) ATM. (E) FXN. (F) PEA15. (G) EEF2. (H and I) Validation of expression levels of key genes in (H) GSE65682 and (I) GSE95233. (J and K) ROC curve analysis of key genes (AUC > 0.8) in (J) GSE65682 and (K) GSE95233. The results are presented as means±SD, ****p<0.0001.

PEA15 and ATM Were Regulated by Multiple Regulatory Factors and Affected the Progression of Sepsis Through Multiple Pathways

In the GSE65682 dataset, biological pathways associated with EM-related key genes were explored through GSEA. ATM and PEA15 were significantly linked to critical pathways, including “ribosome”, “spliceosome”, “glycosphingolipid biosynthesis lacto and neolacto series”, “olfactory transduction”, and “neuroactive ligand-receptor interaction” (Figure 4A and B). These results suggest that EM plays a pivotal role in sepsis progression by modulating these pathways.

Figure 4 Crucial functional pathways and detailed molecular regulatory networks of EM-related key genes. (A and B) Single-gene gene set enrichment analysis (GSEA) of hub genes ATM (A) and PEA15 (B). (C) LncRNA-miRNA-mRNA regulatory network analysis of key genes.

Additionally, upstream regulatory elements of EM-related key genes were predicted. A constructed lncRNA-miRNA-mRNA regulatory network revealed 6 miRNAs (eg, hsa-miR-4262 and hsa-miR-17-5p) targeting specific key genes (ATM and PEA15), along with 13 lncRNAs (eg, NEAT1 and XIST) regulating these miRNAs (Figure 4C). Identifying these regulatory elements is crucial for understanding the pathophysiological mechanisms linking EM to sepsis.

Immune Microenvironment Was Altered, and PEA15 and ATM Were Significantly Correlated with Most of Differential Immune Infiltration Cells in Sepsis

In the GSE65682 dataset, differential immune microenvironment profiles between the sepsis and control groups were delineated (Figure 5A). Notably, 25 differential immune cell types were identified, including activated B cells (P < 0.01) (Figure 5B). Significant correlations were observed between most differential immune infiltrating cell types (Figure 5C), with the strongest positive correlation between activated CD8 T cells and effector memory CD8 T cells (cor = 0.83, P < 0.0001). In contrast, activated CD8 T cells exhibited the strongest negative correlation with activated dendritic cells (DCs) (cor = −0.54, P < 0.0001). Moreover, EM-related key genes showed significant correlations with most differential immune infiltrating cell types (Figure 5D). These genes exhibited strong positive correlations with various immune cells, with the strongest correlation observed between PEA15 and T follicular helper cells (cor = 0.617, P < 0.0001). Conversely, EM-related key genes were strongly and inversely correlated with activated DCs (|cor| > 0.3, P < 0.0001). In summary, the immune microenvironment profiles of sepsis and control groups differed, and these differences were associated with significant correlations between EM-related key genes and many immune infiltrating cell types. These immune infiltrating cell types may provide new targets for personalized therapeutic approaches to sepsis.

Figure 5 continued.

Figure 5 Immune cell infiltration analysis. (A) Heatmap of immune cell enrichment scores in the sepsis and control groups. (B) Box plot showing differences in the enrichment scores of 28 immune cell types between the sepsis and control groups. (C) Heatmap showing the correlation between key genes and differentially expressed immune cells. (D) Heatmap of the correlation among differentially infiltrated immune cells. Red and blue squares indicate positive and negative correlations, respectively. The results are presented as means±SD. *p<0.05, **p <0.01, ***p<0.001, ****p<0.0001.

PEA15 Overexpression Sepsis Models Were Successfully Constructed

PEA15, an EM-related biomarker, was selected for further investigation due to its reported association with sepsis-induced acute respiratory distress syndrome.³² To explore its role, a stable human PEA15 overexpression vector was utilized. Successful transfection of the overexpression vector into BEAS-2B cells was confirmed through sequencing and alignment (Figure S2). Following LPS treatment, PEA15 overexpression sepsis models were established, along with a separate sepsis model. Cell samples from both the PEA15 overexpression sepsis model and sepsis model groups were observed under an inverted microscope. The cells exhibited normal morphology and adhered well to the culture surface, confirming their suitability for subsequent experimental analysis (Figure S3).

To validate the successful transfection of the PEA15 overexpression vector into BEAS-2B cells, RT-qPCR and WB assays were performed. PEA15 expression in the PEA15 overexpression sepsis model group was significantly higher compared to the sepsis model group (P < 0.01) (Figure 6A). Consistent results were obtained from WB analysis (P < 0.05) (Figure 6B), further confirming the successful construction of the PEA15 overexpression vector and its efficient transfection into BEAS-2B cells.

Figure 6 Construction and functional validation of the PEA15 overexpression sepsis cell model. A sepsis cell model was established using the BEAS-2B cell line. The sepsis model group transfected with empty vector served as the negative control (NC), while the sepsis cell model with stable PEA15 overexpression was designated as the PEA15 overexpression group (PEA15-OE). (A) Results of real-time quantitative PCR for PEA15 overexpression. (B) Western blot analysis of PEA15 overexpression. (C and D) Apoptosis analysis via flow cytometry. (E) Cell viability measured by CCK-8 assay. (F) Expression levels of TNF-α, IL-1β, and IL-6 detected by ELISA. The results are presented as means±SD. *p<0.05, **p <0.01, ****p<0.0001.

Overexpression of PEA15 Could Enhance Cell Activity, Reduce Apoptosis, and Suppress Inflammation in Sepsis

The cellular functions of PEA15 were assessed. Compared to the sepsis model group, the PEA15 overexpression sepsis model group demonstrated a significant increase in cell viability (P = 0.0001) (Figure 6E). Furthermore, the PEA15 overexpression sepsis model group showed a notable reduction in the percentage of apoptotic cells compared to the sepsis model group (P < 0.0001) (Figure 6C and D). Collectively, these results indicate that PEA15 overexpression enhances cell viability and reduces apoptosis in sepsis-affected cells, suggesting a potential protective role for PEA15 in the pathophysiology of sepsis. Additionally, compared to the sepsis model group, the PEA15 overexpression sepsis model group exhibited significantly lower levels of TNF-α, IL-1β, and IL-6 (P < 0.0001) (Figure 6F). These results suggest that PEA15 overexpression effectively suppresses the inflammatory response in sepsis, offering a promising therapeutic strategy to mitigate sepsis progression.

Related Functional Pathways of Key DEGs and Recognition of Core Genes

RNA-seq was performed to analyze the gene expression profiles of the PEA15 overexpression sepsis model, sepsis model, and blank control groups. A total of 201 DEGs1 were identified in the blank control (vs sepsis model) group, including 40 upregulated and 161 downregulated genes, and 527 DEGs2 were identified in the PEA15 overexpression sepsis model (vs sepsis model) group, consisting of 203 upregulated and 324 downregulated genes (adjusted P < 0.05) (Figure 7A–D). After excluding genes from DEGs2 that shared the same expression trends as DEGs1, 492 key DEGs resulting from PEA15 overexpression were identified (Figure 7E and F).

Figure 7 Functional pathways related to key DEGs and identification of core genes. (A) Volcano plot of differentially expressed genes (DEGs1) between the control and sepsis groups. (B) Heatmap of DEGs1 between the control and sepsis groups. (C) Volcano plot of differentially expressed genes (DEG2) between the PEA15 overexpression and sepsis groups. (D) Heatmap of DEG2 between the PEA15 overexpression and sepsis groups. (E and F) Venn diagrams showing the 492 differentially expressed genes (DEGs1) in the context of PEA15 overexpression, including 199 upregulated (E) and 293 downregulated (F) genes. (G) GO enrichment analysis results of the 492 DEGs. (H) KEGG enrichment analysis results of the 492 DEGs. (I) Identification of five core genes through protein-protein interaction (PPI) network analysis of DEGs1. Abbreviations: BP, Biological Process; CC, Cellular Component; MF, Molecular Function.

Further analysis revealed that the key DEGs were enriched in GO terms such as “extracellular matrix organization” (P < 0.05) (Figure 7G). They were also associated with KEGG pathways, including “hematopoietic cell lineage” (P < 0.05) (Figure 7H). These findings provided insight into the potential regulatory mechanisms of EM-related PEA15 in cell functions and sepsis progression. The shared GO terms and KEGG pathways are presented in Tables S2 and S3, offering a more comprehensive understanding of PEA15’s role in sepsis and highlighting its potential as a biomarker. The PPI network of key DEGs included 205 interaction pairs involving 122 proteins, with RSAD2, IFI44, IFIT3, ISG15, and MX1 identified as core genes (Figure 7I).

PEA15 Could Affect the Progression of Sepsis Through Extensive Interactions with Core Genes

Correlations between the core genes and EM-related PEA15 were visualized in Figure 8A. Notably, ISG15 and MX1 exhibited strong positive correlations with PEA15 (cor > 0.7, P < 0.05). A GGI network of core genes was constructed, highlighting the complex interactions between core genes and 20 related genes (Figure 8B). The physical interactions among these genes were particularly strong, involving functions such as “response to type I interferon” and “response to virus”. This network suggests that PEA15 may influence sepsis progression by modulating these interactions. Additionally, the Friends analysis revealed that IFI44 exhibited the strongest functional correlation with the other core genes, indicating its critical role in the regulatory function of PEA15 in sepsis (Figure 8C).

Figure 8 Correlation, gene-gene interaction (GGI) network, and functional correlations of core genes. (A) Heatmap showing the correlation between PEA15 and core genes. Red and blue represent positive and negative correlations, respectively. (B) Gene-gene interaction network of the five core genes (RSAD2, IFI44, IFIT3, ISG15, MX1). (C) Results of the Friends analysis for the core genes.

Core Genes Were Located on Different Chromosomes and Played a Role in the Cytoplasm

The genomic regions of core genes were mapped, revealing that RSAD2, IFI44, IFIT3, ISG15, and MX1 are located on chromosomes 2, 1, 10, 1, and 21, respectively (Figure 9A). These findings provide a novel perspective on their functional associations at the chromosomal level. Additionally, the core genes were found to be localized in the cytoplasm (Figure 9B), suggesting a potential link between their functions and the cytoplasmic environment in sepsis. This association may be relevant to the regulatory role of EM-related PEA15 in sepsis progression.

Figure 9 Chromosomal and subcellular localization of core genes. (A) Chromosomal locations of the core genes. (B) Subcellular localization of the core genes.

PEA15 Could Affect the Progression of Sepsis by Regulating the Expression Patterns of Core Genes in Different Cell Types and Influencing Cell Functions

At the single-cell level, the molecular mechanisms by which EM-related PEA15 influences sepsis progression were further explored. Initially, the cell types involved in sepsis progression were analyzed. Following the filtration of scRNA-seq data, 30,319 cells and 19,808 genes were retained for further analysis (Figure S4A and B). Variability across these genes was depicted in Figure S5, highlighting 2000 highly variable genes. Based on a scree plot and PCA replacement test, the top 30 principal components were selected (P < 0.05) (Figure 10A and B). After removing batch effects (Figure S6A and B), cell clustering and annotation were performed, resulting in 16 distinct clusters (Figure 10C and D). These clusters were further annotated into 7 cell types using marker genes, including B cells, monocytes, erythroid precursors, T cells, DCs, natural killer (NK) cells, and platelets (Figure 10E–H and Table S4). This comprehensive characterization of cellular diversity in the sepsis environment was crucial for understanding the complexity of the disease. The proportions of cell types in sepsis and control groups were compared, revealing that monocytes represented the largest proportion in both groups (Figure 11A–C). Notably, B cells showed the most significant difference between groups, with a lower proportion in the sepsis group (P < 0.0001) (Table S5). The distribution and expression of core genes across different cell types were examined (Figure 11D–F). RSAD2 was most highly expressed in erythroid precursors, while IFI44 and MX1 showed the highest expression in DCs. IFIT3 and ISG15 exhibited the highest expression in monocytes. To further understand the functional roles of different cell types, related biological pathways were analyzed. These pathways included “NEIL3-mediated resolution of interstrand crosslinks (ICLs)” and “sterol 12-hydroxylation by CYP8B1” (Figure 11G). In summary, PEA15 likely plays a central role in sepsis progression by modulating the expression patterns of core genes across various cell types and influencing their functions. The distinct expression profiles of RSAD2, IFI44, MX1, IFIT3, and ISG15 highlight the complex interactions between different cell types and EM, providing valuable insights into potential therapeutic targets for sepsis management.

Figure 10 Single-cell sequencing analysis of cell types in sepsis. (A) Principal component scree plot. (B) Principal component line chart. (C) UMAP dimensionality reduction clustering. (D) UMAP dimensionality reduction clustering of different groups. (E) Bubble plots showing the expression levels of marker genes in different cell clusters. (F) Bubble plot of marker gene expression levels in different cell types. (G) UMAP diagram of cell type annotation. (H) UMAP diagrams annotated with cell types in different groups.

Figure 11 Identification of key cell types. (A) Number of cells for different cell types in GSE167363. (B) Proportions of different cell types in the control group, sepsis group, and all samples. (C) Proportions of different cell types across different samples. (D) UMAP graph showing the expression levels of key genes. (E) Bubble plot of expression levels of key genes. (F) Heatmap of expression levels of key genes. (G) Pathway enrichment in different cell types.

Discussion

Sepsis is a common and critical condition in clinical practice, characterized by a high incidence and mortality rate.³³ One of its key features is EM disorder.³⁴ Numerous studies have demonstrated that various EMGRs contribute to the onset and progression of sepsis through diverse mechanisms.^13,35–37 Identifying reliable biomarkers related to EM is essential for the prevention and treatment of sepsis. In this study, public RNA-seq data were utilized to identify the EM-related key genes ATM and PEA15 in sepsis through bioinformatics approaches. This study further explored the potential mechanisms of PEA15 by constructing a PEA15 overexpression cell model and conducting single-cell analysis, confirming its diagnostic and therapeutic value in sepsis.

In our study, pathway enrichment analysis of candidate genes identified multiple significantly enriched pathways, with “regulation of apoptotic signaling pathway” emerging as a key component. Accumulating evidence has firmly established apoptosis as a pivotal mechanism governing immune cell death—including neutrophils, macrophages, T lymphocytes, and dendritic cells(DCs)—which often precedes the onset of multiple organ failure in sepsis.³⁸ Previous studies have shown that sepsis induces apoptosis and mitophagy in DCs, with mitophagy exerting a protective effect against sepsis-induced DC apoptosis.³⁹ Consistent with this established paradigm, prior investigations have demonstrated that inhibition of apoptosis enhances survival rates in animal models of severe sepsis, while anti-apoptotic signaling pathways have been shown to mitigate sepsis-induced myocardial dysfunction.⁴⁰ These results indicate the complexity of energy metabolism gene regulation in sepsis. We infer that it may affect the changes in immune cells and the immune microenvironment through different mechanisms, thereby influencing the progression of the disease.

PEA15 is widely expressed in the nervous system, especially in astrocytes, and plays a role in regulating cellular processes such as signal transduction, proliferation, and apoptosis.⁴¹ PEP-1-PEA15 can inhibit inflammatory responses by modulating macrophages and the MAPK pathway. Mechanistically, PEA15 interacts with mitochondrial hexokinase II (HKII) to form a molecular switch that governs cell fate under varying metabolic conditions, determining whether cells follow anabolic pathways for energy storage or enter catabolic states to produce energy.⁴² Furthermore, accumulating evidence links PEA15 to diverse pathological processes beyond sepsis. For instance, the TNF-α-LPA-LPA receptor-PLCε-PKD-PEA15-RSK-IκB-NF-κB pathway plays a pivotal role in promoting colonic inflammation,⁴³ reflecting PEA15’s conserved involvement in pro-inflammatory cascades. Additionally, PEA15 interacts with hexokinase II (HKII) to inhibit hypoxia-induced apoptosis and regulate cellular metabolism,⁴² underscoring its dual function in cell survival and metabolic homeostasis. Building on these insights, our study addresses a critical research gap by systematically exploring PEA15’s role in sepsis for the first time. Initially identified via energy metabolism-focused screening—consistent with its known metabolic regulatory activity⁴²—PEA15 correlates significantly with sepsis patient prognosis and exhibits robust discriminatory power between sepsis and control samples. Functionally, PEA15 overexpression in cell models effectively reduces pro-inflammatory cytokine levels. Collectively, these findings provide the first comprehensive demonstration of PEA15’s “diagnostic-prognostic-functional” triple value in sepsis, highlighting its potential as a clinically relevant biomarker and therapeutic target for this life-threatening disease.

ATM is pivotal in regulating the cell cycle, DNA repair, apoptosis, and antioxidant responses. In response to DNA double-strand breaks, ATM activates downstream signaling pathways to maintain genomic stability.⁴⁴ In sepsis, ATM influences immune recovery and immunosuppression by regulating oxidative stress, DNA repair, and immune responses.⁴⁵ Previous studies have demonstrated significant downregulation of ATM gene expression in sepsis patients.⁴⁶ Notably, reduced ATM expression in HK-2 cells attenuates lipopolysaccharide (LPS)-induced inflammation and autophagy in sepsis-associated acute kidney injury (SA-AKI) in vitro, supporting that LPS may induce HK-2 cell autophagy via the ATM pathway to promote pro-inflammatory cytokine upregulation.⁴⁷ In our study, both PEA15 and ATM were downregulated in sepsis, consistent with previous findings. Overexpression of PEA15 enhanced cell viability, reduced apoptosis, and suppressed inflammation.

However, this study has several limitations. First, the utilized datasets are limited in sample size, and machine learning approaches have inherent constraints (eg, overfitting); the predictive performance of the ceRNA network also depends on existing data and algorithms, potentially leading to discrepancies between predicted and actual outcomes. Second, only the BEAS-2B cell line was employed, which may not fully recapitulate the biological complexity of sepsis. Third, the functional role of ATM remains experimentally unvalidated, and the association between PEA15 and EM relies solely on gene expression correlations, lacking direct metabolic validation data. Fourth, the clinical applicability of the identified key sepsis-related genes requires further robust evidence, necessitating additional in-depth investigations. In future work, we plan to collaborate with clinical institutions to conduct large-sample, multi-center prospective studies, validate the diagnostic efficacy of PEA15, and establish clinical application guidelines to support its eventual translation into clinical practice.

Conclusions

In conclusion, this study identified two key EM-related genes (ATM and PEA15) with strong diagnostic potential in sepsis. The constructed cell model and subsequent single-cell analysis demonstrated that PEA15 has the potential to serve as a candidate biomarker for sepsis, though further validation is required. These findings provide a new theoretical basis and direction for the early diagnosis and treatment of sepsis, as well as for exploring the role of energy metabolism in sepsis.