Based on the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/), three transcriptomic datasets (GSE54514, GSE65682, and GSE95233) including normal control (NC) and sepsis (SP) samples were retrieved and downloaded. Gene annotation was performed using the corresponding platform annotation files with Perl scripts. To correct for batch effects and normalize expression data across datasets, the “sva” package was employed within the R environment. Following the exclusion of sepsis samples lacking prognostic data, 100 NC samples and 479 SP samples remained for subsequent analyses.

To identify sepsis clinical phenotypes associated with gene modules, the “WGCNA” R package was utilized to build a weighted gene co-expression network. First, we selected the 5000 most significantly ranked genes with the highest variance across all samples from the input gene expression matrix to enhance computational efficiency. Next, we chose a soft-threshold power to ensure the scale-free topology of the network. By following the scale-free topology criterion, an appropriate soft threshold was determined by calculating the network topology fit index (R²) for different power values, selecting the minimum power where R² $>$0.85. The adjacency matrix was computed based on gene expression correlations and subsequently converted into a Topological Overlap Matrix (TOM) to quantify the similarity between genes. Using the dynamic tree cut method, genes were grouped into modules with a minimum module size of 30. Module eigengenes (MEs) were utilized to assess the correlation between each gene module and phenotypic traits, identifying modules highly correlated with the target phenotype. To further refine candidate genes, candidate genes within key modules were identified according to their gene significance and module membership values, defining candidate genes as those with high gene significance (GS) and module membership (MM).

Using the GeneCard database (https://www.genecards.org/), we identified a series of MOSRG for further exploratory analysis (Supplementary Table 1). A differential expression threshold was set at $|$fold change$|$ $\ge$2 and p.adjust 0.05. With the “limma” and “heatmap” R script, differential gene expressions between the NC group and SP group were analyzed and visualized. By applying the “venn” R script, we performed an analysis of genes exhibiting differential expressions, key gene modules, and MOSRG signatures, identifying the intersecting genes as key differentially expressed MOSRG signatures. To explore potential molecular regulatory mechanisms, we used the “clusterProfiler” R script to conduct Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses on the differentially expressed MOSRG.

In addition, using the STRING database, we predicted protein-protein interactions (PPI) between different genes (https://string-db.org/).

With the “survival” R package univariate Cox regression analysis was made, based on the expression patterns of the MOSRG signature and the clinical prognosis data of sepsis samples to identify MOSRG signatures associated with sepsis prognosis. A Least Absolute Shrinkage and Selection Operator (LASSO) function model was constructed using the “glmnet” R script to determine prognostic MOSRG signature variables, followed by multivariate Cox analysis to calculate the risk coefficients for each signature variable. Based on the prognostic variable risk coefficients and expression levels, we built the MOSRG score for each SP and developed the MOSRG scoring model. The MOSRG score formula is as follows: MOSRG score = CX3C chemokine receptor 1 (CX3CR1) $\times$ –2.056 + epoxide hydrolase 2 (EPHX2) $\times$ –1.720 + RNA-binding protein, ribonuclease 2 (RNASE2) $\times$ 2.065. By the optimal cutoff value for the 28-day clinical prognosis outcomes, we classified the samples into low- and high-MOSRG score subgroups and plotted a scatter plot of the subgroup classifications using the “ggplot2” R script. The “survival” and “survivalROC” R scripts were employed to generate 28-day clinical survival curves, as well as time-dependent receiver operating characteristic (ROC) curves for 7-day, 14-day, and 28-day outcomes. Using the “caret” machine learning algorithm, sepsis samples were split into two independent cohorts, using a 7:3 ratio for the training and validation groups. Using the risk coefficients and expression profiles from the MOSRG prognostic signature, we computed MOSRG scores for both cohorts. Prognosis and time-dependent ROC curves were then generated for both cohorts.

Using clinical baseline data of SP and the MOSRG scoring signature, we performed both multivariate and univariate Cox analyses by “survminer” R script to assess the independent prognostic significance in predicting the clinical outcomes. By conducting “rms” R script, we constructed a diagnostic nomogram incorporating clinical-pathological characteristics and the MOSRG score to estimate survival probabilities of the samples. Calibration curves and the concordance index (C-index) were created using the “regplot” and “pec” R scripts to assess the reliability of the nomogram diagnostic model and the accuracy of each variable. Using the “ggDCA” R script, we plotted the decision curve analysis (DCA) for different clinical-pathological characteristics, the nomogram model, and the MOSRG score.

Using the “ConsensusClusterPlus” R script, we classified the sepsis samples based on the “K-means” algorithm with a classification range of k = 2 to 9. The optimal classification was determined based on the delta area curve and consensus cumulative distribution function (CDF) curve results. Clinical survival curves for the 28-day outcomes of the molecular subgroups of SP were plotted using the “survival” R script. We assessed the immune score and Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data (ESTIMATE) score using the “estimate” R package. To quantify immune infiltration status, single-sample Gene Set Enrichment Analysis (ssGSEA) analysis was conducted using the “Gene Set Variation Analysis (GSVA)” R package, with marker genes for 23 immune cell types. Additionally, using the reference gene list for KEGG signaling pathways, we utilized the “GSVA” R script to calculate the differential regulated signaling pathways.

The single-cell RNA sequencing (scRNA-seq) dataset GSE167363, comprising data from sepsis patients, was retrieved from the GEO database for this study [15] (https://www.ncbi.nlm.nih.gov/). Following the exclusion of samples that did not meet inclusion criteria, peripheral blood mononuclear cells (PBMCs) from 2 healthy controls (HCs) and 5 SPs were selected for downstream analysis. First, the raw data underwent quality control processing using the “Seurat” package in R (v4.0.5, Satija Lab, New York, NY, USA), filtering cells with gene counts between 200 and 5000 and mitochondrial gene content below 10%, while removing low-expressed genes. The dataset was normalized via the “LogNormalize” method, followed by selection of the 2000 most variable genes. To control for batch effects across different experimental conditions, we applied the FindIntegrationAnchors and IntegrateData functions to integrate the samples. After data normalization, principal component analysis (PCA) was performed to reduce data dimensionality, selecting the top 20 principal components for cell clustering. A K-nearest neighbor (KNN) graph was constructed using the FindNeighbors and FindClusters functions for clustering, and cell population distribution was visualized using t-distributed Stochastic Neighbor Embedding (t-SNE) and Uniform Manifold Approximation and Projection (UMAP) algorithms. Cell types were annotated using the SingleR algorithm and known marker genes, leading to the cell subpopulation identification such as Natural Killer (NK) cells, platelets, monocytes, and neutrophils. Finally, we further analyzed and displayed their expression distribution of prognostic-related genes using VlnPlot and FeaturePlot.

The cell line RAW264.7 was validated through short tandem repeat (STR) profiling and confirmed to be free of mycoplasma contamination. RAW264.7 murine macrophages (obtained from the Shanghai Cell Bank, Chinese Academy of Sciences, Shanghai, China) were cultured in Dulbecco’s Modified Eagle Medium (DMEM; Cat# C11995500BT, Beyotime Biotechnology, Shanghai, China) supplemented with 10% heat-inactivated fetal bovine serum (FBS, Cat# 164210-50, Tianhang Bio, Shaoxing, China), 100 U/mL penicillin and 100 µg/mL streptomycin (Cat# C0222, Beyotime Biotechnology, Shanghai, China). Cells were incubated at 37 °C in a humidified atmosphere containing 5% CO₂, and the medium was refreshed every 2–3 days to maintain optimal growth conditions.

RAW264.7 cells were cultured in DMEM and seeded in 6-well plates at a density of 2 $\times$ 10⁵ cells per well until reaching approximately 50–60% confluency. siRNA or negative control siRNA (siNC) was diluted in serum-free Opti-MEM medium (Cat# 31985-070, Gibco, New York, USA) to a final concentration of 50 nM and mixed with Lipofectamine 3000 transfection reagent (Cat# L3000-015, Invitrogen, Carlsbad, USA) according to its instructions. The transfection mixture was gently added to each well and incubated for 48 hours. Subsequently, the siRNA-transfected cells were stimulated with 500 ng/mL lipopolysaccharide (LPS, Cat# L2630, Sigma-Aldrich, St. Louis, USA) for an additional 12 hours. The target sequence was as follow: RNASE2 siRNA, 5^′-AGAACAUGUAAGGGCUUAAAUTT-3^′. siRNA for RNASE2 was generated by Sangon Biotech (Shanghai, China).

Cells were seeded in 6-well plates at a density of 1 $\times$ 10⁵ cells per well and incubated overnight to allow for adherence. Subsequently, the culture medium was discarded, and each well was incubated with 1 mL of serum-free medium containing 10 µM 2^′,7^′-dichlorodihydrofluorescein diacetate (DCFH-DA; Cat# S0033S, Beyotime, Shanghai, China) to allow for intracellular ROS detection. Cells were incubated in the dark for half an hour. Following incubation, cells were washed three times with PBS, each for 5 minutes, to remove unreacted probes. Fluorescence microscopy (Olympus IX73, Tokyo, Japan) was used for detection with an excitation and emission wavelength of 488 nm and 525 nm.

Cells were plated in a 6-well plate at a density of 1 $\times$ 10⁵ cells per well and incubated overnight until adherent. The medium was then discarded, and 1 mL serum-free medium containing 5 µg/mL JC-1 (5,5^′,6,6^′-Tetrachloro-1,1^′,3,3^′-tetraethylbenzimidazolylcarbocyanine iodide, Cat# C2006, Beyotime, Shanghai, China) was added, incubating in the dark for 20 minutes. After two washes with pre-cooled PBS, fluorescence microscopy was used for detection. For JC-1 monomers and aggregates, the excitation wavelength and emission wavelength were 485, 525 nm and 530, 590 nm respectively.

Cells were harvested, homogenized with Radio-Immunoprecipitation Assay (RIPA) lysis buffer, centrifuged to collect the supernatant. Detection was performed following the instructions provided with the MDA kit (Cat# S0131S, Beyotime, Shanghai, China). Main steps included mixing samples with working solution, incubating at 37 °C, centrifuging, and measuring absorbance of the supernatant. Absorbance was quantified at a wavelength of 530 nm, and MDA concentration was determined.

The levels of IL-1$\beta$, IL-6, and TNF-$\alpha$ in the culture supernatants were quantified using enzyme-linked immunosorbent assay (ELISA) kits (Beyotime, Shanghai, China). First, the cell supernatant was diluted to an appropriate concentration and added to 96-well plates pre-coated with specific antibodies. Then, standards and samples were added with 1-hour incubation. After thorough washing, a secondary antibody (provided in the ELISA kit) conjugated with an enzyme was added with 30 minutes incubation. Finally, substrate solution was added and the absorbance was measured at 450 nm (BioTek, Winooski, VT, USA) after 30 minutes. The supernatant concentrations were calculated by the standard curve.

Cells were seeded at a density of 5000 cells per well in 96-well plates and incubated overnight to allow adherence. 10 µL of Cell Counting Kit-8 (CCK-8) solution (Cat# C0037, Beyotime, Shanghai, China) was added to each well followed by 2 hours incubation. Absorbance of each well was measured at 450 nm. After the blank control absorbance comparison, cell viability was calculated.

Cells were harvested and lysed using an appropriate volume of RIPA lysis buffer, and centrifuged to collect the supernatant. Protein samples were resolved using electrophoresis on a 10% sodium dodecyl sulfate - polyacrylamide gel electrophoresis (SDS-PAGE) gel, followed by transfer of the separated proteins to a polyvinylidene fluoride (PVDF) membrane (Cat# IPVH00010, Millipore, Darmstadt, Germany). After 1 hour block with 5% non-fat milk in (Tris-buffered saline containing 0.1% Tween-20 (TBS-T), Cat# P1379, Sigma-Aldrich, St. Louis, MO, USA, pH 7.4), the membrane was incubated with RNASE2 antibody (1:1000, A9949, ABclonal, Wuhan, China) and $\beta$-actin antibody (1:2000, ab8226, Abcam, Cambridge, MA, USA) at 4 °C overnight, TBS-T washed three times, and exposed to the appropriate secondary antibody (Rabbit anti-secondary antibody, 1:20,000, ab288151; Mouse anti-secondary antibody, 1:20,000, ab6728, Abcam, Cambridge, USA) at room temperature for one hour. After that, Enhanced chemiluminescence (ECL) substrate (Cat# P0018FS, Beyotime, Shanghai, China) was used for detection, and the imaging system (ChemiDoc, Bio-Rad, Hercules, CA, USA) was used for capturing and analyzing images.

Statistical analyses were conducted using Perl software (v5.30.0; The Perl Foundation, Walnut, CA, USA), R software (v4.3.1; The R Foundation for Statistical Computing, Vienna, Austria), GraphPad Prism (v9.5.0; GraphPad Software, San Diego, CA, USA) and SPSS 17.0 (SPSS Inc., Chicago, IL, USA). All experiments were independently repeated three times to ensure data accuracy and reproducibility. For comparisons between two groups, the unpaired Student’s t-test was used for normally distributed data, while the Wilcoxon rank-sum test was applied for non-normally distributed data. For comparisons among more than two groups, one-way analysis of variance (ANOVA) followed by Tukey’s multiple comparisons test was conducted. The normality test is performed using the “shapiro.test” function in the R language. Correlation analysis was performed using Pearson’s correlation coefficient to assess the relationship between MOSRG scores and immune cell infiltration levels. Multiple testing corrections were applied using the Benjamini-Hochberg method, and adjusted p-values 0.05 were considered statistically significant. Data are presented as mean $\pm$ standard deviation (SD), with significance levels indicated as *p 0.05, **p 0.01, ***p 0.001.

To uncover critical gene modules linked to sepsis, we integrated three GEO datasets and extracted data from 100 NC samples and 479 SP samples with 28-day survival prognosis information to construct the WGCNA model. According to the gene expression matrix, we performed clustering analysis using a dendrogram to identify and exclude outlier samples (Fig. 1A). With a soft-thresholding power ($\beta$) set at 11, we constructed a scale-free network and identified distinct gene modules from the transcriptomic matrix (Fig. 1B). As shown in Fig. 1C, the clustering tree results displayed the distribution of each gene module (height = 0.25). Similar gene modules were merged by algorithm “dynamic tree cut”, resulting in 14 independent gene modules (Fig. 1D). Correlation analysis was employed to evaluate potential associations between genes within each module, revealing that there were no significant correlations between gene matrices within each module (Fig. 1E). Additionally, the correlation results between gene module members demonstrated significant independence among module members (Fig. 1F). Using Pearson’s correlation algorithm, we assessed the associations between each gene module and clinical traits. The brown module was significantly positively correlated with NC traits (r = 0.61, p = 7 $\times$ 10^-61) and exhibited a strong negative correlation with SP traits (r = –0.61, p = 7 $\times$ 10^-61), suggesting it may be the module most relevant to clinical traits (Fig. 1G). Moreover, scatter plot analysis showed a strong positive correlation between module membership and gene significance within the brown module (r = 0.74, p = 4.3 $\times$ 10^-189) (Fig. 1H).

Fig. 1.

WGCNA-based identification of key gene modules correlated with clinical traits. (A) Clustering analysis of samples. (B) Determination of the soft-thresholding power ($\beta$). (C) Cluster dendrogram analysis of different gene modules. (D) Dynamic tree cutting analysis. (E) Correlation analysis of genes within different gene modules. (F) Correlation analysis of gene modules. (G) Correlation analysis between clinical traits and gene modules. (H) Scatter plot analysis of module membership versus gene significance. WGCNA, weighted gene co-expression network analysis; NC, normal control; SP, sepsis.

Using a differential threshold of $|$FC$|$ $>$2 and p.adjust 0.05, we analyzed gene expression patterns between the NC and SP groups, identifying 305 significantly differentially expressed gene signatures (Fig. 2A,B). Venn diagram analysis revealed the intersection of differentially expressed genes, WGCNA and MOSRG signatures, identifying 59 key MOSRG signatures with differential expression (Fig. 2C). Based on enrichment analysis algorithms, we evaluated the potential molecular regulatory mechanisms of the 59 DE-MOSRG in sepsis. GO enrichment analysis indicated that these DE-MOSRG signatures were associated with immune regulatory functions. KEGG enrichment analysis showed that these DE-MOSRG signatures were involved in regulating immune cell receptor signaling pathways (Fig. 2D,E).

Fig. 2.

Screening and enrichment analysis of differentially expressed MOSRG signatures. (A) Differential gene expression analysis between NC and SP groups. Filtering thresholds: $|$fold change$|$ $\ge$2 and p.adjust 0.05; red: upregulated genes, blue: downregulated genes. (B) Heatmap analysis of differentially expressed gene signatures. (C) Identification of differentially expressed MOSRG signatures. (D,E) KEGG and GO enrichment analysis. MOSRG, mitochondrial oxidative stress-related gene; KEGG, Kyoto Encyclopedia of Genes and Genomes; GO, Gene Ontology; DEG, differentially expressed gene; MHC, major histocompatibility complex; BP, biological process; CC, cellular component; MF, molecular function.

We further assessed the prognostic value of the 59 MOSRG signatures in sepsis and constructed a prognostic risk model to indicate clinical survival outcomes. By integrating the clinical data and MOSRG signature expression features from 479 SP samples, the prognostic value of each variable in SP was assessed by univariate Cox regression analysis. This analysis identified 16 MOSRG signatures associated with SP prognosis, including 15 beneficial factors and one risk factor (Fig. 3A). Further analysis using LASSO regression identified 6 MOSRG features relevant to SP prognosis (Fig. 3B). Subsequently, to determine the independent prognostic impact of these six features, multivariate Cox regression analysis was carried out. The results indicated that RNASE2, CX3CR1, and EPHX2 were independent prognostic variables for SP clinical outcomes. Through PPI network analysis, we assessed the interactions between these three independent prognostic variables (Supplementary Fig. 1A). The results of GO enrichment showed RNASE2, CX3CR1, and EPHX2 may be associated with biological functions such as chemotaxis and neuronal cell body membrane. KEGG enrichment analysis revealed that RNASE2, CX3CR1, and EPHX2 might be involved in the regulation of signaling pathways such as arachidonic acid metabolism, peroxisome, and chemokine signaling pathway (Supplementary Fig. 1B,C). By conducting the risk coefficients and expression profiles of RNASE2, CX3CR1, and EPHX2, we developed a MOSRG scoring model for risk stratification in SP and to indicate clinical survival outcomes for different subgroups (Fig. 3C,D). PCA pattern analysis showed a distinct separation in the distribution patterns between MOSRG score subgroups, highlighting significant differences between them (Fig. 3E). Clinical survival curve analysis demonstrated that the 28-day survival rate in the low MOSRG subgroup was significantly higher, suggesting that SP samples in the low MOSRG subgroup may have better clinical survival benefits (Fig. 3F). Time-dependent ROC curve analysis showed area under the curve (AUC) values of 0.703, 0.677, and 0.666 for 7, 14, and 28 days, respectively, indicating a high diagnostic predictive ability (Fig. 3G).

Fig. 3.

Construction of a prognostic risk model based on MOSRG signatures in sepsis. (A) Identification of MOSRG signatures associated with SP prognosis using univariate Cox analysis. (B) Construction of the LASSO model to determine MOSRG feature prognostic variables. (C,D) Development of the MOSRG scoring prognostic system. (E) Principal component analysis (PCA) plot analysis. (F) 28-day clinical survival curves analysis for MOSRG scoring subgroups. (G) Time-dependent ROC curve analysis. LASSO, Least Absolute Shrinkage and Selection Operator; ROC, receiver operating characteristic; AUC, area under the curve; RNASE2, RNA-binding protein, ribonuclease 2; CX3CR1, CX3C chemokine receptor 1; EPHX2, epoxide hydrolase 2.

Using the “caret” algorithm, we divided the 479 SP samples into training set and validation set at a 7:3 ratio and constructed the MOSRG scoring model to confirm the robustness and predictive performance of the MOSRG scoring system for clinical outcomes in sepsis patients. Based on RNASE2, CX3CR1 and EPHX2 risk coefficients and expression profiles, we developed the MOSRG scoring model in both cohorts and classified SP samples into two MOSRG scoring subgroups (Fig. 4A–D). In the training set, we observed that the 28-day clinical survival outcomes for the low MOSRG scoring subgroup were significantly better than those for the high MOSRG scoring subgroup (p = 0.002, Fig. 4E). Time-dependent ROC curve analysis indicated AUC values of 0.706, 0.678, and 0.663 for 7, 14, and 28 days, respectively (Fig. 4F). Notably, in the validation set, the 28-day survival outcomes for the low MOSRG scoring subgroup were also significantly superior to those for the high MOSRG scoring subgroup (p = 0.003, Fig. 4G). Time-dependent ROC curve analysis showed AUC values of 0.697, 0.675, and 0.677 for 7, 14, and 28 days, respectively (Fig. 4H). According to the survival analyses of both independent cohorts, we conclude that the MOSRG scoring system, constructed using MOSRG prognostic signatures, is capable of reliably predicting the 28-day survival outcomes.

Fig. 4.

Independence and stability validation of the MOSRG prognostic scoring model. (A,B) Expression analysis of MOSRG prognostic signatures. (C,D) Construction of the MOSRG prognostic scoring model in the training and validation cohorts. (E,F) Analysis of 28-day clinical survival curves and time-dependent ROC curves for sepsis samples in the training cohort. (G,H) Analysis of 28-day clinical survival curves and time-dependent ROC curves for sepsis samples in the validation cohort.

By integrating the clinical and pathological characteristics of SP samples with the MOSRG scoring system, we evaluated the independent prognostic value. Univariate and multivariate cox regression indicated age and MOSRG score were associated with poor prognosis of SP (Fig. 5A,B). ROC curve analysis demonstrated that the AUC for the MOSRG score (0.714) was significantly higher than that for gender (AUC = 0.485) and age (AUC = 0.570), highlighting the superior prognostic capability of the MOSRG scoring system (Fig. 5C). Integrating the clinical and pathological characteristics with the MOSRG score, we constructed a nomogram to predict 7-, 14-, and 28-day survival probabilities (Fig. 5D). Calibration analysis revealed that the survival probabilities estimated by the nomogram were well-aligned with the actual survival rates of SP samples (Fig. 5E). The C-index curve results showed the accuracy of using the MOSRG scoring system to assess SP sample outcomes was significantly superior to that of age and gender (Fig. 5F). Moreover, the DCA analysis indicated that the nomogram model exhibited significantly higher accuracy in predicting clinical survival outcomes of SP compared to the MOSRG score and clinicopathological variables (Fig. 5G).

Fig. 5.

Independent prognostic analysis and nomogram diagnostic model construction of the MOSRG scoring system. (A,B) Independent prognostic value assessment based on clinical-pathological features and MOSRG scoring. (C) ROC curve analysis for different variables. (D) Nomogram model construction. (E,F) Calibration curves and concordance index (C-index) evaluation. (G) The decision curve analysis (DCA) of different variables. **p 0.01, ***p 0.001.

Using the GSVA algorithm, we assessed potential mechanisms underlying survival outcome differences. The GSVA revealed that subgroup with low MOSRG score, pathways such as O-glycan biosynthesis and terpenoid backbone biosynthesis were significantly downregulated. A significant downregulation of multiple immune function-associated pathways, including the intestinal IgA immune network, T cell receptor (TCR) signaling, and NK cell-mediated cytotoxicity, was observed in the high MOSRG score subgroup (Fig. 6A). ESTIMATE showed subgroup with the high MOSRG score, immune scores and ESTIMATE scores were significantly lower, suggesting that SP samples in the high MOSRG score subgroup may exhibit a pronounced immunosuppressive state (Fig. 6B). Pearson’s correlation analysis assessed potential links between MOSRG scores and immune microenvironment characteristics, revealing a significant positive correlation with macrophages, mast cells, and activated dendritic cells, while showing a marked negative correlation with other immune infiltrating cells such as T follicular helper cells, activated CD8⁺ T cells, natural killer T cells, myeloid-derived suppressor cells (MDSCs), and immature B cells (Fig. 6C). The quantitative assessment of immune cell infiltration showed that in the high MOSRG score subgroup, the infiltration levels of several immune cell types, such as activated B cells, CD4⁺ T cells, CD8⁺ T cells, and neutrophils, were significantly lower (Fig. 6D). Building on these results, we explored the possible interplay between the MOSRG score and immune microenvironmental features in sepsis, suggesting that reduced immune suppression may be linked to worse outcomes in sepsis.

Fig. 6.

Immune microenvironment landscape analysis of MOSRG scoring subgroups. (A) GSVA between MOSRG scoring subgroups. (B) Assessment of immune scores and ESTIMATE scores. (C) Correlation analysis between MOSRG scoring and immune microenvironment features. (D) Quantitative analysis of 23 types of immune infiltrating cells within MOSRG scoring subgroups. GSVA, Gene Set Variation Analysis; ESTIMATE, Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data. *p 0.05; **p 0.01; ***p 0.001; ns, not significant.

We further identified the molecular subtyping landscape of MOSRG based on three MOSRGs prognostic variables (RNASE2, CX3CR1, and EPHX2). Using the optimal model parameters from consensus clustering analysis, we classified the sepsis samples into two distinct MOSRG molecular subgroups (k = 2). A total of 303 samples were assigned to Subgroup A, and 176 samples to Subgroup B (Fig. 7A–C). PCA revealed two clear separation patterns between the MOSRG molecular subgroups (Fig. 7D). The 28-day clinical survival curves indicated that samples in MOSRG Subgroup A had a longer overall survival time compared to those in Subgroup B (Fig. 7E). Additionally, we observed a significant increase in MOSRG scores in the poorer prognostic MOSRG subgroup, highlighting the possible interplay between the MOSRG scoring system and molecular subgroups (Fig. 7F). Using ESTIMATE and ssGSEA immune infiltration algorithms, we explored potential distinctions in immune cell infiltration patterns between MOSRG molecular subgroups. Moreover, significantly higher immune scores and ESTIMATE scores could be observed in the subgroup with better prognosis, suggesting a pronounced immune activation state in this subgroup (Fig. 7G,H). Quantitative results indicated a significant increase in the proportion of various immune cells in MOSRG subgroup A (Fig. 7I). GSVA further revealed potential regulatory KEGG signaling pathways between MOSRG molecular subgroups. In the poorer prognostic MOSRG Subgroup B, immune regulation-related pathways were significantly downregulated, whereas metabolic pathways such as sucrose, starch, sphingolipid and glycerophospholipid metabolism were notably upregulated (Fig. 7J). Based on these results, we found that MOSRG-based prognostic signatures can accurately assess the MOSRG molecular subtype landscape in sepsis and indicate the clinical prognosis outcomes and immune microenvironment infiltration of different molecular subgroups.

Fig. 7.

Identification and profiling of MOSRG molecular subtypes and immune landscape in sepsis. (A–C) Identification of molecular subtype characteristics based on prognostic MOSRG signatures. (D) PCA pattern plot analysis. (E) Clinical prognosis curves for 28-day survival among MOSRG molecular subtypes. (F) Comparison of MOSRG Scores Across Different MOSRG Molecular Subtypes. (G,H) Assessment of immune scores and ESTIMATE scores within MOSRG molecular subtypes. (I) Quantitative analysis of the proportion of 23 types of immune cell infiltration. (J) GSVA elucidating differential regulation of KEGG signaling pathways between MOSRG molecular subtypes. CDF, cumulative distribution function. *p 0.05; **p 0.01; ***p 0.001; ns, not significant.

We further analyzed the MOSRG prognostic signature characteristics within different cellular subpopulations of sepsis using single-cell sequencing analysis. By dataset GSE167363, we obtained scRNA-seq data from two HC and five SP. After quality control, we retained a sum of 38,562 high-quality cells for further analysis (Fig. 8A). Following normalization, we identified 2000 genes with high variability for dimensionality reduction analysis, including JCHAIN, IGKC, HAMP, SERPINB2, and IGHA1 (Fig. 8B). Using the “FindAllMarkers” function, we identified marker genes for each cell subpopulation and visualized the distribution of 21 cell subpopulations using UMAP/tSNE dimensionality reduction plots (Fig. 8C–E). The violin plot results indicated that RNASE2, CX3CR1, and EPHX2 were expressed in all 21 cell subpopulations, with RNASE2 showing the most prominent expression (Fig. 8F). Cell annotation using SingleR revealed five cell subtypes in sepsis: monocytes, nature killer cells, platelets, neutrophils, and erythroblasts (Fig. 8G). The distribution features of these five cell subpopulations were further illustrated by UMAP and tSNE dimensionality reduction plots (Fig. 8H,I). We then analyzed the expression profile of MOSRG prognostic signatures within each cell subpopulation. UMAP results showed that CX3CR1 was observed in monocytes and NK cells, while EPHX2 was expressed in platelets. Notably, RNASE2 exhibited the most significant expression in monocytes (Fig. 8J).

Fig. 8.

Single-cell sequencing analysis of MOSRG prognostic signatures. (A) Quality control and normalization from HC (n = 2) and SP (n = 5) samples. (B) Variable genes identification. (C,D) tSNE and UMAP plots illustrating the distribution of 21 cell subtypes. (E) Expression analysis of marker genes for the 21 cell subtypes. (F) Expression of MOSRG prognostic signatures across the 21 cell subtypes. (G) Single-cell subtype annotation. (H,I) tSNE and UMAP plots depicting the distribution characteristics of 5 cell types. (J) Expression distribution of CX3CR1, EPHX2, and RNASE2 across the 5 cell types.

To examine the impact of RNASE2 on mitochondrial oxidative stress, we cultured macrophage cells and treated them with LPS to establish a sepsis model, utilizing siRNA to knock down RNASE2 expression. Western blotting (WB) results showed that compared to the Normal group, LPS stimulation upregulated RNASE2 levels, while siRNA effectively interfered with RNASE2 expression (Fig. 9A,B). LPS exposure significantly suppressed RAW264.7 macrophage viability, as determined by the CCK-8 assay, whereas interfering with RNASE2 effectively improved cell viability (Fig. 9C). In addition, we further verified the efficiency of the cell model interference at the mRNA level, and the PCR results suggested that the expression of RNASE2 in RAW264.7 macrophage was significantly reduced after RNASE interference compared with the LPS group (Fig. 9D). MDA, mitochondrial membrane potential, and ROS levels reflect cellular oxidative stress. Results indicated that oxidative stress levels in cell mitochondria significantly increased after LPS stimulation, manifested by elevated MDA levels, decreased mitochondrial membrane potential, and increased ROS production. However, interfering with RNASE2 effectively inhibited mitochondrial oxidative stress, reducing MDA and ROS levels and increasing mitochondrial membrane potential (Fig. 9E–I). Furthermore, we further investigated the changes in inflammatory mediators in RAW264.7 macrophages after RNASE2 interference. Compared to the model group, our result showed lower levels of IL-1$\beta$, IL-6, and TNF-$\alpha$ in RAW264.7 macrophages following siRNASE2 interference (Fig. 9J–L). Taken together, our findings suggest that inhibiting RNASE2 expression can significantly ameliorate mitochondrial oxidative stress damage in sepsis and reduce the occurrence of inflammation, highlighting the potential role of RNASE2 in sepsis.

Fig. 9.

RNASE2 inhibition attenuates mitochondrial oxidative stress of sepsis. (A,B) Western blots of RNASE2 in RAW264.7 macrophages following treatment with LPS (n = 3). (C) CCK-8 assay of RAW264.7 macrophages (n = 3). (D) Verification of interference efficiency of RNASE2 in RAW264.7 macrophages following treatment with LPS and siRNASE2 groups. (E,F) Mitochondrial membrane potential test by JC-1 aggregates and monomers, scale bars correspond to 20 µm for low-magnification images and 10 µm for high-magnification views (n = 3). (G,H) ROS immunofluorescent staining revealed oxidative stress levels in the cells, scale bars correspond to 40 µm for low-magnification images and 20 µm for high-magnification views (n = 3). (I) The MDA concentration in RAW264.7 macrophages was measured (n = 3). (J–L) Quantitative analysis of IL-1$\beta$, IL-6, and TNF-$\alpha$ in RAW264.7 macrophages (n = 3). a: p 0.05 as compared with the Normal group, b: p 0.05 as compared with the LPS group, c: p 0.05 as compared with the LPS+siNC group. LPS, lipopolysaccharide; siNC, negative control siRNA; CCK-8, Cell Counting Kit-8; JC-1, 5,5^′,6,6^′-Tetrachloro-1,1^′,3,3^′-tetraethylbenzimidazolylcarbocyanine iodide; ROS, reactive oxygen species; MDA, malondialdehyde; IL-1$\beta$, interleukin-1 beta; IL-6, interleukin-6; TNF-$\alpha$, tumor necrosis factor-alpha.