Identification of Key Genes Associated with Changes in the Host Response to Severe Burn Shock: A Bioinformatics Analysis with Data from the Gene Expression Omnibus (GEO) Database

Introduction

With the development of the economy and improvements in public safety awareness, the incidence of burn in less-developed and developing areas is gradually decreasing.¹ However, burns are still the fourth most common type of trauma worldwide, after traffic injuries, falls, and interpersonal violence, with up to 90% of the burns occurring in low- and middle-income areas.² According to an estimate released by the World Health Organization in 2018, 18,000 burn-related deaths occur annually.³ Burn injuries are associated with a high mortality rate in patients with severe burn shock which is a combination of distributive and hypovolemic shock accompanied by immune and inflammatory responses, and metabolic changes.⁴ Hence, the state of shock experienced following burn injuries can be a challenge in patient management and may result in multiple organ failure.⁵ Notably, the number of deaths within the burn shock period has been greatly reduced through improvements in the available treatment and resuscitation procedures for burn injuries.^6–8 However, the state of burn shock is associated with numerous complications, as patients experiencing an unstable period of shock will have a greatly weakened defense and repair system due to tissue ischemia, hypoxia, and reperfusion injury that accompany the state of shock. As a result, these patients are at a greater risk of experiencing concurrent severe systemic infections and multiple organ failure, which are currently the primary cause of death in severe burn cases.^9,10 Therefore, exploring the precise mechanism that underlies the host response to severe burn shock is crucial.

Numerous studies have indicated that gene expression plays an important role in different burn stages. It has been reported that the expression of high-mobility group box protein 1 (HMGB1) is elevated in burn patients and has an important impact on the immune function of postburn patients.¹¹ Carter found that the melanocortin 1 receptor (MC1R) gene polymorphism was associated with wound infection after burn injury.¹² In the repair stage after burn injury, nucleotide-binding domain leucine-rich repeat (NLR) and pyrin domain containing receptor 3 (NLRP3) promote wound healing through regulation of macrophage polarization and inflammation.¹³ However, there has been little information regarding the important genes expressed during burn shock. Given the popularity of gene detection technology, gene chips are widely used in research on many diseases.¹⁴ Bioinformatic analysis of the data generated from the use of gene chips has become a promising and effective tool for screening the significant genetic or epigenetic variations associated with the development and progression of diseases.^15,16 The differentially expressed genes (DEGs), signaling pathways, and hub genes that are potentially related to the host response to burn shock have not been identified. Moreover, to our knowledge, no studies have been performed to analyze the DEGs, enrichment analysis of Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, and protein–protein interaction (PPI) networks associated with burn shock using datasets obtained from the Gene Expression Omnibus (GEO) database.

This study, to our knowledge, is the first to identify 1059 DEGs by screening the gene expression dataset GSE77791, which included 30 samples from burn shock patients and 13 samples from healthy volunteers. The analyzed biological functions, signal pathways, constructed PPI network, and selected hub genes were used to understand the molecular mechanism underlying the host response to burn shock, which may potentially be insightful for studying burn shock, including early detection and treatment.

Materials and Methods

Gene Expression Microarray Data

We acquired the microarray expression profiles of severe burn shock patients from the National Center for Biotechnology Information GEO database (http://www.ncbi.nlm.nih.gov/geo), which contains high-throughput gene expression data, chips, and microarrays, under the accession number GSE77791. The GSE77791 dataset, based on GPL570 (a platform for Affymetrix Human Genome U133 Plus 2.0 Array), contains 117 samples, including 30 samples collected from patients with burn shock prior to any treatment and 13 samples collected from healthy volunteers that served as controls. These 43 samples were selected for deep analysis in this study in accordance with our aims. The remaining samples collected from patients with different medical conditions were excluded.

Data Preprocessing: Identification of DEGs

First, the probe sets were correspondingly converted into gene symbols by the platform with annotation information. The mean value of multiple probe sets was calculated if they corresponded to the same gene. The data were normalized, and DEGs were subsequently identified using the linear models for microarray data (LIMMA) package in R software (version 3.6.1). Individual p-values were calculated and converted to adjusted p-values (adj. p. val) for comparisons by false discovery rate correction of the Benjamini and Hochberg test. DEGs were selected with the commonly used thresholds of∣log₂fold change (FC)> 1 and adj. p. val < 0.05. The heatmap and volcano plot of the DEGs were constructed using the ggplot2 package in R.

Functional and Pathway Enrichment Analyses of DEGs

The functions of up- and down-regulated DEGs in severe burn shock were predicted based on GO function and KEGG pathway enrichment analyses and visualized using the ClusterProfiler package in R. The GO database (http://www.geneontology.org) primarily comprises the following three categories: biological processes (BP), cellular components (CC), and molecular functions (MF). The KEGG database (http://www.genome.ad.jp/kegg/) collects systematic functional, genomic, and chemical information. The significance criterion to screen the significantly enriched GO terms and KEGG pathways was set at p <0.05.

PPI Network Construction and Module Analysis

Functional interactions between proteins were analyzed to shed light on the exact mechanisms of generation and development of diseases. In this study, the PPI network of DEGs was predicted in an online database (Search Tool for the Retrieval of Interacting Genes database, STRING, version 11.0, http://string-db.org) and subsequently visualized in Cytoscape (version 3.7.2). DEGs were uploaded onto the database, and single protein nodes were removed. A combined score > 0.9 was set as the cutoff criterion. The PPI network was downloaded and further analyzed in Cytoscape. The important modules were identified and visualized by Molecular Complex Detection (MCODE), which is a plugin in Cytoscape for discovering densely connected nodes in a given network. The module selection was in accordance with the following criteria: k-score=2, MCODE score > 10, degree cut-off=2, max depth=100, and node score cut-off=0.2.¹⁷ The nodes in the key modules are presented as highly connected proteins that have important biological functions. Moreover, the pathway enrichment analysis of the DEGs in selected modules was performed using the ClusterProfiler package.

Hub Gene Selection

The hub genes in the PPI network were characterized by the cytoHubba plugin in Cytoscape. The top 10 genes were identified by three methods (Maximal Clique Centrality (MCC), Degree, and Maximum Neighborhood Component (MNC)) and overlapped to identify the hub genes.

Results

Patients in GSE77791

Thirty-two patients were enrolled in Plassais’s study,¹⁸ however, two of them were discarded due to a major batch effect related to a technical issue. The demographics of the patients are shown in Table 1. Briefly, these patients were severely burned with a total burn surface area (TBSA) range from 30 to 98 and registered for hospitalization between 24 and 72 h after injury. Few data about the severity of the burns in the original article were obtained: the mean abbreviated burn severity index (ABSI) score was 11, 12 patients had inhalation injuries, and 3 patients had infections during hospitalization. The exact time points of injury were not recorded. Samples were immediately collected from hospitalized patients with severe burn shock, which was defined by a norepinephrine dosage above 0.5 µg/kg/min.

Table 1 Demographics of the Patients in Jonathan Plassais’s Study

Identification of DEGs

The microarray data of 30 burn shock patients and 13 healthy volunteers from the GSE77791 dataset were analyzed to identify DEGs according to the following criteria:∣log₂fold change (FC)∣> 1 and adj. p. val < 0.05. A total of 1059 DEGs from the samples of patients with burn shock were identified by comparing with the controls. This total included 508 downregulated genes and 551 upregulated genes, which are indicated in the volcano plot shown in Figure 1A. The top 100 DEGs of all samples were clustered and represented as a heatmap (Figure 1B).

Figure 1 DEGs from GSE77791. (A) 1059 DEGs are showed in the volcano plot that containing 551 upregulated genes in red and 508 downregulated genes in green. (B) The heatmap shows the top 100 most significant DEGs. Red indicates a relatively high expression and blue indicates a relatively low expression. DEGs were identified by the criteria of∣log2fold change (FC)∣>1 and adj. p. val <0.05. DEGs, differentially expressed gene.

Analyses of GO Function and KEGG Enrichment of DEGs

Analyses of GO function and KEGG enrichment were performed using the ClusterProfiler package to further understand the functions and mechanisms of the identified DEGs. The top eight enriched terms in BP, CC, and MF are shown in Figure 2A. For BP, the DEGs were enriched in neutrophil activation, neutrophil activation involved in immune response, neutrophil degranulation, neutrophil-mediated immunity, T cell receptor signaling pathway, T cell activation, T cell differentiation, and lymphocyte differentiation. In terms of CC, the DEGs were enriched in specific granules, tertiary granules, secretory granule lumen, tertiary granule lumen, specific granule lumen, secretory granule membranes, vesicle lumen, and cytoplasmic vesicle lumen. Regarding MF, the DEGs were enriched in carbohydrate binding, MHC protein complex binding, immunoglobulin binding, cytokine receptor activity, MHC class II protein complex binding, non-membrane spanning protein tyrosine kinase activity, cytokine binding, and chemokine binding.

Figure 2 The GO terms and KEGG pathway enrichment analysis of all DEGs. (A) The top 8 terms enriched in BP, CC, and MF. (B) The top 8 enriched terms of KEGG pathway. Cutoff value is adj. p. val <0.05. GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; DEGs, differentially expressed gene; BP, biological process; CC, cellular component; MF, molecular function.

Based on KEGG pathway analysis, the DEGs were significantly enriched in Th1 and Th2 cell differentiation, hematopoietic cell lineage, antigen processing and presentation, Th17 cell differentiation, and Staphylococcus aureus infection (Figure 2B). The top six terms associated with the up- and down-regulated DEGs, enriched in BP, CC and MF as well as various KEGG pathways are shown in Tables 2 and 3, respectively. The details of the enriched DEGs in BP category of GO analysis and KEGG pathways are shown in Supplementary Tables 1 and 2, respectively.

Table 2 GO Terms Enrichment Analysis of DEGs (Top 6 of Each)

Table 3 KEGG Pathway Enrichment Analysis of DEGs (Top 6 of Each)

PPI Network Construction and Module Analysis

The interactions between DEGs, which consisted of 439 nodes and 2430 edges (Figure 3A), were downloaded from the STRING database and visualized using Cytoscape. Four significant modules (Figure 3B) were identified by cluster analyzing the PPI network in MCODE based on the aforementioned criteria. KEGG pathway enrichment analysis was then performed for the DEGs in the four modules, and the results are listed in Table 4. The DEGs in module 1 were significantly enriched in S. aureus infections. The DEGs in module 2 were enriched in chemokine signaling pathways, viral protein interactions with cytokines, and cytokine receptors. The DEGs in module 3 were mainly enriched in ubiquitin-mediated proteolysis. Finally, the DEGs in module 4 were enriched in Th1 and Th2 cell differentiation, Th17 cell differentiation, primary immunodeficiencies, and T cell receptor signaling pathways.

Table 4 KEGG Pathway Enrichment Analysis of the DEGs in the Four Modules (Top 6 of Each)

Figure 3 PPI network analysis. (A) The PPI network consists 439 nodes and 2430 edges were visualized in Cytoscape. Red notes the upregulated genes and green notes the downregulated genes. The size of each node is positive correlated with log2fold change (FC). The width of each edge is positively correlated with the combine score. (B) It shows the most significant modules in this PPI network. PPI, protein-protein interaction.

Hub Gene Selection

Hub genes were identified by CytoHubba in the present study. The top 10 hub genes, which were selected based on the 3 most commonly used classification methods in cytoHubba, are displayed in Table 5. By overlapping the first 15 genes, 7 central genes (ARG1, CKAP4, C3AR1, ELANE, GGH, ORM1, and QSOX1) were consequently identified as shown in Figure 4A. After being marked on selected modules, we observed that the hub genes were upregulated and gathered in module 1 (Figure 4B).

Table 5 List of the Top 10 Hub Genes Selected by MCC, MNC, and Degree Methods in cytoHubba

Figure 4 Identification of the hub genes. (A) Hub genes were identified by the overlapping of the three methods in cytoHubba. (B) From the module analysis we find that seven of the hub genes upregulated and gathered in module 1 in which showed by shape of square.

Discussion

Severe burn injuries usually result in distributive shock,¹⁹ an abnormal physical state that is marked by capillary leakage of fluid from the intravascular to the interstitial space, and is commonly accompanied by multiple host responses.²⁰ In the initial stages following a severe burn injury, burn shock presents complex pathophysiological changes under the influence of chemical transmitters, injury factors, and toxins due to the effects of injury, stress, and acute phase and immune-inflammatory reactions.^21,22 Severe burn injury is often accompanied by smoke inhalation injury, and there were 12 of 30 patients have inhalation injury in this study. In most cases, inhalation injury aggravates burn shock by increasing vascular permeability, and edema formation.²³ Burn shock is the earliest serious systemic complication of severe burns and could lead to multiple organ failure, especially in instances of delayed resuscitation for severe burns which not only has serious consequences in the early stages of the injury, but also has a significant impact on the entire course of the disease.²⁴ However, sufficient fluid resuscitation alone does not ameliorate the effects of burn shock.²⁵ Therefore, the exact molecular mechanism underlying the changes in host response during burn shock needs to be understood.

In the present study, the gene expression data of the GSE77791 dataset were obtained from the GEO database and bioinformatically analyzed to explore the key pathways and genes involved in the host response to burn shock. The GSE77791 microarray dataset contained 117 samples, including 13 samples from healthy volunteers, 30 samples from patients without any treatment administration, and the remaining 74 samples from 1 day, 5 days and 7 days after hydrocortisone or NaCl 0.9% administration. The dataset was utilized for assessing the hydrocortisone-induced transcriptional modulation of the immune response in severe burn shock without a deep microarray analysis in a previous study, which found that hydrocortisone administration may worsen the immunosuppression associated with severe burn shock.¹⁸ The present study, to our knowledge, is the first to extract data from the 13 healthy volunteer samples and the 30 burn shock patient samples before any treatment administration for further analysis and we subsequently identified the key genes associated with the host response to severe burn shock.

A total of 1059 DEGs, included 508 downregulated genes and 551 upregulated genes, were identified. The GO analysis annotates each DEG and enriches the DEGs with the same attribute into one term. As shown in Figure 2A, the DEGs enriched terms in the BP category were neutrophil activation, neutrophil activation involved in immune response, neutrophil degranulation, neutrophil-mediated immunity, T cell activation, T cell differentiation, lymphocyte differentiation, and T cell receptor signaling pathway. Interestingly, Table 2 shows that the upregulated genes were mostly enriched in the GO terms of neutrophil activation, neutrophil activation involved in immune response, neutrophil degranulation, and neutrophil-mediated immunity, which are associated with neutrophil activation, while the downregulated genes were mostly enriched in terms of T cell activation, T cell differentiation, lymphocyte differentiation, and T cell receptor signaling pathway, all of which are associated with T cell activation. There were many DEGs enriched in neutrophil activation and T cell activation and Supplementary Table 1 was added to show the specific DEGs. The key genes (ARG1, CKAP4, C3AR1, ELANE, GGH, ORM1, and QSOX1) identified by cytoHubba were all enriched in terms associated with neutrophil activation. Further in vivo and in vitro studies are warranted to determine the role of these key genes in burn shock. These GO analysis results are consistent with the postburn immune response, in which the neutrophil is over activated and T cell activation is suppressed. Neutrophils are immediately activated and migrate to the injury site to remove debris and pathogens as well as promote wound healing after burn injury.^26,27 However, they also accumulate in other remote organs, such as the lung and small intestine and eventually lead to the injury of these organs.²⁸ Overactivated neutrophils can secrete free radicals and inflammatory mediators such as interleukin (IL)-1, tumor necrosis factor (TNF), IL-6, and IL-8 after severe burn injury and result in systemic inflammatory response syndrome and multiple organ failure.²⁹ Meanwhile, T cell proliferation is also suppressed in severe burn injuries.³⁰ According to the p value in Figure 2B, the top three pathways with the highest reliability were Th1 and Th2 cell differentiation, Th17 cell differentiation, and hematopoietic cell lineage, in which the downregulated DEGs were enriched. There were many DEGs enriched in the various pathways and Supplementary Table 2 was added to show the specific DEGs. The signal transduction in T cell activation in burn shock was suppressed according to KEGG pathway analysis.

After severe burn injury, immune functions are severely compromised, including the decline of neutrophil killing of invading pathogens and macrophage antigen presentation as well as suppression of T cell activation.⁵ Together, these events lead to the development of immune suppression during the burn stage.⁵ According to the results of GO function analysis and KEGG pathway enrichment analysis, the most downregulated DEGs were enriched in T cell activation. Neutrophils are activated in burn shock, and the granules released by overactivated neutrophils can suppress T cell activation.³¹ Therefore, immunity may be impaired by the suppression of T cell activation.³¹ Our findings are consistent with the pathophysiology of burn shock in terms of neutrophil overactivation and lymphocyte suppression.^32,33 This can explain the increased risk of severe systemic infection and multiple organ failure in the severe burn patients with an unstable shock period.

Four significant modules were selected from the PPI network using MCODE in Cytoscape. Next, KEGG pathway enrichment analysis was performed for the DEGs in these modules. As shown in Table 4, the DEGs in modules 1, 2, and 4 were mostly enriched in the pathways associated with immune response, and the DEGs in module 3 were mostly enriched in the pathways associated with protein metabolism. This is consistent with the hypermetabolic state in severe burn injury.³⁴ Jeschke suggested that the continued release of TNF, IL-6, and other inflammatory mediators could further contribute to the hypermetabolic state of severe burn injury patients.³⁵ Furthermore, immune dysregulation could also affect the metabolic state of the body.

Seven hub genes (ARG1, CKAP4, C3AR1, ELANE, GGH, ORM1, and QSOX1) were identified using cytoHubba and gathered in module 1 with a high degree. ARG1, which is contained in cytoplasmic azurophil granules and the gelatinase grains of neutrophils,^36,37 hydrolyzes L-arginine to L-ornithine and urea in the cytosol. ARG1-induced depletion of arginine, suppression of T lymphocyte proliferation, and secretion of cytokines via down-regulated expression of the CD3 zeta chain in T cells³¹ have emerged as fundamental mechanisms underlying inflammation-associated immunosuppression.³⁸ In addition, ARG1 competes with inducible nitric oxide synthase for substrate L-arginine to inhibit the production of nitric oxide, which can also suppress the inflammatory response.³⁹ C3AR1 is a G protein-coupled receptor protein that is predominantly distributed on neutrophils, eosinophils, basophils, mast cells, and monocytes/macrophages.⁴⁰ Previous studies have described C3a as a proinflammatory mediator in some chronic inflammatory responses.^41,42 However, several studies have indicated that C3AR1 exerts a potent anti-inflammatory response with a high expression of C3AR1 on neutrophils by inhibiting the mobilization of neutrophils from the bone marrow to circulation following acute injury.^43,44 Neutrophil elastase (NE) is coded by ELANE and secreted by activated neutrophils to kill bacteria and promotes the secretion of pro-inflammatory mediators, such as TNFα and IL-1β. Numerous studies have indicated that NE is a marker of neutrophil activation.⁴⁵ ORM1, also known as alpha–1–acid glycoprotein (AGP1), encodes an acute-phase protein in activated neutrophil that mediates neutrophil migration failure and inhibition of neutrophil activation in the acute inflammatory response.^46,47 The latest study indicates that ORM1 serves as a potential diagnostic marker for sepsis.⁴⁸

CKAP4 is a transmembrane protein residing predominantly in the endoplasmic reticulum of eukaryotes that can promote cellular proliferation through the phosphoinositide 3-kinase (PI3K)/protein kinase B (AKT) pathway.⁴⁹ Research has shown that PI3K plays a vital role in the regulation of neutrophil activation through AKT activation.^50,51 Therefore, CKAP4 may activate neutrophils through the PI3K/AKT pathway in burn shock patients. GGH is a lysosomal enzyme involved in the metabolism of folates and the disposition of methotrexate that plays a role as a tumor marker in breast and lung cancer.^52,53 GGH also catalyzes the hydrolysis of folylpoly-gamma-glutamates to generate glutamate, which is a substrate of glutamine.⁵⁴ Glutamine, a substrate of enterocytes, plays a vital role in the normal immunological function and the structure of the gut.⁵⁵ The intense stress and persistent inflammation caused by severe burn injury may result in intestinal mucosal barrier damage which is sensitive to infection.⁵⁶ Wang and colleagues found that glutamine improved tissue perfusion and increased energy synthesis in enterocytes, thus alleviating the intestinal injury after burn injury.⁵⁷ QSOX1 is mainly expressed in several tumor cells and can promote tumor growth and metastasis.^58,59 A recently study reported that the peripheral blood level of QSOX1 related to neutrophil infiltration is a novel independent predictor of left ventricular dysfunction.⁶⁰

To date, the clinical value of ARG1, CKAP4, C3AR1, ELANE, GGH, ORM1, and QSOX1 in burn shock has not been documented. From the above description, we know that ARG1, C3AR1, ELANE, and ORM1 are released from activated neutrophils and have different effects on immune response. Neutrophil activation is considered the most dominant host response to burn shock; thus, it is crucial to conduct further studies to explore the diagnostic or therapeutic value of these four genes in burn shock. Although CKAP4, GGH, and QSOX1 play crucial roles in tumor proliferation and metastasis, these genes were also significantly upregulated in neutrophil in burn shock and had close interactions with ARG1, C3AR1, ELANE, and ORM1 in the PPI network constructed using the DEGs. Furthermore, an increasing amount of research indicates that CKAP4, GGH, and, QSOX1 are associated with the immune response of neutrophils. Therefore, we speculate these three genes play crucial effects on the immune response to severe burn injury, but more studies are needed to confirm this hypothesis. Our study revealed the host response to burn shock and identified the key genes that provide a new direction for further in vivo and in vitro studies to explore the mechanism of occurrence and development of burn shock.

There are no precise methods to identify hub genes. Some researchers select the higher-degree nodes in the PPI network as hub genes,⁶¹ while others identify the hub genes via cytoHubba in Cytoscape using a single method.⁶² In addition, some researchers believe that the overlapping results from several methods in cytoHubba is a more reasonable method for detecting hub genes.⁶³ In this study, we overlapped the top 10 genes selected by the MCC, MNC, and degree methods in cytoHubba to identify the hub genes.

Several previous studies have used the bioinformatic analysis of microarray datasets to explore gene expression after burn injury. Xu⁶⁴ and Wu⁶⁵ reanalyzed the GSE19743 dataset and determined the time course change in gene expression after burn injury. The time points of sample collection in the GSE19743 dataset were separated into two groups: early-stage for < 11 days postburn and middle-stage for 11–49 days after thermal injury. Bioinformatic analysis of the GSE19743 dataset was preformed to reveal the gene changes over the whole burn stage which ranged from 0 to 49 days after thermal injury. However, the 30 samples in GSE77791 were collected within 72 h after burn injury, and the analysis of these samples was utilized to identify gene changes in the shock stage after burn injury. The GSE77791 dataset, to our knowledge, contains unique burn shock data in the GEO database and this study is the first to use GSE77791 to explore the different gene expression levels during the burn shock stage after burn injury.

Compared to previous studies, the present study has several advantages. First, our study successfully identified novel hub genes associated with severe burn shock. Second, we further analyzed and presented the functions and pathway enrichment of the main DEGs. Our findings can help to elucidate the molecular mechanism underlying the changes in host response during burn shock as well as provide potential targets for early detection and treatment of burn shock. However, some limitations are also present in our study. We only made use of a single dataset from the GEO database, which is the only one we could find that was associated with burn shock. We will repeat the experiment in the future should other datasets, associated with burn shock, become available. Moreover, we did not verify the expression of these hub genes in patients with severe burn injury because samples from patients with severe burn injury are limited. However, we will verify the expression levels of the hub genes in severe burn injuries and perform further studies on the mechanisms of these genes in the occurrence of severe burn injury shock once the samples are available.

Conclusions

In conclusion, following the bioinformatic analyses of DEGs, GO terms, KEGG pathway enrichment, and the PPI network, our study found that the immune response of the host is altered during burn shock, which may be mediated by specific hub genes. The present study provides a new direction for future research to investigate the underlying mechanism of burn shock occurrence and development.