Identification of Potential Biomarkers of Septic Shock Based on Pathway and Transcriptome Analyses of Immune-Related Genes

Immunoregulation is crucial to septic shock (SS) but has not been clearly explained. Our aim was to explore potential biomarkers for SS by pathway and transcriptional analyses of immune-related genes to improve early detection. GSE57065 and GSE95233 microarray data were used to screen differentially expressed genes (DEGs) in SS. Gene Ontology and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway enrichment analyses of DEGs were performed, and correlations between immune cell and pathway enrichment scores were analyzed. The predictive value of candidate genes was evaluated by receiver operating characteristic (ROC) curves. GSE66099, GSE4607, and GSE13904 datasets were used for external validation. Blood samples from six patients and six controls were collected for validation by qRT-PCR and western blotting. In total, 550 DEGs in SS were identified; these genes were involved in the immune response, inflammation, and infection. Immune-related pathways and levels of infiltration of CD4 + TCM, CD8 + T cells, and preadipocytes differed between SS cases and controls. Seventeen genes were identified as potential biomarkers of SS (areas under ROC curves >0.9). The downregulation of CD8A, CD247, CD3G, LCK, and HLA-DRA in SS was experimentally confirmed. We identified several immune-related biomarkers in SS that may improve early identification of disease risk.


Introduction
Septic shock (SS) is defned as an infection-related circulatory dysfunction and metabolic disorder and is clinically associated with myocardial dysfunction and decreased ejection fraction [1,2]. Patient signs and symptoms include fever, rapid heart rate, shortness of breath, weakness and sweating, hypoxia, and altered mental status [3]. Kidney failure, malignancy, diabetes, chronic lung disease, congestive heart failure, and immunosuppression may also increase the risk of SS [4]. Currently, SS is the leading cause of death among hospitalized patients, with a fatality rate of up to 40% [5,6]. Although tissue perfusion can be quickly regained in patients with SS after a positive fuid resuscitation and symptomatic treatment with vasoactive drugs and anti-infection agents, the risk of recurrence after hospital discharge is high [7]. Furthermore, more than half of the survivors have impaired physical or neurocognitive function, mood disorders, and poor quality of life [8,9]. To reduce the mortality rate and improve the quality of life of patients after SS resuscitation, clinical studies are ongoing; however, early untreatable multiorgan failure makes the recovery process difcult. Terefore, in addition to improving programmed management, increasing the specifcity of early recognition can also efectively improve survival in SS.
Innate immunity and adaptive immunity play a key role in the response to SS. During the induced infammatory response, monocytes, macrophages, and neutrophils are activated, exacerbating vascular damage by producing cytokines, proteases, kinases, and reactive oxygen species [10]. Immunosuppressive responses are active in patients with SS [11], and the apoptosis of B cells and follicular dendritic cells is involved in this process [12]. T lymphocytes contribute markedly to the immune system, and T-cell abnormalities have been found in patients with SS [13]. Furthermore, immune dysfunction has been shown to impair the ability to clear primary infections and increase secondary infections [14,15]. Tolsma et al. found that neutropenia and specifc immunodefciency are independently associated with an increased risk of death in patients with SS [16]. However, the mechanism by which immune dysfunction contributes to the progression of SS remains unclear despite the important implications for the development of biomarkers.
In this study, we analyzed the immune-related pathways involved in SS progression and identifed candidate biomarkers. First, we screened diferentially expressed genes (DEGs) between patients with SS and healthy controls based on gene expression data from public databases, followed by enrichment analyses of DEGs. An immune cell deconvolution analysis and pathway-immune cell correlation analysis were performed to identify key immune-related genes. Finally, the expression of key genes was validated by quantitative real-time polymerase chain reaction (qRT-PCR) and western blotting. A fowchart of the study is displayed in Figure S1. Our fndings are expected to contribute to the development of personalized immunotherapy regimens for patients with SS.

Data Acquisition.
Two microarray datasets, GSE57065 and GSE95233, were downloaded from the Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/) [17]. GSE57065 comprised data for 28 patients whose blood samples were collected within 30 min, 24 h, and 48 h after SS and 25 healthy controls. Te GSE95233 included 22 healthy controls and 51 patients with SS, sampled twice at admission and once at D2 or D3. All samples from the GSE57065 and GSE95233 datasets were detected using the GPL570 [HG-U133_Plus_2] Afymetrix Human Genome U133 Plus 2.0 Array. GSE4607 and GSE13904 were used for the external expression validation of candidate genes between healthy controls and SS. GSE66099, which includes data for 30 systemic infammatory response syndrome (SIRS) and 181 SS samples, was used to compare the expression diferences between SS samples and the samples of the diseases with similar or related phenotypes to SS.

Data Preprocessing and Gene
Annotation. Based on the two datasets, the probe expression matrix after normalization and log2 transformation was downloaded. Tereafter, the annotation fle from the detection platform was obtained to match the gene symbol by the probe number inside; probes that did not match a gene symbol were removed. If diferent probes were mapped to the same gene symbol, the mean value of the probes was used as the fnal expression value for the gene.

Screening for DEGs between SS Samples and Controls.
Te limma package version 3.10.3 in R [18] (http://www. bioconductor.org/packages/2.9/bioc/html/limma.html) was used to analyze diferences in gene expression between SS samples and healthy controls. Te Benjamini & Hochberg (BH) method [19] was used to calculate the adjusted p value. Genes with an adjusted p value <0.05 and |logfold change (FC)| > 1 were selected as DEGs.

Protein-Protein Interaction (PPI) Network Construction.
Te intersecting DEGs between the GSE57065 and GSE95233 datasets were obtained to predict the PPI using STRING version 10.0 [23] (http://www.string-db.org/). During this analysis, the species was set to Homo, and the highest confdence score was set to 0.9. Cytoscape version 3.4.0 [24] (http://chianti.ucsd.edu/cytoscape-3.4.0/) was used to visualize the PPI network, and CytoNCA plug-in version 2.1.6 [25] (http://apps.cytoscape.org/apps/cytonca) was used to analyze the node degree with parameter setting "without weight." 2.6. Gene Set Enrichment Analysis (GSEA). Te R package clusterProfler version 3.16.0 [26] (http://bioconductor.org/ packages/release/bioc/html/clusterProfler.html) was used to perform GSEA. Te KEGG pathway gene set in the molecular signature database [27] (MSigDB, http://software. broadinstitute.org/gsea/msigdb/index.jsp) was used as the background gene set. DEGs between SS samples and controls with adjusted p < 0.05 and |logFC| > 0.263 were selected as input genes and were sorted in descending order by logFC. BH was fnally used to adjust p values, and pathways with adjusted p values less than 0.05 were considered statistically signifcant.  Table S1. Te reaction conditions were as follows: 50.0°C for 2 min, 95.0°C for 10 min, and 40 cycles at 95.0°C for 15 s and 60.0°C for 60 s. Te relative expression level was detected in triplicate and was normalized using the 2 −ΔΔct method.

Statistical Analysis.
Tanon Image was used to analyze the gray level of the western blot data. Te qRT-PCR and western blot results were analyzed using GraphPad Prism 5 (GraphPad Software, San Diego, CA, USA) and visualized using histograms. Diferences in the mRNA and protein expression levels of candidate genes between the SS and control groups were analyzed by t-tests. Statistics p < 0.05, p < 0.01, and p < 0.001 represent signifcant, highly signifcant, and extremely signifcant, respectively.

Screening of DEGs between SS Samples and Healthy
Controls. Based on the expression matrix of GSE57065 and GSE95233, principal component analysis (PCA) was performed (Figures 1(a) and 1(b)), which revealed a sharp distinction between the SS and control groups in both datasets. Tere were no signifcant outliers, and all samples could be used for further analyses. In total, 763 and 933 DEGs were obtained from GSE57065 and GSE95233, respectively, as shown in Figures 1(c) and 1(d).
DEGs with |logFC| > 2 were selected to construct a heatmap (Figures 1(e) and 1(f )); these genes were identifed as signifcantly diferentially expressed between the SS and control groups.

Construction of a PPI Network.
Te intersection of DEGs in GSE57065 and GSE95233 included 550 genes (Figure 3(a)). Tese DEGs were used to construct a PPI using STRING. A total of 1,104 interactions involving 243 proteins were obtained (Figure 3(b)). Nodes with more connections had larger contributions to the PPI network.

GSEA.
With GSEA, 5 upregulated and 13 downregulated pathways were obtained in GSE57065, whereas 7 upregulated and 15 downregulated pathways were obtained in GSE95233. According to the normalized enrichment scores (NES), the top six upregulated pathways, including Alzheimer's disease, complement and coagulation cascades, oxidative phosphorylation, and Parkinson's disease, and the top six downregulated pathways, such as graft versus host disease, primary immunodefciency, and T-cell receptor signaling pathway, were selected in both datasets and are displayed in Figure 4.

Diferential Analysis of Enriched Pathways.
Te enrichment scores for KEGG pathways for each sample were calculated using the GSVA algorithm. Based on the diferential analysis of enriched pathways, both GSE57065 and GSE95233 had 15 signifcantly downregulated KEGG pathways. As shown in heatmaps (Figures 5(a) and 5(b)), enrichment scores for these pathways difered signifcantly between the SS and control groups.

Deconvolution Analysis of Immune and Stromal Cells.
To evaluate diferences in immune cell enrichment between the two groups, we obtained the enrichment scores for 64 types of immune cells and stromal cells in each sample and compared SS and healthy controls by a Wilcoxon test (Figures 6(a) and 6(b)). In total, 11 and 14 types of cells had signifcantly diferent enrichment scores between the two groups in GSE57065 and GSE95233, respectively. By considering the intersection (Figure 6(c)), six types of cells, including megakaryocytes, CD8 + T cells, CD4 + TCM, preadipocytes, osteoblasts, and epithelial cells, difered signifcantly between SS samples and controls in both datasets.

Correlations between Pathways and Immune Cells.
Te intersection of signifcant pathways obtained by GSEA and GSVA included 12 KEGG pathways. Spearman correlation coefcients were obtained to evaluate the relationships among these 12 key pathways and frequencies of six signifcant cell types in each sample. Accordingly, 41 and 27 cell-pathway relationships were obtained in GSE57065 and GSE95233, respectively. Finally, 21 overlapped cellpathway relationships with signifcant positive correlations were identifed. Tese 21 cell-pathway relationships contained 10 key pathways and 3 key cells, including CD8 + T cells, CD4 + T cells, and preadipocytes; their relationships in GSE57065 and GSE95233 are summarized in Table 1.

Verifcation with External Datasets.
To verify the expression diferences in the 17 candidate genes between healthy controls and SS samples, GSE4607 and GSE13904 datasets were used for external verifcation. In the GSE4607 dataset (Figure 8(a)), all genes showed signifcant diferences in expression between the two groups, of which only MAPK14 was upregulated in SS. In the GSE13904 dataset (Figure 8(b)), MAPK14 was still overexpressed in SS, and the other 15 genes were signifcantly downregulated in the SS samples, except HLA−DQA1. Te GSE66099 dataset was subsequently included to verify whether these genes have expression specifcity in SS that enables their distinction from other types of infammation. Expression levels of  CD3E, CD3G, FYN, HLA-DPA1, HLA-DPB1, and HLA-DRA difered signifcantly between SIRS and SS (Figure 9), indicating that these six genes may contribute to SS infammation.

Experimental Expression Validation.
A total of 12 whole blood samples (6 cases and 6 controls) were collected to validate the expression of key genes. Among the 17 key genes, those with the highest degrees of connections in the PPI network (LCK and HLA-DRA) and genes with the highest fold change values in the GSE57065 or GSE95233 datasets (CD8A, CD247, and CD3G) were selected for validation by qRT-PCR. As depicted in Figure 10 in SS samples than in the control group. As determined by western blotting (Figure 10(b)), the protein expression levels of LCK and HLA-DRA were signifcantly lower in patients with SS than in healthy controls.

Discussion
Sepsis and SS are life-threatening diseases caused by a dysregulated immune response to infection and may lead to tissue and organ damage and even death [31,32]. Currently, there is no efective treatment for SS. Accordingly, the disease burden can only be reduced by early detection, resuscitation, and the prompt administration of appropriate antibiotics [9]. Unlike sepsis, SS leads to uncontrolled and intensifed infammatory responses, but the exact timing of triggering this process is elusive, which underscores the importance of identifying gene expression diferences between SS and normal controls, rather than sepsis, for diagnosis in patients early in disease progression [33]. Terefore, this study compared mRNA expression profles between SS and normal samples and then identifed 17 immune-related genes involved in the progression of SS    [34,35], but our study difers in that we focus more on the gene set involved in SS-related immune regulation. Terefore, we also proposed three key immune cells including CD4 + T cells, CD8 + T cells, and preadipocytes that may be regulated by immune-related candidate genes and involved in disease progression in SS. Tese potential immune regulatory mechanisms are important clues for understanding the role of candidate genes in disease diagnosis and for developing new drugs to prevent SS. In this study, DEGs between SS samples and controls were mainly enriched in biological functions and pathways related to immunological and infammatory responses. Using the GSVA algorithm, signifcant diferences in immune-related functions and pathways were also found between them, including the T-cell receptor signaling pathway, autoimmune disease, and primary immunodefciency, among others. A related bioinformatics analysis supported our fndings and demonstrated that DEGs in SS were involved in immune response, chemokine-mediated signaling, neutrophil chemotaxis, and chemokine activity [36]. Additionally, immunodefciency is commonly observed in patients with severe sepsis and SS and is associated with an increased risk of short-term mortality [16].
Considering the role of immunodefciency in SS development, we carried out deconvolution and correlation analyses to determine the efect of immune cells on SS. Based on our results, three key immune cell types, i.e., CD4 + T cells, CD8 + T cells, and preadipocytes, difered signifcantly between SS samples and controls. In terms of adaptive immunity, sepsis-induced apoptosis leads to lymphocytopenia in patients with SS, and this process involves all types of T cells, including T regulatory cells, CD4 + T cells, CD8 + T cells, and natural killer cells, which are conducive to immunosuppression [14]. Immunosuppression is a compensatory anti-infammatory response that explains the short-term death of SS patients, while survivors may experience a prolonged state of immunosuppression, which could be reactivated by pathogenic infection [37]. During the immunosuppression in SS, the loss of T-cell function is associated with reduced resistance to secondary infections in patients with SS [38]. Furthermore, decreased expression of cytotoxic molecules weakens the lytic activity of CD8 + T cells [39]. In this study, the enrichment scores for CD4 + T cells and CD8 + T cells were found to be signifcantly lower in SS cases than in controls. Roger et al. further supported our fndings and proposed that rates of CD4 + T cell and CD8 + T cell apoptosis were higher in patients with SS than in controls [40]. Te above evidence indicated that DEGs identifed in this study may trigger immunosuppression and lead to SS by down-regulating CD4 + T cell and CD8 + T cell levels. In addition, we also proposed a relationship between preadipocytes and immune defciency and T-cell receptor-related pathways. Several factors secreted by preadipocytes have pro-infammatory and anti-infammatory efects and can contribute to diseases associated with immune system dysfunction [41]. Some immunodefciency virus protease inhibitors inhibit preadipocyte diferentiation and promote adipocyte death [42]. In this study, we found several genes involved in the Tcell receptor signaling pathway, of which FYN is expressed in human preadipocytes and is induced after the initiation of        diferentiation [43]. Terefore, we hypothesized that preadipocytes participate in the activation of immune-related pathways during SS development by inducing the expression of FYN; however, further mechanistic investigations are still needed. Furthermore, 17 key genes were identifed in immunerelated cell-pathway pairs, and core genes from the PPI network were selected for expression validation. Relevant results suggested that the mRNA expression levels of CD8A, CD247, and CD3G were downregulated in SS samples, while LCK and HLA-DRA were decreased at both the mRNA and protein levels. As a member of the Src family of kinases, LCK is involved in changes in the activity of CD4 + T cells and CD8+ T cells during T-cell development [44]. LCK plays a crucial role in T-cell diferentiation, survival, and activation [45]; however, the contribution of LCK to the development of SS has not been explored. With regard to HLA-DRA, studies have found that the reduced expression of HLA-DR mRNA is correlated with increased mortality after SS [46]. Winkler et al. reported that HLA-DR expression is decreased in sepsis [47], and the reduced HLA-DR expression may be a characteristic feature of septic monocytes

16
Genetics Research [14]. Furthermore, the dynamic changes in HLA-DRA gene expression and helper T cell subsets in patients with sepsis are indicative of immunosuppression [48]. Combined with the results of this study, we speculated that the loss of LCK and HLA-DRA expression may lead to the failure of T-cell diferentiation and dynamic changes of helper T cell subsets, thus resulting in immunosuppression and the onset of SS.
Te main limitation of this study is that patients with autoimmune diseases were excluded from the independent clinical validation cohort, which to some extent afects the extrapolation of the results. Furthermore, we only performed a preliminary exploration of the potential roles of these genes in disease. Te regulatory mechanisms by which these candidate genes contribute to immunosuppression during SS are not clearly established. In future studies, bioinformatic analyses will be used to predict the upstream regulatory mechanisms, and experimental approaches will also be carried out to confrm the regulatory efects of the candidate genes. Additionally, due to limited sample size, we did not further screen the most robust biomarkers out from the 17 identifed key genes by LASSO regression and/or multivariate Cox analyses. In the future, a most valuable predictive signature should be developed based on more convincing methods and larger sample size.
In conclusion, we found that DEGs between SS cases and controls were mainly enriched in immune-and infammation-related functions and pathways. In addition, CD4 + T cells, CD8 + T cells, and preadipocytes were proposed as key immune cells to involve in the SS progression. Tese immune cells were also associated with 17 key immune-related genes, among which the downregulation of CD8A, CD247, CD3G, LCK, and HLA-DRA in SS samples was further experimentally validated. Our fndings reveal several novel biomarkers for the early identifcation of SS.

Data Availability
Microarray datasets including GSE57065, GSE95233, GSE4607, GSE13904, and GSE66099 used in this study can be downloaded from the GEO database at http://www.ncbi. nlm.nih.gov/geo/. Te raw data of qRT-PCR and western blot are available at https://doi.org/10.4121/19074482.