Screening and identification of the core immune‐related genes and immune cell infiltration in severe burns and sepsis

Abstract Severe burns often have a high mortality rate due to sepsis, but the genetic and immune crosstalk between them remains unclear. In the present study, the GSE77791 and GSE95233 datasets were analysed to identify immune‐related differentially expressed genes (DEGs) involved in disease progression in both burns and sepsis. Subsequently, weighted gene coexpression network analysis (WGCNA), gene enrichment analysis, protein–protein interaction (PPI) network construction, immune cell infiltration analysis, core gene identification, coexpression network analysis and clinical correlation analysis were performed. A total of 282 common DEGs associated with burns and sepsis were identified. Kyoto Encyclopedia of Genes and Genomes pathway analysis identified the following enriched pathways in burns and sepsis: metabolic pathways; complement and coagulation cascades; legionellosis; starch and sucrose metabolism; and ferroptosis. Finally, six core DEGs were identified, namely, IL10, RETN, THBS1, FGF13, LCN2 and MMP9. Correlation analysis showed that some core DEGs were significantly associated with simultaneous dysregulation of immune cells. Of these, RETN upregulation was associated with a worse prognosis. The immune‐related genes and dysregulated immune cells in severe burns and sepsis provide potential research directions for diagnosis and treatment.


| INTRODUC TI ON
Burns are the fourth most common trauma worldwide, 1 and they often lead to sepsis and high mortality. Sepsis is an excessive immune response induced by severe infection that can trigger systemic multiple organ dysfunction and ultimately endanger the patient's life. 2 After sepsis occurs in the body, the production of proinflammatory cytokines increases, the number of lymphocytes decreases, and immune function is suppressed, 3 leading to dysfunction of blood vessels, nerves, hormones, metabolism and coagulation. 4,5 F I G U R E 1 Research design flow chart. In this study, we used WGCNA to explore the correlation between module signature genes and the occurrence of sepsis and burns. At the same time, differentially expressed genes (DEGs) of burns and sepsis were obtained with GEO2R. The genes of the overlapping part between the module genes and DEGs were used for further analysis. Subsequently, core gene identification, immune infiltration analysis and clinical relevance studies were performed to explore the potential link between burns and sepsis, thus providing new references for the diagnosis and treatment of patients with secondary sepsis after burns.
In severe burns, the adaptive immune function represented by T cells is usually inhibited, and the inflammatory response represented by neutrophils and dendritic cells is activated. 6 The disordered immune state causes the body to easily induce severe infection.

Inflammatory factors, cytokines and immune cells have significant
prognostic value in severe burns. 7 Numerous studies have shown that patients with burn-induced sepsis are characterized by opportunistic infections and immunosuppression. 3,[8][9][10][11] The mechanisms of sepsisinduced immunosuppression include apoptosis depletion of immune cells, increased expression of negative costimulatory molecules, increased expression of regulatory T cells (Tregs) and programmed cell death of CD4 + T cells. 3 In the early stage of sepsis, macrophages and neutrophils initiate a proinflammatory cascade, which is mainly achieved by the release of cytokines, such as tumour necrosis factor (TNFα) and interleukin (IL1). 11 In addition, anti-inflammatory mediators, such as IL4 and IL10, are also released during sepsis, a process that not only suppresses the proinflammatory cascade but also leads to a state of relative immunosuppression in patients with sepsis. 12 In the late stage of sepsis, excessive apoptosis of macrophages leads to immune dysfunction and organ damage, 13,14 which increases the probability of secondary host infection and eventually leads to patient death.
Currently, there are no tests or diagnostic criteria to identify sepsis in critically ill patients. Although many patients are affected by burn sepsis, the treatment of burn sepsis still faces many difficulties due to the lack of clear diagnosis and treatment guidelines. To overcome the delay of diagnosis and treatment, we used weighted gene coexpression network analysis (WGCNA) to explore the correlation between module signature genes (highly interconnected gene set) and the occurrence of sepsis and burns. At the same time, DEGs of burns and sepsis were obtained with GEO2R. The genes of the overlapping module genes and DEGs were used for further analysis. Subsequently, core gene identification, immune infiltration analysis and clinical relevance studies were performed to explore the potential link between burns and sepsis, thus providing new references for the diagnosis and treatment of patients with secondary sepsis after burns.
F I G U R E 2 (A,D) Clustering dendrogram of samples based on their Euclidean distance in GSE95233 and GSE77791. (B,E) Determination of soft-thresholding power for GSE95233 and GSE77791. (C,F) Heatmap of the correlation between module eigengenes and the occurrence of sepsis and burns. Modules with different colours represent different gene modules, and the numbers in the modules represent the correlation between the module and the phenotype.

| Raw data
The datasets of GSE95233 15 and GSE77791 15 (GPL570 platform, Affymetrix Human Genome U133 Plus 2.0 Array) were downloaded from the Gene Expression Omnibus (GEO) database (http://www. ncbi.nlm.nih.gov/geo). 16 The GSE95233 dataset contains 51 sepsis patients and 22 healthy volunteers. We selected blood samples without any treatment on the first day for comparative analysis.
The GSE77791 dataset contains 30 severely burned patients (a total burns surface area (TBSA) range from 30% to 98%) and 13 healthy volunteers. Similarly, we selected pre-treatment blood samples for subsequent analysis. Some research methods in this paper have been applied in our previously published manuscripts. 17

| WGCNA network construction and module identification
Coexpression networks were constructed to identify coexpression modules via the WGCNA package in R. 18 First, the samples were subjected to hierarchical clustering analysis to assess whether there were any significant outliers; second, the soft threshold power β was calculated by the picksoftworthold algorithm in R to construct a biologically meaningful scale-free network; third, the topological overlap matrix establishes network interconnectivity and identifies gene modules based on a dynamic tree-cutting algorithm; fourth, modules were linked to clinical features by calculating gene significance and module membership. The genes involved in the corresponding modules were used for subsequent analysis. Finally, the eigengene network was visualized.

| Identification of common DEGs
GEO2R (www.ncbi.nlm.nih.gov/geo/ge2r) is an online tool based on Limma package for analysing DEGs. 19 'Adjusted p-values < 0.01 and |logFC| ≥ 1' were defined as the thresholds for the screening of DEGs. We extracted the genes in the WGCNA modules most positively associated with sepsis and burns and cross-analysed these two gene lists to obtain common DEGs.

| Enrichment analyses of common DEGs
To further understand which biological functions these common DEGs are involved in, we performed gene enrichment analysis using the KOBAS 3.0 database, including gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. 20 The species was set as Homo sapiens. Adjusted p-values < 0.05 was considered significant.
We downloaded the currently known immune-related genes from the Immunology Database and Analysis Portal database (ImmPort; https://www.immpo rt.org) (Table S1). Then, Venn diagrams were intersected to obtain core immune-related DEGs in sepsis and burns.

| Immune infiltration analyses
We assessed the corresponding immune cell infiltration status of these two datasets by Immune Cell Abundance Identifier (ImmuCellAI) (http://bioin fo.life.hust.edu.cn/ImmuC ellAI/ #!/), 25 which provides comprehensive predictions of immune cell abundance by assessing the abundance of 24 immune cell types. We further compared the differences between the disease group and the control group by t-test. Finally, we explored the connection between immune cells and core immune-related DEGs by Spearman correlation analysis.

F I G U R E 4 (A) Enrichment result of common DEGs GO term; (B) enrichment result of common differentially expressed genes (DEGs)
KEGG pathway. Adjusted p-value < 0.05 was considered significant. The ordinate represents the enriched term, and the abscissa represents the proportion of genes involved in the term. The size of the dots represents the number of genes, and the colour of the dots represents the p value.
F I G U R E 5 (A) PPI network constructed using the STRING database. The modules in the box are the core genes screened out by the MCODE algorithm. (B) The core module genes were intersected with immune-related genes from the database to obtain six core immunerelated differentially expressed genes (DEGs). (C) Core immune-related DEGs and their co-expression genes were analysed via GeneMANIA.

TA B L E 1
The details of the core immune-related DEGs.
No. This protein plays a key role in extracellular matrix protein degradation in normal physiology and in a variety of pathologies, including those associated with inflammation

| Clinical significance of core immune-related DEGs in sepsis and burns
To verify our results, the expression of core immune-related DEGs was extracted from GSE95233, GSE77791, GSE19743 and GSE57065, and the difference between disease group and normal group was analysed by t-test. Receiver operating characteristic curve analysis was used to examine the diagnostic value of core immune-related DEGs.
Since there is survival information for the burns and sepsis cohort in the GSE19743 and GSE95233, we compared the differences in core immune-related DEGs between the survival and non-survival groups.

| WGCNA network construction and module identification
The flowchart of the present study is shown in Figure 1. WGCNA was used to explore core modules significantly associated with disease. Potential outlier samples marked with red bars in the <outlierC> row were eliminated in the GSE95233 dataset, and no samples were eliminated in the GSE77791 dataset (Figure 2A,D).
In the WGCNA method, b = 16 was the optimal soft threshold for the GSE95233 dataset ( Figure 2B), and b = 12 was the optimal soft threshold for the GSE77791 dataset ( Figure 2E). Ten and seven modules were identified in the GSE95233 and GSE77791 datasets, respectively. Correlation analysis showed that the turquoise module was strongly positively correlated with sepsis (r = 0.735), while the blue module was strongly negatively correlated with sepsis (r = −0.729) ( Figure 2C). For burns, the turquoise module showed a strong positive correlation (r = 0.896), while the blue module had a strong negative correlation (r = −0.879) ( Figure 2F).

| Identification of common DEGs
According to the GEO2R online analysis results, 1546 DEGs and 1752 DEGs were identified in the GSE95233 and GSE77791 datasets, respectively ( Figure 3A,B). In addition, we extracted genes from the most prominent modules for burns and sepsis. As shown in the flow chart, we used WGCNA to obtain the module signature genes F I G U R E 6 Differential expression of six core immune-related genes in GSE95233 and GSE77791. Differences in expression of the six core genes in the original data were verified by t-test. The numbers on the ordinate represent probe values. *p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001.
closely related to diseases. The module signature genes were then intersected with the DEGs identified by GEO2R to further obtain the common DEGs. Through Venn diagram calculation combined with WGCNA, we obtained 282 common DEGs between sepsis and burns ( Figure 3C and Table S2).

| Analysis of the functional characteristics
GO analysis showed that the common DEGs were significantly enriched in neutrophil degranulation, protein binding, extracellular exosome, plasma membrane and cytosol ( Figure 4A). KEGG pathway analysis showed that the DEGs were mainly concentrated in the following pathways: metabolic pathways; complement and coagulation cascades; legionellosis; starch and sucrose metabolism; and ferroptosis ( Figure 4B).

| Protein-protein interaction (PPI) network construction and analysis of core immunerelated DEGs
Through the STRING database, we obtained a PPI network with 205 nodes and 348 interaction pairs ( Figure 5A). Subsequently, we obtained a core module containing 10 DEGs and 30 interaction pairs using the MCODE plugin ( Figure 5A). After taking the intersection of the Venn diagrams, we identified six core immunerelated DEGs, namely, IL10, RETN, THBS1, FGF13, LCN2 and MMP9 ( Figure 5B and Table 1). Coexpression analysis revealed that these genes showed a complex PPI network with physical interactions of 57.02%, coexpression of 29.44%, prediction of 10.34% and colocalization of 3.2% via the GeneMANIA database ( Figure 5C). Finally, a t test confirmed that the expression of the six core immune-related DEGs was significantly upregulated in the disease group ( Figure 6).

| Immune infiltration analyses
Evaluation of the immune cell composition of sepsis and burn pa-  Figure 9A). In burn patients, FGF13 was strongly positively correlated with CD8_T cells but significantly negatively correlated with neutrophils ( Figure 9B).

| Clinical relevance of core immune-related DEGs in sepsis and burns
According to a t test, the six immune-related DEGs conformed to the expression trends in the aforementioned sepsis and burn datasets ( Figure 10A,B). In addition, we analysed the immune cell infiltration profiles of the GSE57065 and GSE19743 datasets, which were consistent with previous results (Figure 11A,B). ROC analysis showed that these genes had good diagnostic value for burns and sepsis ( Figure 12A,C). Because the GSE95233 and GSE19743 datasets contain survival information, we compared the differences in core immune-related DEGs between survival and nonsurvival groups to explore the prognostic potential of these immune-related DEGs for sepsis and burn patients. The results showed that the expression levels of RETN and THBS1 were lower in the survival group, which indicated that high expression of RETN and THBS1 was associated with a poor prognosis in the sepsis group. In the burn group, however, the expression levels of RETN, LCN2 and MMP9 were lower in the survival group, which suggested that high expression of RETN, LCN2 and MMP9 predicted a poorer prognosis ( Figure 12B,D).

F I G U R E 8 (A) Comparison of immune cell fractions between sepsis patients and healthy controls. (B) Comparison of immune cell fractions between burns patients and healthy controls. The horizontal axis represents different immune cells, and the vertical axis
represents the proportion of immune cells. t test was used to compare the disease group and the control group. *p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001.

| DISCUSS ION
Burn-induced sepsis can affect both the innate and adaptive immune systems of the body. In burn sepsis patients, the immune function of NK cells and neutrophils is impaired, [26][27][28] while the phagocytic ability of macrophages is weakened. 29 Evaluation of the immune cell composition of sepsis and burn patients identified significant differences in immune cell profiles between the diseased and control groups. In sepsis and burn patients, monocytes, macrophages, neutrophils and NKT cells were increased, whereas NK cells, CD4 + T cells and CD8 + T cells were decreased.
Due to their widespread presence in various tissues, macrophages play roles in phagocytosis, sterilization, antigen presentation and inflammatory cytokine secretion in various stages of sepsis. [30][31][32] After the body suffers from burns, macrophages are overactivated into proinflammatory (M1) macrophages in the early stage of sepsis, which release inflammatory factors (IL1, IL6 and TNFα). Decreased phagocytic capacity has also been shown after severe burns. 33,34 Studies have shown that the release of IL4 and IL10 significantly inhibits the antigen-presenting ability of macrophages as well as inhibits the bactericidal activity of NK cells and neutrophils. [35][36][37] Therefore, there is a correlation between elevated IL10 levels and infectious events. 38,39 Interestingly, the present study identified six In addition, the respective proportions of CD4 + T cells and CD8 + T cells decreased after severe burns. 45 In the event of sepsis, naive

| CON CLUS IONS
In conclusion, the present study identified core immune-related DEGs between severe burns and sepsis, including IL10, RETN, THBS1, FGF13, LCN2 and MMP9. Among these DEGs, increased expression of RETN predicts poor prognosis. Furthermore, there is some association between these core genes and dysregulated immune cells. Collectively, these dysregulated core genes and immune cells provide potential research avenues for severe burns and sepsis.

CO N FLI C T O F I NTER E S T S TATEM ENT
The authors declare no conflict of interests.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data analyzed in this article comes from Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo). The accession number can be found in the section 2 of the article.
F I G U R E 1 2 (A,C) ROC curve analysis of core immune-related DEGs in sepsis and burns. Area under the curve (AUC) varies from 0.5 to 1.0. The closer the AUC gets to 1, the more accurate the diagnostic test is (0.5 -0.7: low-moderate accuracy, 0.7 -0.9: moderate accuracy, 0.9 -1: high accuracy). (B, D) Differences in expression of core immune-related DEGs between survival and non-survival of sepsis and burns patients. t-test was used for analysis. *p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001.