The GSE65682, GSE13904, and GSE26378 datasets were retrieved from the Gene Expression Omnibus (GEO) database (www.ncbi.nlm.nih.gov/geo/). The GSE65682 dataset contained 802 blood samples: 760 sepsis and 42 healthy volunteers based on the GPL3667 platform. The GSE13904 and GSE26378, covering pediatric sepsis gene expression data, were used for external validation. The former dataset included 18 healthy children, 52 sepsis subjects, and 106 subjects with septic shock based on the GPL570 platform. The latter comprised 21 normal children and 82 sepsis cases, also based on the GPL570 platform (Table 1). The profiles of apoptosis-related genes were obtained from the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database (https://www.genome.jp/entry/map04210). Then, data normalization and variation analysis were conducted using the “limma” R package. Heatmap was generated under p < 0.05 using the “pheatmap” R package. All data were collected from public databases; hence, ethics approval and informed consent were not required. The overview of the workflow steps is shown in Supplementary Figure S1.

To study the molecular functions (MFs), biological processes (BPs), and cellular components (CCs) of apoptosis-related genes, “clusterProfiler,” “org.Hs.e.g.,.db,” and “enrichplot” R packages were used for Gene Ontology (GO) annotation. KEGG pathway enrichment analysis was also conducted. The threshold for each analysis was set at p < 0.05.

We applied least absolute shrinkage and selection operator (LASSO) regression and support vector machine-recursive feature elimination (SVM-RFE) model for disease-specific gene selection. The optimal λ selected returned the minimum cross-validation error under 10-fold cross-validation in the LASSO model. The SVM-RFE model was executed using the “SVM” R package. Each feature’s score was sorted after removing the minimum value and iterations were sustainably performed until the best features were selected. Accordingly, differentially expressed genes (DEGs) for sepsis were identified by overlapping both results. Moreover, the diagnostic ability of DEGs was assessed using receiver operating characteristic (ROC) curves with the “pROC” R package. The volcano plot and heatmap presented the detailed expression of DEGs with “ggplot” and “pheatmap” R package.

The univariate Cox and LASSO regressions with Candidate DEGs were conducted using the “glmnet” R package to construct a prognostic signature. Then, the “timeROC” R package was employed to perform 1-, 2- and 3-week ROC analyses, and results were quantified by the area under curve (AUC). Meanwhile, the Kaplan-Meier survival analysis was conducted using ‘‘survival” and ‘‘survminer” R packages. Using the formula Risk score = (each gene’s expression×corresponding coefficient), septic patients were classified into the high-risk group (risk score > median risk score) and low-risk group (risk score ≤ median risk score). The entire dataset was divided into train and test set, keeping a ratio of about 1:1. These analyses were carried out on the whole dataset, train dataset, and test dataset, respectively.

The nomogram was built using the multivariate Cox regression with the “survival” and “rms” R packages. The concordance index (C-index) was applied to assess the discriminatory performance of hub genes. The calibration curve was applied for nomogram calibration. This algorithm can be an efficient graphical description in which each risk factor responds to the risk of mortality at 1-, 2-, and 3- week for an individual patient.

To study transcriptome differences and personalized analysis, GSEA was carried out using the “GSEABase” R package. According to the log₂ [Fold Change (FC)] value of the differential analysis, genes were ranked from high to bottom to be defined as the tested gene set. The KEGG pathway set was applied to conduct pathway-level analysis of each hub gene using the “GSVA” R package, in which the relative pathway activity represented as the t-value. Up- and downregulated pathways were determined using the “limma” R package.

CIBERSORT is a classic method to quantify the relative abundance of various cell composition from enormous gene expression profile (Chen et al., 2018). By means of CIBERSORT in the GSE65682, the proportion of 22 types of infiltrating immune cells was analyzed between sepsis and control groups, further, the correlation degree between immune cells and hub genes was calculated using Pearson correlation. The results were shown in visual diagrams via “ggplot” R language.

Based on GEO database, GSE167363 raw data was downloaded. Single cells were extracted under the criteria: nFeature_RNA>500, percent. mt<20%, percent. HB < 1, nCount_RNA>1000. The doublet and non-viable cells were removed by “DoubletFinder” R language and Scrublet algorithm. Normalization, scaling, clustering of cells was implemented by “Seurat” R package, which generated dominant immune cell types. The results were distinguished using t-distributed stochastic neighbor embedding (tSNE). The average gene expression levels in each identified cell type were estimated using “Seurat” R package. The Variation analyses were conducted using the “limma” R package.

For high quality evidence in external validation, we not only employed GSE13904, GSE26378 form GEO database, also harvested 158 adult samples from ZhongNan hospital, Wuhan University, composed of 15 healthy volunteers and 143 septic patients. All patients were adequately informed, and the study received approval from the Ethics Committee of ZhongNan Hospital, Wuhan University (2017004). The demographic and clinical characteristics of recruited samples were summarized in Supplementary Table S1. The peripheral blood mononuclear cells (PBMCs) were collected from patients on the first day after the diagnosis of sepsis using leucocyte cell separation medium kit (TBD Science, Tianjin, China), further, the total RNA was extracted. Based on thorough evaluation of its integrity and quality, we selected RNA with the RNA Integrity Number (RIN)≥7.0 and 28S/18S ratio>1.0 for library preparation and mRNA transcriptome sequencing. Subsequently, the library was constructed through the TruSeq Stranded mRNA Library Prep Kit (Illumina), and sequencing was performed using NovaSeq 6000 (Illumina). Regarding the external datasets validation, following the calculation of expression level of hub genes using RSEM (Li and Dewey, 2011) in two separate groups, differential expression analysis was performed using the “DESeq2” R package with p < 0.05, and the diagnostic value of hub genes was presented with the fitting ROC curves. Plus, the volcano plot and heatmap revealed the intuitive expression of hub genes. Besides, we explored the prognostic significance of hub genes shown in 1-, 2- and 3-week ROC analyses of Zhongnan Hospital dataset.

C57BL/6 mice (males and 8 weeks old) were housed under a standard raising condition. Ten mice were randomly divided into control and sepsis groups with five in each group. The sepsis model was established using the cecal ligation and puncture (CLP) method according to guidelines (Rittirsch et al., 2008). After anesthetizing mice with pentobarbital sodium solution (0.01 mg/g) by intraperitoneal injection, the caecum was exposed and ligatured during laparotomy. The major step was the caecal wall puncture with a 20-gauge needle and gently squeezing the stool (length: 2 mm). The surgery ended with abdomen closure layer by layer. The heating pad and normal saline injection were used for better post-operation resuscitation. In the control group, the abdomen was only opened and closed without extra operation. After 12 h, blood samples were collected and centrifuged to obtain serum. Subsequently, the mice were humanely euthanized with pentobarbital sodium (0.1 mg/g), then, the heart, lung, liver, and kidney tissues were harvested for further experiments. The Animal Care and Use Committee of the Zhongnan Hospital of Wuhan University approved the animal experiment.

Following the extraction of mice serum in control and CLP groups, the creatinine (Cr) level was determined by Creatinine Assay Kit (Sigma-Aldrich, Shanghai, China) and blood urea nitrogen (BUN) level were determined by BUN Quick Test Strips (Sigma-Aldrich, Shanghai, China) according to the manufacturer’s instructions.

Following the samples harvested from mice models, total RNA was extracted using the FastPure Plant Total RNA Isolation Kit (Vazyme, Nanjing, China). Then, cDNA was synthesized using the Hifair^® first Strand cDNA Synthesis SuperMix (Yeasen Biotechnology, Shanghai, China) with the following procedure: 25°C/5 min, 42°C/30 min, and 85°C/5 min. Next, we performed qPCR using ChamQ SYBR qPCR Master Mix (Vazyme Biotech, Nanjing, China). The amplification process was implemented in the Roche LightCycler^®96 Real-Time PCR Detection System (Hoffmann-La Roche Ltd., Shanghai, China). The primer pairs used were presented in Table 2. With GAPDH serving as reference gene, the expression levels of hub genes in serum and target organs (heart, lung, liver and kidney) were analyzed using the 2^−ΔΔCT method (Livak and Schmittgen, 2001). In general, GAPDH had an identical presentation despite distinct treatments in two groups, the relative mRNA abundance of each gene was estimated as the fold change of its expression level in CLP groups over that in control groups.

As described above, the heart, lung, liver, and kidney tissues were acquired and made into paraffin tissue slides. Sections were dewaxed through xylene, ethanol, and distilled water. Sodium Citrate Antigen Retrieval Solution (Solarbio Science and Technology, Beijing, China) was used to recover antigens. After 30-min blocking of Bovine Serum Albumin (Solarbio Science and Technology, Beijing, China), primary antibodies were added to sections and incubated overnight at 4°C. The source of primary antibodies was listed as following: anti-BCL2 antibody: Proteintech, Cat No. 68103 at 1/800 dilution; anti-FASLG antibody: Proteintech, Cat No. 60196 at 1/200 dilution; anti-IRF9 antibody: Proteintech, Cat No. 14167 at 1/500 dilution; anti-JAK3 antibody: Proteintech, Cat No. 80331 at 1/500 dilution. After that step, the sections were washed three times in Phosphate Buffered Solution (PBS) and incubated with secondary antibody at room temperature for 20 min. The final chromogenic reaction was conducted using the DAB Substrate kit (Solarbio Science and Technology, Beijing, China), and the sections were dehydrated and sealed. The target proteins were stained in brown and further histological analysis was completed by an experienced pathologist using a light microscope in a blinded method. The percentage of positive expression of target proteins was independently assessed in four kinds of organs. All the images were shown at 30* magnification.

Apoptosis of heart, lung, liver and kidney were evaluated through the terminal deoxynucleotidyl transferase-mediated fluorescein-dUTP nick-end labeling (TUNEL) technique using the Click-iT™ Plus TUNEL Assay Kit (ThermoFisher Scientific, Shanghai, China). The results were observed under Nikon Eclipse 50i Fluorescence Upright Microscope (Tokyo, Japan).

Statistical analyses were performed using R software. Comparisons between two groups were assessed with the Student’s t-test. Unless stated otherwise, statistical tests were bilateral, with p < 0.05 as significant. The Kaplan-Meier method was applied to compare the overall survival (OS) time. The RT-qPCR results were interpreted using GraphPad Prism 8. The IHC and TUNEL images were rendered with Adobe Photoshop 2022.

With access to the apoptosis-related genes, 82 genes were differently expressed in sepsis in the GSE65682 cohort (Figure 1A). Then, we conducted GO functional and KEGG pathway enrichment analysis to screen predominant signaling pathways. The top 10 significantly enriched terms in BPs, CCs, and MFs are presented in Figure 1B, and the top 30 KEGG pathways in Figure 1C. The preliminary results suggested that apoptosis-related pathways were enriched and upregulated in sepsis, and apoptosis-related genes were linked to apoptosis pathway, necroptosis pathway, NOD-like receptor signaling pathway, JAK-STAT signaling pathway, NF-kappa B signaling pathway, et al., reflecting the strong relationship between apoptosis and sepsis.

We performed two machine-learning algorithms: LASSO and SVM-RFE. The LASSO regression showed 18 apoptosis-related features (Figures 2A, B), the SVM-RFE algorithm filtered 19 genes as the optimal combination of feature genes (maximal accuracy = 0.988, minimal root mean square error = 0.0125) (Figures 2C, D). After intersecting both results, 11 DEGs were selected: CASP8, VDAC2, CHMP1A, CHMP5, FASLG, IFNAR1, JAK1, JAK3, STAT4, IRF9, and BCL2 (Figure 2E). Then, we established ROC curves indicated the AUC of nine genes were >0.85 (Figure 2F). The volcano plot and heatmap indicated that JAK3, CHMP5 and IFNAR1 were upregulated while CASP8, VDAC2, FASLG, JAK1, STAT4 and BCL2 were downregulated in sepsis (Figures 2G, H). These results showed the superiority of filtered DEGs in distinguishing sepsis from healthy cases.

To broaden the clinical significance of DEGs, we integrated survival data from the GSE65682 dataset, then we employed the univariate COX regression (Figure 3A) and LASSO regression (Figures 3B, C) to establish a survival model, ultimately, 4 hub genes: BCL2, FASLG, IRF9 and JAK3 were screened out. The AUC for 1, 2, and 3 weeks was calculated in the time-ROC curves (Figures 3D–F). Then, septic patients were stratified into low- and high-risk groups using the median cut-off value based on risk scores and the survival comparison showed that low-risk patients have greater survival time than high-risk ones (Figures 3G–I). The risk survival status charts and expression patterns of four hub genes were presented in Supplementary Figure S2, showing the survival time and rate decreased as the risk scores increased. The algorithms for the time-ROC curve, survival analysis, and expression profiles as described above were run three times separately in the whole, train and test datasets. Additionally, a nomogram model was constructed for overall survival (OS) prediction incorporating FASLG, JAK3, IRF9, BCL2 and other predictors (gender and age) (Supplementary Figures 3A, D). These results revealed the prognostic capability of hub genes for sepsis, while acknowledging their limitations.

Besides, the GSEA-KEGG analysis indicated that hub genes were closely correlated to certain functional pathways participated in the occurrence of sepsis (Supplementary Figures 3B–E, F). According to GSVA terms (Supplementary Figure S4), we observed that hub genes were enriched for multiple immune response pathways, including antigen processing and presentation, T cell receptor signaling pathway, et al.

Previous studies have demonstrated the close connection between sepsis and immune microenvironment (Delano and Ward, 2016; Liu et al., 2022). Therefore, we assessed the concrete changes in immune microenvironment using CIBERSORT (Figure 4A). The proportions of memory B cells and naive CD4⁺ T cells were lower in sepsis than in normal samples, while monocytes, neutrophils, and macrophages were higher in sepsis samples. The results of correlations between hub genes and immune cells showed that NK cells and CD8⁺ T cells had strong positive links with BCL2 and FASLG, IRF9 was negatively related to CD4⁺ T cells, whereas JAK3 exhibited a strong positive influence on B cells and CD4⁺ T cells (Figure 4B). Hence, our findings exhibited that apoptosis-related hub genes were inseparable from the changes in immune system in sepsis.

GSE167363 was downloaded from GEO databases for normalization, scaling, clustering, and highly variable genes screening. The dimensionality-reduced clusters were showed on the 2D map produced with the t-distributed t-SNE (Figure 5A). The clusters predominantly consisted of six immune cell types: T cells, B cells, monocytes, NK cells, platelets and neutrophils, with the former four types accounting for large proportions. Based on the analysis of differentially expressed genes in distinct cell types (Figure 5B), further investigation revealed hub genes showed differential expression across distinct immune cell types. In septic patients, a decrease in the expression levels of BCL2 was observed in B cells and monocytes compared to healthy samples, the expression levels of FASLG in B cells, monocytes, and T cells were lower in septic patients than in healthy subjects. In contrast, higher expression levels of IRF9 were detected in B cells and T cells, while JAK3 expression was elevated in B cells and monocytes among septic patients, comparing to healthy volunteers (Figure 5C). Subsequently, hub genes were visually distinguished in cell clusters from healthy cases and sepsis cases using multi-element expression diagrams (Figures 5D, E). These results provided a clearer understanding of the relationships between hub genes and immune cells in sepsis.

To enhance the reliability of results, the expression levels of BCL2, FASLG, IRF9 and JAK3 were investigated in three external datasets. In the GSE13904 dataset (Figures 6A–J), the expression trends of BCL2 and FASLG showed a significant decrease in sepsis, while JAK3 exhibited elevation. However, there was no significant difference in the expression level of IRF9 between two groups. The ROC curves for BCL2, FASLG and JAK3 demonstrated good AUC values. The volcano plot and heatmap yielded similar results. The verification results of the GSE26378 dataset (Figures 6K–T) and Zhongnan Hospital dataset (Figure 7) closely aligned with the findings from the GSE13904 dataset, confirming the diagnostic advantages of hub genes for septic patients.

On the other hand, we attempted to estimate the prognostic value of individual hub gene. Whereas, the time-ROC curves showed that hub genes did not exhibit prognostic nature in either the Zhongnan Hospital dataset or the GSE65682 dataset (Supplementary Figure S5).

To start with, the serum levels of Cr and BUN of sepsis group were significantly higher than that of control group (Supplementary Figure S6). Subsequently, the serum levels of hub genes were depicted in Figures 8A, F, K, P, in which we observed elevated levels of JAK3 and decreased levels of BCL2 and FASLG in sepsis group compared to controls, providing further compelling evidence corresponding to the results from external datasets. Considering the general impact of sepsis on various organs, the expressions of hub genes were assessed in four major organs. In the heart (Figures 8B, G, L, Q), IRF9 expression was upregulated in sepsis, while BCL2 was downregulated. However, there was no intergroup difference in the expression levels of FASLG and JAK3. In the lung (Figures 8C, H, M, R), all four hub genes exhibited significant upregulation in the CLP group. In the liver (Figures 8D, I, N, S), IRF9 and JAK3 were elevated in the sepsis model, while no difference was observed in the expression of BCL2 and FASLG between two groups. In the kidney (Figures 8E, J, O, T), BCL2, IRF9, and JAK3 displayed upregulation in the CLP group. The IHC results manifested the representative positive expression of hub genes as brown staining (Figure 9A). The quantification of IHC results (Figure 9B) corroborated the differences in hub gene expression between two groups, consistent with the RT-qPCR analysis. Besides, the TUNEL assay was performed on the four organs as routine detection for apoptosis (Supplementary Figure S7), in which the apoptotic cells significantly increased in the sepsis group.