Six-gene Signature Based on Metabolism Is Closely Related to the Prognosis and Immune Infiltration of Patients With Sepsis


 Background Sepsis is a leading cause of mortality and morbidity in the intensive care unit. Current studies indicated that metabolism associated genes (MAGs) have critical roles in sepsis T his study aim s to construct a gene signature for sepsis and explore its possible mechanism ss. Methods: Differentially expressed metabolism associated genes (DEMAGs) were identified in sepsis patients based on the Gene Expression Omnibus (GEO) database. Univariate and multivariate Cox regression analyses were performed to identify and construct the prognostic gene signature. Quantitative real time PCR was applied to examine the mRNA level of sig nature genes in the sepsis mice model established by cecalligation and puncture (CLP). Next, Gene S et Enrichment Analysis (GSEA) was performed to further understand the underlying molecular mechanisms of gene signature We then assessed the abundance of infiltrating immune cell s in each sample to explore the immune microenvironment of sepsis patient s Finally, the Tumor Immune Dysfunction and Exclusion (TIDE ) and SubMap algorithms were used to explore the immunotherapy response of sepsis patients .Results: A novel six gene signature (including ELANE , NQO2 , TLR2 , PTGDS , SMAD3 , CD3E ) was established for sepsis prognosis prediction , which could accurately predict the overall survival ( of sepsis patients. In the sepsis mice model, except PTGDS, the mRNA expression of the other five genes was consistent with that of sepsis patients. T he result s of GSEA analysis indicate that the prognostic signature w as closely related to immunity function Besides , we found that the risk score was strikingly negatively correlated with the tumor microenvironment (TME) immunecells infiltration and expression of critical immune checkpoints The patients in the low risk group were more likely to benefit from immune therapy through the TIDE and SubMap analysis Conclusion : In conclusion, our signature can predict the OS of sepsis patients and provide potential guidance for exploring patients who may benefit from immunotherapy.


Background
Sepsis is a syndrome of the immune system out of control caused by infection, which can be accompanied by multiple organ failure and systemic inflammatory response [1,2]. The incidence of sepsis is high, with more than 18 million cases of severe sepsis worldwide each year, and this number is increasing at a rate of 1.5%-8.0% per year [3,4]. The mortality rate of sepsis has surpassed that of myocardial infarction as the leading cause of death among non-cardiac patients in the intensive care unit. With advances in critical care treatment technology, the hospital mortality rate of sepsis patients is still as high as 17%, although it has decreased significantly [4][5][6].
Therefore, early diagnosis of sepsis and its effective prevention are the keys to improve patient survival.
The metabolic process is necessary to maintain life, including lipid metabolism, glucose metabolism, amino acid metabolism and so on. A growing body of research suggests that changes in metabolic processes are a hallmark of sepsis [7,8], especially the changes of glucose and lipid metabolism that play important roles in sepsis [9,10]. In general, glucose metabolism mainly affects the function of immune cells, while lipid metabolism affects the homeostasis of plasma lipoproteins, both of which have a significant influence on the progression and survival of sepsis patients. PKM2 is a key gene in glycolysis pathway, which affects the function of immune cells by regulating Warburg effect and plays an important role in sepsis [11]. Pharmacological inhibition of PKM2-related pathways has been shown to protect mice from fatal sepsis [12]. Studies have shown that the up-regulation of PCSK9, a gene involved in cholesterol metabolism, is associated with acute organ failure in sepsis [13,14]. The high expression of PCSK9 in sepsis reduces the clearance of pathogen lipids, and PCSK9 inhibition may be an effective strategy in the treatment of sepsis [15]. These genes, which are associate with metabolic processes, are key markers of sepsis progression and treatment. However, the prognostic value and potential mechanisms of metabolism associated genes (MAGs) in sepsis has not been fully explored.
In this study, we screened the DEMAGs based on the GEO database and evaluated the prognostic value of these genes. In addition, we analyzed the biological functions related to prognosis by GSEA to further analyze the immune microenvironment and immunotherapy in sepsis.
Furthermore, we further validated the expression of prognosis signature in sepsis mice model by RT-qPCR. These results provide a research direction for exploring the molecular mechanism of MAGs in immunotherapy of sepsis.

Data collection
The transcriptome array data of sepsis patients were downloaded from the GEO database， inlcuding the GSE69528 and GSE65682 datasets. The GSE69528 dataset comprised of 111 samples, including 28 normal and 83 sepsis samples. The GSE65682 dataset contained 802 samples, including 42 normal and 760 sepsis samples. After removing 292 patients from the GSE65682 dataset due to a lack of survival data, a total of 468 sepsis patients and their information were utilized to identify prognostic DEMAGs, which was used to built a prognostic risk model. We downloaded the gene set (h.all.v7.2.entrez.gmt) from the GSEA website (https://www.gsea-msigdb.org/gsea/index.jsp), and a total of 1811 MAGs were obtained from metabolism-associated pathways by enrichment analysis. The detailed MAGs were listed in Supplementary materials Table S1.

Identification of differentially expressed genes (DEGs)
The DEGs were determined between normal and sepsis samples in the GSE65682 and GSE69528 datasets using the 'limma' R package. The |log2 fold change (FC)| ＞ 1 and P-value < 0.05 were regarded as the cut-off criteria. Volcano plot was drawn to visualize the differences in gene expression levels between the normal and sepsis samples. After overlapping with 1811 MAGs, we obtained the 75 DEMAGs.

Construction and verification of the prognostic signature
We randomly divided patients (648 samples) of GSE65682 datasets randomly into a training set (n=328) and a validation set (n=140) at a ratio of 7:3. The univariate Cox regression model was used to identify DEMAGs that were significantly associated (P < 0.05) with the OS of the sepsis patients. Subsequently, the candidate DEMAGs were subjected to multivariate Cox regression analysis to evaluate their contribution as prognostic factors in patient survival. The risk score calculating formula was: Risk score = ExpGene1*Coef1 + ExpGene2*Coef2 + ExpGene3*Coef3…where Coef means the regression coefficients of genes, Exp is the normalized expression values of each gene signature.
To stratify patients into high-and low-risk groups, the optimum cutoff value for the risk score was determined using the 'survminer' package in R. Next, the K-M survival curve and log-rank test were performed to evaluate the OS between high-and low-risk groups. The area under the curve (AUC) of receiver operating characteristic (ROC) was calculated using the 'survival ROC' package in R. In addition, the risk plot was illustrated using the 'pheatmap' package in R.

Functional enrichment analysis (GSEA)
The 468 sepsis patients of the GSE65682 dataset were divided into high-and low-risk groups based on the median value of risk score. Then, we performed GSEA on all genes between the two risk groups using the clusterProfiler of R. Biological functions of GO and KEGG terms with the threshold of P-value < 0.05 was considered as significant.

Evaluation of immune microenvironment and immunotherapy efficiency
The immune-, stromal-and estimate scores from each sepsis patients were calculated by applying the 'Estimation of Stromal and Immune cells in MAlignant Tumours using Expression data' (ESTIMATE) algorithm. Spearman correlation analysis was used to evaluate the correlation between risk score and the immune-, stromal-and estimate scores, respectively.
Microenvironment cell population-counter (MCP-counter), a methodology based on gene expression profile data, was used to evaluate absolute abundance of 8 immune and 2 nonimmune stromal cell populations. In addition, the SubMap analysis was applied to predict the clinical response to the immune checkpoint blockade (CTLA-4, PD-1 and PD-L1). TIDE (http://tide.dfci.harvard.edu/) was used to estimate TIDE prediction scores with normalized transcriptome data from each patient. A low TIDE prediction score represents weak potential immune escape, and these patients would potentially exhibit a greater immune therapy response [16]. Thus, TIDE was used to predict the efficacy of immunotherapy in two risk groups.

Statistical analysis
Statistical analysis was performed using the R programming language and GraphPad Prism version 8.0.2 software. Student's t-test was used to analyze real-time gene expression analysis. In all analyses, P-value less than 0.05 were considered as statistically significant.

Identification of DEMAGs in Sepsis
To identify genes that may play critical roles in the development and progression of sepsis, the datasets (GSE65682 and GSE69528) were analyzed respectively to identify DEGs between normal and sepsis samples. A total of 1165 DEGs were identified in the GSE65682 dataset, of which 400 were upregulated and 765 were downregulated (Fig. 1A, Supplementary materials  Table S4).

Selection of prognosis-related DEMAGs
We conducted a univariate Cox regression to investigate the correlation of the 75 DEMAGs with the OS of sepsis patients in the GSE65682 training set and identified 15 DEMAGs significantly related to OS in sepsis patients ( Fig.2A, Supplementary materials Table S5, P < 0.05).
Subsequently, we conducted a stepwise multivariate Cox regression analysis. Six of the 15 DEMAGs were extracted and considered as significantly connection with the OS of sepsis patients.

Evaluation and validation of the prognostic risk model
We calculated the prognostic risk score for each patient in the training set. Patients were divided into the high-and low-risk groups according to the median value of risk score. K-M survival analysis demonstrated that patients with high risk scores had significantly poorer OS than patients with low risk scores (Fig.3A, P < 0.0001). We drew a ROC curve to indicate the accuracy of the risk model, and the AUC values of the risk model for 28 days was 0.696, which indicated that the risk model had good sensitivity for survival prediction (Fig.3B). The survival status of the patients and the expression of the 6 DEMAGs were shown in Fig.3C-E. With the increase of risk score, the mortality rate of patients in the high-risk group was significantly higher than that in the low-risk group. To validate the prognostic value of the 6-gene signature, the patients of validation set was divided into high-and low-risk groups according to the same cut-off value of training set, and the results were consistent with the training set, as shown in Fig.3F-J. Taken together, these data suggest that the 6 key DEMAGs could effectively predict the prognosis of sepsis patients.

Identification of the gene signature-related biological processes and pathways
The GSEA analysis was used to elucidate the biological functions and pathways associated with  Table S7). The KEGG analysis revealed that significant enrichment pathway were Natural killer cell mediated cytotoxicity, T cell receptor signaling pathway, and Antigen processing and presentation  Table S8). All these demonstrated the prognostic signature is highly positive associated with the immune related biological pathway.

Comparison of the immune microenvironment between high-and low-risk groups
In order to evaluate the discrepancy of immune microenviroment between the high-and low-risk groups, the stromal-, immune-, and estimate scores were compared between high-and low-risk groups. The results showed that stromal-, immune-, and estimate scores in the high-risk group were significantly lower than that in the low-risk group (Fig.5A-C). At the same time, we found that risk score has a negative correlation with immune-and estimate scores (Fig.5D). Furthermore, we analyzed the proportion of 10 immune cells in the high-and low-risk groups, as shown in

Immune checkpoint receptors including programmed death-1 (PD-1) and cytotoxic T lymphocyte antigen-4 (CTLA-4) have been shown to be increased on immune cells during sepsis
and hypothesized to be one of the major contributors causing sepsis induced immune cell dysfunction [17]. Therefore, we used the SubMap algorithm to estimate the clinical response of the high-and low-risk groups to immune checkpoint blockade (CTLA-4 and PD-1). Interestingly, we found that the low-risk group patients had a sensitive respond for anti-PD-1 therapy (P < 0.05; Fig.6B). In addition, the TIDE algorithm was utilized to calculate the immunotherapy response rates in the two risk groups. And the TIDE value was significantly higher in the high-risk group compared with that in low-risk group, indicating that patients with low risk score had more sensitive to respond to immune therapy (P < 0.05; Fig.6C).

Experimental verification of 6 DEMAGs in mRNA level
To confirm the reliability of the bioinformatics analysis, we validated the expression level of the six DEMAGs in sepsis mice model by quantitative real-time PCR. Consistent with the results of bioinformatics analysis, the expression of NQO2, ELANE and TLR2 were significantly upregulated in the(CPL), while the expression of CD3E and SMAD3 were significantly downregulated in CPL group. The expression of PTGDS was not consistent with the results of bioinformatics analysis, which showed upregulated in CPL group (Fig.7A). In addition, the expression of the six genes in GSE65682 and GSE69528 are shown in the figure (Fig.7B-C).

Discussion
Despite advances in treatment strategies for sepsis in recent years, it remains a leading cause of death worldwide due to the absence of specific treatments and drugs [3]. The pathologic process of sepsis is complex and subject to the interplay of multiple pathways, such as hypoxia [18], immune activation, and changes of metabolic processes [19]. Among them, the changes of metabolic processes on sepsis has been well understood in recent years. Therefore, it seems logical to develop new predictive or therapeutic targets based on genes involve in metabolic processes. In this study, we focused on the MAGs and firstly constructed a six-gene signature in sepsis. Then, we demonstrated that the 6-gene signature has a good predictive value for patient prognosis, immune microenvironment and responsiveness to immunotherapy. This study has some implications for the mechanism of sepsis research and clinical practice.
In the present study, combining the two independent GEO datasets and metabolism associated gene set, we obtained a 6-gene signature (including ELANE, NQO2, TLR2, PTGDS, SMAD3, CD3E). In the 6-gene signature, ELANE, NQO2 and TLR2 were significantly up-regulated, while PTGDS, SMAD3 and CD3E were significantly down-regulated in sepsis. All the above results were verified in the sepsis mice model. The expression levels of five genes, ElANE, NQO2, TLR2, CD3E and SMAD3, in sepsis mice model were consistent with in sepsis patients. However, the expression level of PTGDS was contrary to our silico mRNA expression in sepsis patients. This result suggests that PTGDS may play different roles in sepsis patients and mice model. In addition, we found that of these six genes, ELANE was associated with poor outcomes, while the remaining five genes were associated with favorable outcomes. Combined with previous studies, we found that ELANE [20][21][22], TLR2 [23][24][25], and SMAD3 [26,27] have been clearly reported in sepsis, while the roles of CD3E, NQO2 and PTGDS in sepsis remain unclear. Previous studies have shown that the expression of ELANE is able to reflect the leukocyte function during sepsis, which is able to provide more prognostic information than leukocyte counting [21]. Brittney Williams et al found that TLR2 and TLR7 are involved in coagulation dysfunction in murine sepsis, which may be an important reason for the high mortality of sepsis [28].The regulatory axis of GDF / SMAD2 / SMAD3 in sepsis that can promote macrophage polarization toward M2 phenotype leading to immunosuppression, and activation of this regulatory axis may be associated with high lethality in sepsis [26]. CD3E, as part of the TCR-CD3 complex, is present on the surface of T-lymphocytes and plays an important role in adaptive immunity [29]. In agreement with the findings of Raquel Almansa et al., the expression of CD3E was inversely associated with sepsis mortality [30]. The related research reports on PTGDS showed that PTGDS was mainly related to nervous system diseases [31,32]. In the inflammatory response, PTGDS is involved in eicosanoid synthesis and arachidonic acid metabolism [31,33]. NQO2 (N-Ribosyldihydronicotinamide:Quinone Reductase 2) encodes a member of the thioredoxin family of enzymes and diseases associated with NQO2 include breast cancer [34][35][36] and alcoholic pancreatitis [36]. Previous studies [37] have shown that NQO1 and NQO2 combined protect mice against lung inflammation as well as lung injury induced by hyperoxia, which was consistent with our result that NQO2 was a protective factor (Hazard ratio=0.59) for sepsis patients.
Sepsis is characterized by the immune system out of control, specifically the early over activation of the immune system and subsequent immunosuppression [38]. In fact, studies have shown that more than 60% of patients die in the immunosuppressive stage of sepsis, consistent with defective host immunity [38,39]. In the present study, we divided the patients into high-or low-risk groups based on the score of the six-gene signature, and performed GSEA enrichment analysis to explore the biological function and pathways related to the prognostic signature. Interestingly, the results suggested that prognostic signature was significantly associated with immune-related biological functions. The results of GSEA enrichment analysis are consistent with conclusion of previous studies, that is, immunosuppression in sepsis is closely related to the risk of poor outcomes [38].
In addition to this, we found that the immune microenvironment of sepsis patient in high-risk and low-risk groups was significantly different, and the composition of immune cells was significantly correlated with the risk score. The results of MCP-counter further distinguished the different types of immune cells in each group. Consistent with previous studies, we found that the proportion of T cells [40], B cells [41,42], NK cells [43] and neutrophils [44] was significantly lower in patients with higher risk. We speculate that this phenomenon may also be closely related to the immunosuppression of sepsis. In addition, the proportion of myeloid dendritic cells, endothelial cells and fibroblasts was higher in the high-risk group. Among the three types of cells, the role of myeloid dendritic cells in sepsis has been poorly studied, and the functions of endothelial cells [45,46] and fibroblasts [47] have only been partially explored. Therefore, the changes in the proportion of these immune cells in different prognostic septic patient need to be further explored.
Based on the conclusion that there was significant differences of immune microenvironment between high-and low-risk groups, we further explored the responsiveness of the two risk groups to immunotherapy. We compared the expression levels of some immune checkpoint molecules commonly used in the field of tumor immunotherapy, hoping to provide enlightenment for immunotherapy of sepsis [48]. The results shown that the expression levels of PD-1 and PD-L2 were significantly up-regulated in the high-risk group, which was contrary to the trend of CD27, HAVCR2, ICOS, LAG3 and TIGIT. Immunosuppression is considered to be an important cause of high mortality in sepsis [49]. Numerous studies have shown that immune checkpoint molecules in sepsis are significantly upregulated, which may be an important reason for immunosuppression [48][49][50]. The regulatory axis formed by PD-1 and its ligand PD-L1 / PD-L2 is the most characteristic interaction in immunopathology of sepsis. In the field of sepsis treatment, there are also an increasing number of preclinical and clinical studies showing that targeting both PD-1 and PD-L1 can improve host anti infectivity [51][52][53][54]. Recently, the first clinical evaluation of a PD-1/PD-L1 inhibitor showed that the drug was well tolerated, found no evidence of cytokine storm and showed signs of recovery of immune status within 28 days at high doses [55,56]. In addition, studies have shown that LAG3 can resist the role of PD-1 in blocking tumor therapy [57], and Niu et al. also observed the opposite expression trend of LAG3 and PD-1 in sepsis and its important role in T cell exhaustion [58]. However, the role of the other immunocheckpoint molecules plays in the immunotherapy of sepsis are still limited and needs to be further explored.
Furthermore, the results of TIDE and SubMap algorithm also showed that low-risk group may be more sensitive to immune checkpoint blockade treatment. Interestingly, we found that although the expression of PD-1 was lower in the low-risk group, while this group was more sensitive to the treatment of PD-1 inhibitor. Thus, we speculated that the expression of immune checkpoint molecules may not be directly related to the sensitivity of immunotherapy to this molecule, and the correlation between them remains to be further explored.
In this study, we successfully constructed a 6-gene signature based on MAGs. The signature could be used to predict survival and responsiveness to immunotherapy of sepsis, which may contribute to the development or optimization of new treatment for sepsis. However, there are some limitations in our study. Firstly, TIDE is an algorithm developed based on the data of melanoma and non-small cell lung cancer to predict anti-PD-1 and anti-CTLA4 treatment reactivity, which may produce bias when applied to sepsis [16,59]. However, considering that the mechanism of immune escape in tumor is similar to the mechanism of immune suppression in sepsis, and there is no better method to predict the response of sepsis immunotherapy to our best knowledge, we apply this algorithm to sepsis in order to make some exploratory contributions to the field of immunotherapy of non-neoplastic diseases [60][61][62]. Secondly, we only verified the expression of these six genes in the sepsis mice model, and the potential biological mechanism remains to be further explored. Thirdly, other types of sepsis, such as urinary tract infection or purulent meningitis, were not included in this study due to limited dataset. These parts of patients not included in the study need further analysis.

Conclusion
In conclusion, we identified and validated a prognostic signature for patients with sepsis, which likely reflects the immune dysregulation in the tumor microenvironment and is a potential prognostic biomarker and therapeutic target.