Identification of Potential Early Diagnostic Biomarkers of Sepsis

Introduction

Sepsis causes life-threatening organ dysfunction due to dysregulated host response to infection.^1,2 It is the most important cause of morbidity and mortality in the intensive care unit (ICU), and it causes substantial health care costs worldwide.^3,4 Recent reports revealed that 48.9 million cases of sepsis were recorded worldwide in 2017, with the 11 million sepsis-related deaths accounting for 19.7% of total global deaths.⁵ The incidence of sepsis varies considerably among different geographical regions.⁶ Different aspects of society, economy, politics, health, and even climate affect the epidemiology of sepsis.⁷ Community infections account for approximately 60–70% of all cases; nosocomial infections outside the ICU, 20–30%; and primary cases in the ICU, 5–9%.^8–10

Pediatric sepsis has a distinct epidemiological profile from adult sepsis and is associated with lower rates of mortality and progressive organ failure.¹¹ Nevertheless, adult and pediatric sepsis have the same bacterial and viral causes.¹² Interventions applied to adults with sepsis are often translated into clinical guidelines for pediatric sepsis patients as well.¹³

Sepsis survivors often suffer from multiple sequelae that significantly compromise their life expectancy and quality of life.¹⁴ Survivors are at increased risk of long-term cognitive impairment, new cardiac events, and new renal failure.¹⁵ Survivors often show progressive atherosclerosis, neurocognitive dysfunction and low-grade inflammation.^16,17

The acute onset of sepsis is characterized by massive release of cytokines and other mediators, which leads to dysregulated immune responses, resulting in organ damage or death.^18,19 Chronic inflammation, coupled with chronic catabolism of drugs, contributes to a state of immune dysregulation that promotes infectious complications and leads to chronic exacerbations and death.²⁰ A striking feature of immunosuppression in sepsis is that apoptosis leads to early and massive loss of lymphocytes,²¹ which if severe or persistent can increase risk of mortality.²² Clinical signs of early sepsis are often subtle, nonspecific, and easy to miss. Knowledge of early sepsis symptoms may help patients and clinicians identify the condition and initiate treatment sooner.^23,24

The present study analyzed gene expression profiles of adult and pediatric sepsis patients from public databases in order to understand early molecular changes as well as biological mechanisms. For the early diagnosis of sepsis patients, we established a least absolute shrink and selection operator (LASSO) model to develop a gene signature for predicting sepsis, which we validated using external data. Analysis of the area under receiver operating characteristic curve (AUC) allowed us to identify potential biomarkers of sepsis. The goal of the present study was to identify early diagnostic molecules and improve the survival of sepsis patients.

Materials and Methods

Data Sources

We collected sepsis data from gene expression omnibus (GEO) databases (http://www.ncbi.nlm.nih.gov/geo/). The GSE54514 dataset included gene expression profiling based on arrays of whole blood from 26 sepsis survivors (days 1–5), 9 sepsis patients who died, and 18 healthy controls. The GSE25504 dataset included gene expression profiling based on arrays of peripheral blood from 37 control neonates and 25 neonates with sepsis. Raw data were background-subtracted and normalized using a robust spline algorithm with the lumi package in Bioconductor.²⁵ Institutional Review Board approval was not required for this study since all data had been released into the public domain.

Analysis of Differential Gene Expression

First, the expression profiles of genes in each sample of sepsis and control groups were constructed. Then, genes differentially expressed between the sepsis and control groups in the GSE54514 and GSE25504 datasets were identified using the Limma R package.²⁶ Differentially expressed genes (DEGs) were defined as those whose expression differences were associated with P < 0.05.

Coexpression Analysis

Weighted gene coexpression network analysis (WGCNA) was used to identify the coexpression network of DEGs in GSE54514. The WGCNApackage in R was used to cluster genes into different modules.²⁷ Briefly, connections among pairs of genes were identified and weighted based on correlated expression levels across multiple samples. We further converted the adjacency matrix into a topological overlap matrix (TOM) to detect gene connectivity in the network. Then, hierarchical clustering was used to identify clusters (modules) of highly interconnected genes on the basis of their connectivity and covariance coefficients. The minimum size of the gene dendrogram was 30. The power value was selected to define a high degree of positive correlation among genes in the same module. A screening threshold crosstalk greater than 900 was set for quantitating crosstalk between modules.

Enrichment Analysis

Module genes were analyzed using Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) in the clusterProfiler package in R.²⁸ biological processes (BPs)was used to identify enrichment in GO. The clusterProfiler package was also used to conduct gene set enrichment analysis (GSEA).²⁹ Results associated with P < 0.05 were considered statistically significant.

The fgsea package in R was used to display the enrichment results of GSEA. Gene set variation analysis (GSVA) was carried out using the GSVA package in order to display enrichment results.³⁰ For each sample, a score was obtained describing the enrichment of a set of genes based on the gene expression profile. Then, the Limma package in R was used to compare GSVA scores between sepsis and controls. This approach identified signaling pathways that were up- or down-regulated in sepsis patients relative to controls.

LASSO Regression

Receiver operating characteristic curves were drawn and their AUCs determined to evaluate the diagnostic value of the key genes using the pROC package in R.³¹ Genes whose AUC > 0.7 and that were differentially regulated (up- or down-regulated) in the same direction in both GSE54514 and GSE25504 datasets were used to build a binomial LASSO model using the glmnet package in R. As a training set, 75% of samples in the GSE54514 dataset were randomly selected. As λ increases, LASSO tends to reduce the regression coefficient to zero. As a validation set, 25% of samples in the GSE54514 and GSE25504 datasets were selected. Genes with AUC > 0.75 in both the GSE54514 and GSE25504 datasets were defined as key genes.

Short Time-Series Expression Miner (STEM) Analysis

We performed STEM analysis³² of the GSE54514 dataset in order to cluster gene modules in healthy controls, sepsis survivors and sepsis patients who died. Clustering was significant if it was associated with P < 0.05. Significantly clustered genes showed a trend of greater up- or down-regulation in the order: healthy controls < sepsis survivors < sepsis patients who died.

Single-Sample Gene Set Enrichment Analysis (ssGSEA)

Marker genes for different types of immune cells were obtained from Bindea et al.³³ The ssGSEA method was used to derive an enrichment score for each immune-related cell type. The infiltration level of immune cell was calculated using ssGSEA function of the GSVA package in R. We analyzed differences in immune cell infiltration between sepsis and controls. We also analyzed differences in immune cell infiltration between sepsis survivors and controls, as well as between survivors and sepsis patients who died. Correlations between different infiltration of immune cells and genes were assessed using Pearson’s correlation.

Results

Coexpression Network of Sepsis-Associated Genes

The flowchart of this study is shown in Figure 1. To identify sepsis-associated genes, we analyzed genes differentially expressed between sepsis and control samples. A total of 6316 differentially expressed genes (DEGs) in adults were obtained in the GSE54514 dataset (Figure 2A, Table S1), and 7832 DEGs in neonates in the GSE25504 dataset (Figure 2B, Table S2). A total of 3355 DEGs were common to adult and pediatric sepsis patients, while 4477 genes may be specific to pediatric sepsis and 2961 genes may be specific to adult sepsis (Figure 2C, Table S3).

Figure 1 Study flowchart. Sequencing data from sepsis patients and controls in GSE54514 and GSE25504 datasets were analyzed by bioinformatics in order to identify early potential biomarkers of sepsis. Abbreviations: LASSO, least absolute shrink and selection operator; STEM, short time-series expression miner; WGCNA, weighted gene coexpression network analysis.

Figure 2 Coexpression network of differentially expressed genes. (A) Genes differentially expressed between sepsis patients and controls in the GSE54514 dataset. Green nodes represent down-regulation in sepsis; red nodes, up-regulation; and grey nodes, no significant difference from controls. The five genes most significantly up- or down-regulated, based on log2(fold change), are marked using symbols. (B) Genes differentially expressed between sepsis patients and controls in the GSE25504 dataset. Green nodes represent down-regulation in sepsis; red nodes, up-regulation; and grey nodes, no significant difference from controls. The five genes most significantly up- or down-regulated, based on log2(fold change), are marked using symbols. (C) Intersection of differentially expressed genes (DEGs) in the GSE54514 and GSE25504 datasets. The count on the left refers to DEGs unique to GSE54514; the count in the middle, DEGs common to both datasets; and the count on the right, DEGs unique to GSE25504. (D) Correlation between soft threshold power and scale-free topology model. (E) Cluster tree of coexpression modules of significantly different gene expression. Different colors represent different modules. (F) Crosstalk between modules. The more crosstalk between module genes and other genes, the greater the proportion of the ring is occupied by that module. Different colors represent different modules.

To identify the coexpression network of DEGs in the GSE54514 dataset, we performed WGCNA. The optimal soft power threshold for WGCNA was set to 10 to preserve scale-free topology and effective connectivity (Figure 2D). Ten coexpression modules of DEGs were then established (Figure 2E). Moreover, we found significant crosstalk between coexpression modules (Figure 2F).

Enrichment of Module Genes

Based on analysis of enrichment and GSVA of module genes, we found that interferon - γ secretion, adaptive immune memory response, and T cell proliferation of biological processes were significantly up-regulated in sepsis patients compared with controls, whereas neutrophil extravasation, pentose metabolic processes, and cytosine metabolic processes were significantly down-regulated (Figure 3A). Among KEGG pathway results, antigen processing and presentation, primary immunodeficiency, as well as Th1 and Th2 cell differentiation were significantly up-regulated in sepsis patients compared with controls. Conversely, the VEGF signaling pathway, cell cycle, and HIF-1 signaling pathway were significantly down-regulated (Figure 3B). Calculation of KEGG pathways in GSEA showed that the top 10 pathways were enriched mainly in sepsis, while the latter 10 pathways were enriched mainly in controls (Figure 3C). This revealed that immune response was significantly up-regulated in sepsis, while progression of metabolic processes was down-regulated.

Figure 3 Biological functions and KEGG pathways enriched for module genes. (A) Significant up- or down-regulated biological processes in module genes of sepsis patients relative to controls, as quantified by gene set variation analysis (GSVA). FC, fold change. (B) Significant up- or down-regulated KEGG signaling pathways in module genes of sepsis patients relative to controls, as quantified by GSVA. (C) Up- or down-regulated KEGG pathways of gene set enrichment results in sepsis patients relative to controls. A P value < 0.05 was considered statistically significant.

Identification of Key Genes Associated with Sepsis

We randomly split all the samples in the GSE54514 dataset into a training set (75%) and validation set (25%). A LASSO model was constructed in the training set by screening 51 genes with AUC > 0.7 from common genes that showed concordant up- or down-regulation in both the GSE54514 and GSE25504 datasets. A λ value of 25 was chosen to determine the genes that could predict sepsis most accurately (Figure 4A). We then obtained a 25-gene signature with a non-zero coefficient (Figure 4B). The AUC for the gene signature was 0.983 in the training set (Figure 4C) and 0.964 in the validation set (Figure 4D). Importantly, we also used external datasets to validate the gene signature: it showed good diagnostic ability, with an AUC of 0.841 (Figure 4E). Genes with AUC > 0.75 in both the GSE54514 and GSE25504 datasets were further screened, ie, CDKN1C, TMEM169, SLC2A6, DERL2, PSMB7, C1ORF55, DYNC1LI2, DUSP5, and RHOB (Figure 4F). Among them, SLC2A6, C1ORF55, DUSP5, and RHOB in the gene signature were recognized as key genes and may be able to diagnose sepsis. These key genes were up-regulated in sepsis patients compared to controls (Figure 4G).

Figure 4 Potential key genes for the diagnosis of sepsis. (A) The gene signature selection of optimal parameter (lambda) in LASSO model. (B) LASSO coefficient profiles of the 25 differentially expressed genes selected by the optimal lambda. (C) The receiver operating characteristic (ROC) curves of the gene signature in the training set of GSE54514. (D) The ROC curves of the gene signature in the validation set of GSE54514. (E) The ROC curves of the gene signature in GSE25504. (F) The genes in GSE54514 and GSE25504 with an area under the ROC curve (AUC) greater than 0.75 are indicated together with their mean AUCs. (G) Differential expression of key genes between sepsis patients and controls in GSE54514. ***P < 0.001.

Persistently Dysregulated Genes Associated with Sepsis Death

STEM software was used to screen module genes for persistently dysregulated genes. Among the 293 genes identified, SLC2A6 also showed a trend towards up-regulation in the order: healthy < sepsis survivors < sepsis patients who died (Figure 5A). Immune cell infiltration differed between sepsis patients and controls (Figure S1A). Infiltration by most types of immune cells was significantly higher in sepsis survivors than in controls, and lower in sepsis patients who died than in those who survived (Figure S1B). In these immune cells, CD8+ T cells are a subtype of cytotoxic T cells.³³ Among the correlations between immune cells and key genes, SLC2A6 showed the strongest positive correlation with Th1 cells (Figure 5B).

Figure 5 Persistently dysregulated gene expression during sepsis development. (A) The expression of SCL2A6 was persistently elevated during sepsis development. (B) Pearson correlation of immune infiltrating cells with the key genes. Red nodes indicate positive correlation, and blue nodes indicate negative correlation. *P < 0.05, **P < 0.01. (C) Heatmap of gene sets showing persistent up- or down-regulation that increased in the trend: healthy controls < sepsis survivors < sepsis patients who died. Gene sets were arranged based on cluster assignment to generate simplified expression profiles. We graphically depict only 4 modules with >40 genes. (D) The box plots of STEM genes in 4 clusters. Line plots and box plots were used to display, respectively, fold changes (log2FC) or absolute expression levels based on fragments per kilobase per million reads (log2 fragments per kilobase million). Representative genes were highlighted using red lines. The key genes were located on the right side of the box map. *P < 0.05, **P < 0.01. (E) Signaling pathways persistently up- or down-regulated as sepsis develops. Red nodes in the heatmap represent up-regulated signaling pathways, while blue nodes represent down-regulated signaling pathways.

Cluster analysis revealed that 253 persistently dysregulated genes significantly clustered into four modules (Figure 5C and D). Subtype GSEA results showed that in sepsis patients who died, the Toll-like receptor signaling pathway and cholesterol metabolism were up-regulated, while asthma and cell adhesion molecules were consistently down-regulated (Figure 5E).

Discussion

Sepsis causes life-threatening organ dysfunction and is a medical emergency associated with high mortality and long-term disability in survivors.³⁴ The present study analyzed differences in gene expression between sepsis patients and healthy controls and identified changes in immune and metabolic functions associated with these genes. Importantly, a combination of the LASSO model and AUCs identified SLC2A6, C1ORF55, DUSP5, and RHOB as having a diagnostic role in sepsis, especially SLC2A6, which may be a marker for early diagnosis and prognosis.

Identification of disease-associated modules by coexpression analysis has emerged as a powerful approach to gain new insights into disease biology.³⁵ Gene sets showing synergistic expression patterns were screened using WGCNA. The enrichment results showed that the immune response was significantly up-regulated in sepsis, while metabolism was significantly down-regulated. Interferon - γ (IFN - γ), which is produced by activated T and natural killer (NK) cells, is an important immunoregulatory factor.³⁶ In a rat model of sepsis, the survival time of rats is prolonged by increasing the number of T cells and NK cells as well as up-regulating IFN - γ.³⁷ In fact, administering IFN - γ to patients with sepsis to strengthen the immune response has emerged as a plausible treatment and has been associated with better prognosis.³⁸ Researchers have suggested that therapeutic interventions that suppress excessive inflammation, shift catabolism to anabolism, and enhance immune function may be beneficial once the initial symptoms of sepsis resolve.^39,40 Indeed, persistent immunosuppression can lead to sepsis and consequently worse outcomes and higher long-term mortality.⁴¹ Sepsis dramatically decreases the number of pre-existing memory CD4 T cells and impairs their function, such that sepsis survivors are at higher risk of secondary infections.⁴² The finding that the immune response was significantly up-regulated in our analysis may be related to the larger proportion of survival patients with sepsis and the fact that all survivors were in early stages of sepsis.

Restoration of immune system function is extremely important for the prognosis of patients with sepsis.⁴³ Th1 cells release mainly the proinflammatory cytokine IFN - γ, which induces cell phagocytosis and intracellular killing of microorganisms, and Th1 cells interact with Th2 cells to jointly regulate infection clearance.⁴⁴ The results of this study showed that Th1 cells were up-regulated in sepsis, while Th2 cells were down-regulated, and that sepsis patients with persistently higher ratios of Th2 to Th1 cell numbers had the highest mortality.⁴⁵ During the inflammatory process, VEGF promotes vascular dilation and permeability, disrupting the endothelial barrier.⁴⁶ VEGF-A levels are increased in patients with sepsis in the ICU, and VEGF-A levels correlate with the severity of sepsis.⁴⁷ HIF-1 α expression is elevated in patients with sepsis, and inhibition of HIF-1 α activity may be a novel therapeutic target for sepsis.⁴⁸

We obtained a 25-gene signature and further identified SLC2A6, C1ORF55, DUSP5 and RHOB as key genes based on the LASSO model and AUCs. The AUCs were validated in the GSE54514 and GSE25504 datasets, suggesting that they are potential diagnostic markers for sepsis not only in adults but also in children. Studies have shown that SLC2A6 is associated with lipopolysaccharide (LPS)-induced inflammatory responses.⁴⁹ SLC2A6 expression is up-regulated in cells infected with Gram-positive bacteria.⁵⁰ A common up-regulated gene after infection, SLC2A6, which is regulated by nuclear factor - κ B (NF - κ B), encodes a protein transporter involved in metabolic responses.^51,52 Expression of C1ORF55 has been associated with prognosis of breast cancer patients.⁵³ Although we are unaware of studies suggesting a link between C1ORF55 and sepsis, the results of our analysis suggest that it may be beneficial for the early diagnosis of sepsis. Dual specificity phosphatase 5 (DUSP5) is induced during LPS-mediated inflammation and inhibits the activity of NF-κ B.⁵⁴ In addition, DUSP5 specifically targets the inactivation of the ERK 1/2 signaling pathway.⁵⁵ DUSP5 can be significantly induced after LPS stimulation and mycobacterial infection, which is beneficial for granulocyte development.^56,57 Knocking DUSP5 out in mice promotes T cell proliferation and cell death.⁵⁸ RhoGTPases are important for maintaining the endothelial barrier, and during inflammation, RHOB negatively regulates Rac-1 membrane translocation and prevents Rac-1-mediated restoration of the integrity of injured endothelium.^59,60 It is clear that inflammatory and anti-inflammatory responses, as well as innate and adaptive immune systems, are equally important and are likely to be targets for future immunotherapy to improve long-term sepsis outcomes.^61,62

The potential limitations of this study need to be taken into account when interpreting the findings. First, the GSE54514 and GSE25504 datasets contain different populations of sepsis patients and controls, which may affect interpretation of the results. In addition, we were able to compare immune cell infiltration between sepsis patients who survived or died in the GSE54514 dataset, but further studies should analyze the two groups for additional differences. Second, because the data we analyzed were obtained from public databases, further experimental studies are necessary to validate the findings of this study.

Conclusion

In summary, our four key genes may have important implications for the early diagnosis of patients with sepsis. The analysis suggests that in the early stage of sepsis, the immune response is enhanced, while the metabolic response is weakened. Further analysis is needed to investigate the molecular mechanisms by which the four key genes affect the prognosis of patients with sepsis.