Unlocking the potential of senescence-related gene signature as a diagnostic and prognostic biomarker in sepsis: insights from meta-analyses, single-cell RNA sequencing, and in vitro experiments

Data collection and processing

We obtained 866 cellular senescence regulatory genes from the CellAge database (https://genomics.senescence.info/cells/) and have listed them in Supplementary Table 1. The training cohort, GSE65682 [14], consisted of transcriptome sequencing data from whole blood samples isolated from 42 healthy control and 760 sepsis subjects, along with their corresponding follow-up duration and survival statuses. We downloaded this dataset using the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) and performed a comprehensive query of the GEO database using “sepsis” as the keyword, based on sample size in both control and sepsis groups being at least ten in “whole blood” type samples where transcriptome sequencing matrices included genes in the predictive model. To ensure accuracy, JS and QL independently performed manual queries using these criteria. Any discrepancies were resolved through discussion involving JC. Based on the search results, GSE4607 [15] and GSE95233 [16] datasets, which included follow-up information, were adopted to verify the prognostic value of the predictive model. Additionally, another 11 sepsis-related datasets (GSE9692 [17], GSE13904 [18], GSE26378 [19], GSE26440 [20], GSE28750 [21], GSE54514 [22], GSE57065 [23], GSE67652 [24], GSE69063, GSE69528 [25], and GSE131761 [26]) were selected to confirm the diagnostic value of the model. The sva package in R was used to minimize batch effects across these cohorts where possible.

To investigate the underlying mechanisms at a higher resolution, we downloaded GSE175453 [27] from the GEO, which contained scRNA-seq data from whole blood samples collected from 5 healthy control and 4 sepsis donors. The scRNA-seq data were processed using Seurat, which included data loading, quality control, and dimension reduction. Cell type annotation was performed using the SingleR package. Supplementary Material 1 contains additional details on the processing methods and filtering thresholds used. Supplementary Table 2 provides detailed information on selected GEO datasets.

Cell culture and treatment

The human monocytic cell line THP-1 and the human immortalized B cell line IM-9 were obtained from American Type Culture Collection (USA) and were maintained in RPMI-1640 media (Gibco, USA). The media were supplemented with 1% penicillin-streptomycin and 10% fetal bovine serum (FBS, Gibco, USA), and the cells were incubated in a humidified atmosphere containing 5% CO₂ at 37° C. To simulate cellular senescence induced by reactive oxygen species (ROS), THP-1 and IM-9 cells were treated with 100 μM and 60 μM H₂O₂, respectively, for a period of 24 hours [28, 29].

Gene expression difference detection and protein-protein interaction (PPI) network construction

We utilized the “limma” package in R to identify senescence-related genes that exhibited significant expression differences between control and sepsis samples. Specifically, we filtered for genes with |log fold change (FC)| > 1 and a false discovery rate (FDR) of < 0.05. Afterwards, we uploaded these differentially-expressed genes to the STRING database (https://cn.string-db.org), with a confidence level set to 0.4, in order to construct a PPI network.

Feature selection and predictive model construction

We utilized the “glmnet” package to conduct LASSO regression and identify differentially-expressed genes significantly associated with the mortality of sepsis patients. Additionally, we adopted the “randomForestSRC” package to screen for hub genes associated with mortality using the random forest algorithm. To conduct the univariate Cox regression, we utilized the “survival” package, and considered P < 0.01 as significant. The genes co-determined by LASSO, random forest, and univariate Cox regression were then included in multivariate Cox regression analyses with stepwise selection to construct a prognostic model. We defined the risk score calculated by this multivariate Cox regression model as senescence-related score (SRS). To compute SRS, we used the following formula: $\text{SRS}={\sum }_{i=1}^{n}Coefi{f}_i\ast Gen{e}_i$ where “coeff” represents the coefficient of the gene in the multivariate Cox regression model. More details and parameters regarding these algorithms can be found in Supplementary Material 1.

Unsupervised clustering

Consensus clustering was performed using the “ConsensusClusterPlus” package in R. The optimal cluster number was determined by analyzing the cumulative distribution function (CDF) curve. Next, principal component analysis (PCA) was carried out to confirm the reliability of the clusters using R’s “prcomp” function.

Functional annotation and gene set enrichment analysis (GSEA)

The genes were functionally annotated using the Metascape database (https://metascape.org/gp/index.html#/main/step1). GSEA was conducted with version 4.3.2 of the GSEA software, which was downloaded from its official website (https://www.gsea-msigdb.org/gsea/index.jsp). The Hallmark and Reactome gene sets, acquired from the Molecular Signatures Database (https://www.gsea-msigdb.org/gsea/msigdb/index.jsp), served as reference. Significance was determined by terms with nominal P < 0.05 and FDR < 0.25.

Meta-analyses

To enhance the clarity of diagnosis and prognosis values, we utilized the meta package in R to combine odds ratios (ORs) and hazard ratios (HRs) obtained from several studies. The data from each study were extracted and pooled using either fixed-effect or random-effect models depending on the degree of heterogeneity determined by I-squared (I^2) statistics.

Kaplan-Meier survival analyses

Kaplan-Meier survival analyses with log-rank tests were performed using the survival package in R. The cut-off value used to divide the GSE65682, GSE95233, and GSE4607 cohorts, which were employed to clarify the prognosis value of SRS in this study into low- and high-SRS subgroups was determined by the median SRS in the training cohort (GSE65682).

Real-time quantitative PCR (RT-qPCR)

The TRIzol reagent (Invitrogen, USA) was utilized to perform total RNA isolation of the cell samples in accordance with the manufacturer’s protocol. Subsequently, cDNA synthesis was conducted using the PrimeScript RT Reagent Kit (Takara, China), and RT-qPCR experiments were carried out through the use of the Fast SYBR Green Master Mix kit (Takara, China) on the Lightcycler 480 II system (Roche, USA). GAPDH was selected as the internal reference gene to normalize gene expressions. The primer sequence can be found in Supplementary Table 3.

Statistical analyses

The statistical analyses for the entire study were performed using R software (version 4.2.0), which can be obtained from the official website (https://cran.r-project.org/). The R code used in this study can be found in Supplementary Material 1. Unless otherwise stated, Wilcoxon signed-rank tests were used to compare continuous variables, while student’s t-tests were employed for data obtained from RT-qPCR experiments. Fisher’s exact tests and Pearson Chi-square tests were conducted to evaluate differences in categorical variables across different groups. The pROC package was utilized for receiver operating characteristic (ROC) analyses and calculation of areas under curve (AUCs). A P-value of <0.05 was considered significant in this study. The significance levels were denoted as *P < 0.05; **P < 0.01; ***P < 0.001.

Availability of data and materials

The raw data can be downloaded from the Gene Expression Omnibus (https://ncbi.nlm.nih.gov/geo/), and the R code used in this study is displayed in Supplementary Material 1.

The senescence-related genes exhibiting expression difference

The workflow of the present study is illustrated in Figure 1. Firstly, 866 senescence-related genes were collected from CellAge and the GSE65682 dataset was chosen as the training cohort. A total of 80 senescence-related genes showed differential expression between control and sepsis samples (Supplementary Table 4). The heatmap (Figure 2A) and volcano plot (Figure 2B) display the expression levels and the corresponding the fold change (log2) and statistical significance (-log10 FDR) of 80 senescence-related genes, respectively. Subsequently, a protein-protein interaction network was constructed to reveal the underlying interactions among these genes (Figure 2C). Functional annotation demonstrated that these genes were primarily associated with cellular senescence, cell proliferation, cell cycle regulation, DNA damage response and cell apoptosis. Furthermore, some immune-related pathways such as TGF-beta signaling and interleukins signaling pathways were also enriched (Figure 2D), suggesting the potential functions of these genes in sepsis pathogenesis.

Figure 1. The workflow of this study.

Figure 2. Differential expression of senescence-related genes in sepsis vs control samples. (A) Heatmap showing the level of expression of 80 senescence-related genes showing expression difference in sepsis and control samples. (B) Volcano plot displaying the fold change (log2) and statistical significance (-log10 adjusted p-value) for each gene. The red dots indicate up-regulated genes, while the blue dots indicate down-regulated genes. (C) A PPI network of the 80 differentially expressed genes associated with senescence. (D) Functional annotation of the 80 differentially expressed genes. Abbreviations: PPI, protein-protein interaction.

TGFBI and MAD1L1 were significantly associated with the mortality of sepsis

LASSO regression identified 15 out of 80 differentially expressed genes as significant predictors for sepsis mortality (Figure 3A and Supplementary Table 5). In addition, four genes, namely ABI3, TGFBI, MAD1L1, and WIPI1, were identified by random forest analyses (Figure 3B). 14 out of the total 80 genes were selected through univariate Cox regression with P < 0.01 filtering (Supplementary Table 6). Ultimately, ABI3, TGFBI, and MAD1L1 were co-determined by LASSO analysis, random forest analysis, and univariate Cox analysis (Figure 3C). By using multivariate Cox regression with stepwise function, TGFBI and MAD1L1 were included in the predictive model (Figure 3D), from which SRS was calculated as follows: SRS = -0.791*expression (TGFBI) – 1.036*expression (MAD1L1). The high expression levels of TGFBI and MAD1L1 indicated favorable prognoses in the training cohort (both P < 0.001, Figure 3E). This conclusion was then validated through meta-analyses across the training cohort as well as GSE4607 and GSE95233 cohorts (Figure 3F, 3G; TGFBI, HR = 0.33, 95% confidence interval [CI] = 0.24-0.46; MAD1L1, HR = 0.40, 95% CI = 0.28-0.56).

Figure 3. Identification of TGFBI and MAD1L1 as significant predictors of sepsis mortality. (A) LASSO regression identified 15 out of 80 genes as significant predictors of sepsis mortality. (B) Random forest analysis identified four genes, including ABI3, TGFBI, MAD1L1, and WIPI1, as significant predictors of sepsis mortality. (C) LASSO, random forest, and univariate Cox analyses identified ABI3, TGFBI, and MAD1L1 as co-determined predictors of sepsis mortality. (D) Multivariate Cox regression with stepwise selection ultimately included TGFBI and MAD1L1 in the predictive model for sepsis mortality. (E) High expression levels of TGFBI (up) and low expression levels of MAD1L1 (down) were associated with favorable prognoses in the training cohort. (F, G) Meta-analyses indicated the prognostic value of TGFBI (F) and MAD1L1 (G) in predicting sepsis mortality. Abbreviations: LASSO, Least Absolute Shrinkage and Selection Operator; TGFBI, transforming growth factor-beta induced protein; MAD1L1, mitotic spindle assembly checkpoint protein.

Unsupervised clustering based on the expressions of TGFBI and MAD1L1

The 802 samples in the GSE65682 cohort were divided into two clusters, Cluster 1 (C1) and Cluster 2 (C2), according to the TGFBI and MAD1L1-based consensus clustering method (Figure 4A and Supplementary Table 7). The clustering results were confirmed by PCA (Figure 4B). The expression levels of TGFBI (P < 0.001) and MAD1L1 (P < 0.001) were significantly downregulated in the C1 subgroup (Figure 4C). Additionally, cases in the C2 subgroup exhibited more sepsis characteristics (P < 0.001, Figure 4D), worse clinical outcomes (P < 0.01, Figure 4E), and a higher level of cellular senescence response (Nominal P < 0.01, FDR < 0.25, Figure 4F). Similar conclusions can be drawn from the analysis of GSE4607 (Supplementary Figure 1) and GSE95233 cohorts (Supplementary Figure 2).

Figure 4. Unsupervised clustering based on TGFBI and MAD1L1 expression. (A) The consensus clustering algorithm divided 802 samples in the training cohort into two subgroups, C1 and C2. (B) PCA was performed to validate the robustness of the clustering. (C) The expression levels of TGFBI (left) and MAD1L1 (right) were compared between C1 and C2 subgroups. (D) The clustering was associated with sepsis characteristics. (E) Subjects in the C2 subgroup exhibited worse prognoses than those in C1 subgroup. (F) The clustering was associated with cellular senescence levels. (G) Signaling pathways enriched in C1 and C2 samples were identified. Abbreviations: PCA, principal component analysis; C1, cluster 1; C2, cluster 2; ***P < 0.001.

To investigate the differences in pathological mechanisms between C1 and C2 subjects, GSEA was conducted using Hallmark gene sets as reference. Some signaling pathways such as reactive oxygen species (ROS) pathway, IL6 JAK STAT3 signaling, and TNF-alpha signaling via NFkB were enriched in C2 subjects, partly accounting for the unfavorable prognosis for these cases (Figure 4G).

SRS was a reliable predictor for the mortality of sepsis

The sepsis subjects in the training, GSE95233, and GSE4607 cohorts were classified into the low- and high-SRS subgroups according to the median SRS level (0.99) in the training cohort (Figure 5A). The high level of SRS exhibited unfavorable clinical outcomes in the training (P < 0.01, Figure 5B), GSE4607 (P < 0.05, Figure 5C), and GSE95233 (P < 0.01, Figure 5D) cohorts. Additionally, more deaths can be observed in the high-SRS subjects from the training (P < 0.01, Figure 5E), GSE4607 (P < 0.05, Figure 5F), and GSE95233 (P < 0.05, Figure 5G) cohorts. Finally, meta-analyses were performed to clarify the predictive ability of SRS on mortality using both continuous and binary SRS measures. The results showed that SRS was a significant prognosis predictor, with a pooled HR of 4.27 (95% CI = 1.63-11.17) for continuous SRS (Figure 5H), and a pooled HR of 2.05 (95% CI = 1.47-2.85) for binary SRS (Figure 5I).

Figure 5. SRS was a reliable predictor of prognosis in sepsis. (A) Cases in the training, GSE95233, and GSE4607 cohorts were divided into high- and low-SRS subgroups based on the same cut-off value. (B–D) Kaplan-Meier survival analyses indicated the prognostic value of SRS for sepsis in the training (B), GSE4607 (C), and GSE95233 (D) cohorts. (E–G) SRS was associated with survival status in the training (E), GSE4607 (F), and GSE95233 (G) cohorts. (H, I) Meta-analyses were conducted to combine effect values using continuous SRS (H) and binary SRS (I). Abbreviations: SRS, senescence-related score.

The clinical association analyses indicated that SRS was associated with ICU infection status (P < 0.05, Figure 6A) and age (P < 0.01, Figure 6B) in the training and GSE95233 cohorts, respectively. Furthermore, SRS was found to be an independent predictor of mortality through univariate and multivariate Cox analyses in both the training (both P < 0.01, Figure 6C, 6D) and GSE95233 cohorts (both P < 0.05, Figure 6E, 6F), after transforming continuous variables into binary variables.

Figure 6. Association between SRS and clinical features in sepsis. (A) SRS was associated with ICU infection condition in the training cohort. (B) SRS was linked with age in the GSE95233 cohort. (C, D) SRS was an independent predictor of prognosis in both univariate (C) and multivariate (D) Cox analyses in the training cohort. (E, F) SRS was an independent predictor of prognosis in both univariate (E) and multivariate (F) Cox analyses in the GSE95233 cohort. *P < 0.05, **P < 0.01.

SRS was a robust diagnostic biomarker of sepsis

After observing the satisfying performance of SRS in predicting prognosis, we proceeded to investigate whether SRS could also function as a diagnostic biomarker for sepsis using 14 cohorts collected from the GEO database. Our meta-analyses revealed that both TGFBI (pooled OR = 0.14, 95% CI = 0.08-0.24, Figure 7A) and MAD1L1 (pooled OR = 0.11, 95% CI = 0.04-0.35, Figure 7B) had significant diagnostic value. Interestingly, their combination, SRS, showed impressive diagnostic abilities with a pooled OR of 5.88 (95% CI = 3.54-9.77), as displayed in Figure 7C. Supplementary Table 8 presents the AUCs that indicate the diagnostic value of TGFB1, MAD1L1, and SRS in these selected cohorts.

Figure 7. Meta-analyses revealing the diagnostic ability of TGFBI (A), MAD1L1 (B), and SRS (C) in sepsis.

SRS was associated with the senescence of B cell and monocyte in sepsis

We utilized a total of 27,808 cell samples from five healthy control subjects and 21,644 cell samples from four sepsis patients to investigate the mechanisms underlying SRS at a higher resolution. The cell samples obtained from healthy control subjects were segregated into seven main cell types, namely HSC, monocyte, platelets, T cells, NK cells, BM, and B cells (Figure 8A). The cell samples collected from sepsis patients were grouped into HSC, monocyte, neutrophils, platelets, T cells, NK cells, and B cells (Figure 8B). We analyzed the expression levels of TGFBI and MAD1L1 in both the control and sepsis samples (Figure 8C, 8D). Our analysis revealed that TGFBI was predominantly expressed in monocytes; however, there was no significant difference in the expression levels between monocytes obtained from control and sepsis subjects (P > 0.05, Figure 8E). On the other hand, MAD1L1 was found to be highly expressed in B cells as well as NK and T cells. However, only B cells exhibited significant differences in MAD1L1 expression between control and sepsis subjects (P < 0.05, Figure 8F), while no significant differences were observed for MAD1L1 expression between NK cells (P > 0.05, Figure 8G) or T cells (P > 0.05, Figure 8H). Consequently, we investigated the relationship of TFGBI with monocyte senescence while exploring the association of MAD1L1 with B cell senescence.

Figure 8. TGFBI and MAD1L1 were associated with monocyte and B cell senescence, respectively. (A, B) Clustering and cell type annotation of cells isolated from healthy control (A) and sepsis (B) subjects. (C, D) Levels of TGFBI and MAD1L in different cells isolated from control (C) and sepsis (D) subjects. (E–H) Levels of TGFBI in monocytes from control and sepsis subjects (E), and levels of MAD1L1 in B cells (F), NK cells (G), and T cells (H) from control and sepsis subjects. (I, J) Positive association between TGFBI with oxidative stress-induced senescence in monocytes (I), and MAD1L with oxidative stress-induced senescence in B cells (J). (K, L) Up-regulation of TGFBI in THP-1 cells treated with H2O2 (K), and up-regulation of MAD1L1 in IM-9 cells treated with H2O2 (L). Abbreviations: scRNA-seq: single-cell RNA sequencing; *P < 0.05; **P < 0.01; NS: not significant.

We assumed that TGFBI and MAD1L1 were linked to oxidative stress-induced cellular senescence in monocytes and B cells, respectively. Subsequently, we conducted GSEA after dividing the monocytes or B cells from sepsis samples into low- and high-expression subgroups based on the median expression level of TGFBI or MAD1L1. Our analysis revealed a positive association between TGFBI and oxidative stress-induced senescence in monocytes (Nominal P < 0.001, FDR < 0.001, Figure 8I). Additionally, we found a similar positive association between MAD1L1 and B cells (Nominal P < 0.05, FDR < 0.05, Figure 8J). We then proceeded to validate these findings experimentally by studying THP-1 and IM-9 cells that were treated with H₂O₂. Our observations demonstrated an upregulation of both TGFBI (P < 0.01, Figure 8K) and MAD1L1 (P < 0.05, Figure 8L) expressions in these treated cells.

Figure 9A depicts the distribution of SRS in cell samples isolated from sepsis patients. Additionally, Figure 9B displays the levels of SRS in each cell type. SRS was significantly associated with cellular senescence and oxidative stress-induced senescence in B cells (both Nominal P < 0.05, both FDR < 0.05, Figure 9C, 9D) and monocytes (both Nominal P < 0.01, both FDR < 0.01, Figure 9E, 9F). Furthermore, molecular mechanisms were investigated through GSEA in B cells (Figure 9G) and monocytes (Figure 9H). Several pathways such as MYC target V1 were found to be significantly enriched, indicating the underlying biological processes linked with SRS.

Figure 9. SRS was associated with cellular senescence in B cells and monocytes. (A) Distribution of SRS levels in cell samples isolated from sepsis subjects. (B) Levels of SRS in each cell type. (C, D) Negative association between SRS and cellular senescence (C) and oxidative stress-induced senescence (D) in B cells. (E, F) Negative association between SRS and cellular senescence (E) and oxidative stress-induced senescence (F) in monocytes. (G, H) Gene set enrichment analyses showing signaling pathways associated with SRS in B cells (G) and monocytes (H).