Identification of hub genes and prediction of the ceRNA network in adult sepsis

Datasets and sample selection

Gene expression datasets related to sepsis were searched via the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/; Barrett et al., 2013). The present study systematically retrieved microarray datasets based on the following terms: ‘Sepsis’, ‘Adult’, ‘Homo sapiens’, and ‘Microarray’. Datasets were acquired based on the following eligibility criteria: (1) containing ≥5 sepsis patients and ≥5 healthy controls; and (2) containing raw data or gene expression by array accessible through GEO. Seven GEO datasets (GSE137340, GSE69063, GSE69528, GSE54514, GSE57065, GSE95233, and GSE28750) were included (Table 1). The series matrix file(s) or raw data were downloaded from the GEO website, along with the corresponding annotation documents from GEO or Bioconductor (https://bioconductor.org/).

Robust rank aggregation analysis

Statistically changed genes and differentially expressed genes (DEGs) were screened using the robust rank aggregation (RRA) method, which is an unbiased and effective approach for integrating data from diverse chip platforms using a comprehensive ranking list algorithm (Jia & Zhai, 2019; Yan et al., 2018). After identifying up-ranked and down-ranked gene lists of each dataset based on the fold-change in expression levels between patients with sepsis and healthy controls, the R ‘Robust Rank Aggregation’ package was used to integrate all of the ranked gene lists. The adjusted P-value in the analysis informed the ranking of each gene in the final list. Finally, the ‘limma’ (v.3.54.2) and ‘RobustRankAggreg’ (v.1.2.1) R packages were used to screen statistically changed genes (adjusted (adj.) P < 0.05) and DEGs (adj. P < 0.05 and —log fold-change— > 1.2).

Function and pathway enrichment analyses

To explore the impact of these genes on the development of sepsis, gene set enrichment analysis (GSEA), Gene Ontology (GO) analysis, and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were performed via the ‘clusterProfiler’ (v. 4.6.2) R package (Wu et al., 2021). The list ‘c2.cp.reactome.v2022.1.Hs.symbols.gmt’ from the Molecular Signatures Database (https://www.gsea-msigdb.org/gsea/msigdb/) was used as a reference dataset (Liberzon et al., 2015). Ridgeline plots and normalized enrichment scores (NES) were used to present the results of the GSEA analysis. GO contained the biological process (BP), cellular component (CC), and molecular function (MF) categories. Adjusted P < 0.05 was set as the threshold criterion.

Evaluation of immune cell infiltration in peripheral blood

Based on the gene expression matrix, cell-type identification via the estimation of relative subsets of RNA transcripts (CIBERSORT) was used to analyze the percentage of immune cells in peripheral blood. The expression matrix was standardized by quantile normalization, and the operational mode was configured as batch mode, concerning the known set LM22 (Newman et al., 2015). P < 0.05 was considered to indicate a statistically significant difference. The Mann–Whitney U test was used to evaluate significant differences in infiltrating immune cell types between individuals with sepsis and healthy controls. The ‘ggpubr’ (v. 0.6.0) package was used to generate boxplots showing immune cell infiltration differentiation.

Weighted gene co-expression network analysis

Weighted gene co-expression network analysis (WGCNA) is a data exploration tool that can identify key gene modules using unsupervised clustering without a priori-defined gene set (Liu et al., 2017). The R package ‘WGCNA’ (v. 1.72.1) was used in the present study to construct a co-expression network based on all expressed genes in GSE69528. Genes were divided into different modules using the topological overlap matrix-based dissimilarity measure. Key modules were identified by setting the soft-thresholding power (β) as 7 (R² = 0.85), cut height as 0.25, minimal module size as 30, and threshold as 0.2. Genes with similar expression profiles were aggregated into the same cluster via a dynamic tree cut algorithm. The gene modules most strongly associated with sepsis were selected through the computation of their association with phenotype. A gene set was obtained via Venn analysis of DEGs from the RRA analysis and genes in the module for subsequent evaluation.

Protein–protein interaction network analysis

The STRING database (https://cn.string-db.org/) is an online resource for accessing established proteins and predicting protein–protein interactions (PPIs), encompassing both direct physical associations between proteins and their indirect functional connections (Szklarczyk et al., 2021). In the present study, a PPI network was constructed using the STRING online tool with default filter conditions (combined score > 0.4). Subsequently, Cytoscape 3.9.0 (https://cytoscape.org/) software was used to optimize the PPI network for enhanced visualization (Majeed & Mukhtar, 2023). Plug-in minimal common oncology data elements (MCODE; https://apps.cytoscape.org/apps/mcode) was used to screen important gene clusters (Set 1). Plug-in cytoHubba-maximal clique centrality (MCC; https://apps.cytoscape.org/apps/cytohubba) with default parameters was also used to identify key genes (Set 2) in this network (Chin et al., 2014). Plug-in iRegulon (v.1.3; http://apps.cytoscape.org/apps/iRegulon) was used with default parameters to identify the transcription factor (TF) of Set 2 (Janky et al., 2014).

Peripheral blood samples from patients with sepsis and health controls

Patients were included in the study if they were 20 years of age or older and had been diagnosed with sepsis according to the Third International Consensus Definitions for Sepsis and Septic Shock (Sepsis-3; Singer et al., 2016). The present study included nine septic patients and nine healthy controls. Healthy controls composed of laboratory and hospital employees were recruited, whereas patients with sepsis were admitted to the Intensive Care Unit (ICU) of The Third Affiliated Hospital of Southern Medical University. Peripheral blood samples were collected through venipuncture. The period of recruitment for this study was from November 1 to December 15, 2023. Written informed consent was obtained from the patient or the patient’s family. The research protocol and use of samples were both approved and formally authorized by the Clinical Trial Ethics Committee of The Third Affiliated Hospital of Southern Medical University (approval no. 2023073; Guangzhou, China) on October 30, 2023.

Animal care and housing

SPF-grade C57BL/6 healthy mice were purchased from Southern Medical University. Upon arrival, the animals were acclimatized for five days prior to the start of the experiments. Mice were provided ordinary Specific Pathogen Free (SPF) food and water freely throughout the experiment. All mice were housed in units of equal size and material and were exposed to 12 h light-12 h darkness-cycles. The feeding temperature was 22(± 2) °C and the humidity was 40∼70%. All experiments were conducted using male adult mice (6∼8 weeks old, weighing 20∼22 g) during the light cycle from 10:00 am to 5:00 pm. All surviving animals at the conclusion of the experiment were to be euthanized by rapid neck dislocation, but this was unnecessary as all mice died after blood collection.

Animal model of cecal ligation and puncture

Animal experiments were approved by the Animal Welfare and Ethics Committee of Southern Medical University, Guangzhou, China (approval no. L2018235) and were conducted in accordance with the Guide for the Care and Use of Laboratory Animals. Sixteen SPF-grade male C57BL/6 mice were randomly divided into two groups, cecal ligation and puncture (CLP) and sham, with eight mice in each group. CLP modeling of sepsis was performed according to a previously established protocol (Li et al., 2022). In brief, C57BL/6 mice were anesthetized with isoflurane at 5% for induction and 2.5% for maintenance. After confirming loss of consciousness, a 2-cm midline incision was made on the anterior abdomen to expose the cecum, and it was then ligated at a site 1 cm from the end of the cecum. The cecum was twice punctured with an 18-gauge needle between the ligation site and the end of the cecum, and a tiny amount of cecal content was then extruded. Following the cecum’s return to the abdomen, wound clips were used to seal the incision. With the exception of the cecum puncture and ligation, the sham group received the same care as the CLP group. During the operation, if any signs of pain were observed, appropriate care was immediately provided and the necessary amount of anesthesia was adjusted accordingly. For fluid replacement, each animal received a subcutaneous injection of 1 mL of sterile saline. Afterwards, the mice were placed on thermal blankets to facilitate their recovery from anesthesia. At twenty-four hours after surgery, all surviving mice were again anaesthetized and peripheral blood was collected from the inferior vena cava. The samples were immediately placed in an ice box for storage and were processed for the subsequent experiments within one hour.

RNA extraction and reverse transcription-quantitative polymerase chain reaction

RNA was extracted from peripheral blood samples of patients with sepsis and healthy controls using a GeneJET Stabilized and Fresh Whole Blood RNA Kit (cat. no. K0871; Thermo Fisher Scientific, Inc., Waltham, MA, USA), and single-stranded cDNA was synthesized from 1 µg of RNA using a ReverTra Ace qPCR RT Kit (cat. no. FSQ-101; Toyobo Life Science, Osaka, Japan). The reaction condition was as follows: 37 °C for 15 min and 98 °C for 5 min. The primers of target genes were acquired from PrimerBank (https://pga.mgh.harvard.edu/primerbank/) and were synthesized at BGI Genomics (Table S1). Quantification of gene expression was conducted by reverse transcription-quantitative polymerase chain reaction (RT-qPCR) using SYBR qPCR Mix (cat. no. QPK-201; Toyobo Life Science) on a 7500 Fast Real-Time PCR System (Applied Biosystems; Thermo Fisher Scientific, Inc.). The polymerase chain reaction (PCR) cycling parameters were as follows: initial denaturation at 95 °C for 3 min followed by 40 cycles of PCR reaction at 95 °C for 10 s, 60 °C for 10 s, and 72 °C for 30 s. The reaction system was set up as follows: SYBR Green Mix (5 µL), forward primer (0.2 µL), reverse primer (0.2 µL), cDNA template (0.1 µL), and ddH2O (4.5 µL), resulting in a total volume of 10 µL. PCR amplification was performed once for each sample, and the mRNA expression of target genes was normalized to 18S rRNA. Relative expression was determined using the 2^−ΔΔCT method (Zang et al., 2023): ΔΔCT = (CT_Target−CT_18srRNA)_Sample−(CT_Target−CT_18srRNA)_Control. Peripheral blood mononuclear cells were extracted from patients with sepsis and mice after CLP using the Peripheral Blood Mononuclear Cell Isolation Kit (cat. no. P8680; Beijing Solarbio Science & Technology Co., Ltd., Beijing, China) according to the manufacturer’s instructions, and subsequent steps were executed as described above.

Receiver operating characteristic analysis

To evaluate the clinical prognostic significance of the identified hub genes in sepsis, the GSE54514 dataset was used, which includes survivor and non-survivor data. Receiver operating characteristic (ROC) curves were plotted using the R ‘pROC’ (v. 1.18.4) package (Jiang et al., 2022). The area under the curve (AUC) was calculated for each hub gene to evaluate its performance. A set of guidelines was established to determine the accuracy of different diagnostic criteria, with each set of criteria having excellent (0.9 ≤ AUC < 1), good (0.8 ≤ AUC < 0.9), moderate (0.7 ≤ AUC < 0.8), poor (0.6 ≤ AUC < 0.7) or insufficient (0.5 ≤ AUC < 0.6) accuracy. If AUC > 0.8, the hub genes were regarded to have excellent sensitivity and specificity in differentiating between sepsis survivors and non-survivors.

Correction analysis

The ‘Immunity’ gene list was downloaded from the GeneCards database (https://www.genecards.org/), which is a searchable, integrative database that provides comprehensive, user-friendly information on all annotated and predicted human genes. Next, the ‘corrplot’ (v. 0.92) R package was used to visually represent the correlation matrix (Pei et al., 2020). P < 0.05 was considered to indicate a statistically significant difference.

CeRNA network construction

Four online microRNA (miRNA or miR) websites, namely miRWalk 3.0 (http://mirwalk.umm.uni-heidelberg.de/), microT v5 (http://diana.imis.athena-innovation.gr/DianaTools/index.php?r=MicroT_CDS/index), miRDB (https://mirdb.org/), and TargetScan Human 8.0 (https://www.targetscan.org/vert_80/), were used to predict target miRNAs (Agarwal et al., 2015; Chen & Wang, 2020; Dweep, Gretz & Sticht, 2014; Paraskevopoulou et al., 2013). The target miRNAs were selected from the intersection of ≥3 databases. Interactions between long noncoding RNAs (lncRNAs) and the selected miRNAs were predicted using StarBase v2.0 (https://rnasysu.com/encori/; Li et al., 2014). Cytoscape software was used to construct and visualize the interaction networks of mRNAs, miRNAs, and lncRNAs.

Statistical analysis

R version 4.2.0 (https://www.r-project.org/) was used to conduct statistical analyses. Paired t-tests were used to compare differences between two groups. Data are presented as the mean ± standard error of the mean. Statistical analysis and image construction were performed using GraphPad Prism 8.0.2 (GraphPad; Dotmatics, Boston, MA, USA). Figures were edited using Adobe Illustrator (AI CC2018) software (Adobe Systems, Inc., San Jose, CA, USA). Statistical significance was detected at the 0.05 level.

Information on included microarrays

According to the established inclusion criteria, the GSE137340, GSE69063, GSE69528, GSE54514, GSE57065, and GSE95233 expression datasets were used for analysis. Prior to the analysis, these datasets underwent batch effect removal followed by normalization and log2 transformation. Next, the gene expression profiles were visualized through boxplots, while significant differences between septic patients and healthy controls were shown by scatter plots in principal component analysis (Figs. S1A–1F). This study included 232 patients with sepsis and 116 healthy controls across six datasets. The basic workflow of the present study is shown in Fig. 1.

Functional enrichment analysis of the trained differential expression gene set based on RRA

A total of 361, 673, 674, 603 590, and 912 DEGs were identified in the GSE137340, GSE69063, GSE69528, GSE54514, GSE57065, and GSE95233 datasets, respectively (Fig. 2A). Subsequent analysis of six cohorts using RRA identified 95 upregulated and 69 downregulated DEGs. A smaller P-value in the RRA analysis was reflected by a higher gene rank of differential gene expression. The heatmap shown in Fig. 2B displays the top 10 upregulated and downregulated genes. Enrichment analyses, including GO and KEGG, were used to explore the biological function of the 164 DEGs. GO analysis revealed that the upregulated DEGs were mainly enriched in neutrophil migration and chemotaxis, positive regulation of NF-κB TF activity, and RAGE and Toll-like receptor binding. The downregulated DEGs were primarily enriched in Lymphocyte and T cell differentiation as well as T cell receptor complex and binding (Fig. 2C). KEGG pathway enrichment analysis showed that the predominant upregulated terms were IL-17, Toll-like receptor, and TNF signaling pathways, while the primary downregulated terms were T cell differentiation and T cell receptor signaling pathway (Fig. 2D). Collectively, these results suggest that the activation of pro-inflammatory signaling pathways may be concurrent with the inhibition of T cell differentiation and associated pathways in patients with sepsis.

Gene set enrichment analysis of statistically changed genes based on RRA

To overcome the limitations of differential enrichment analyses, such as the possibility of missing many moderate but meaningful changes due to an inappropriate threshold being chosen (Geng et al., 2022; Xia & Wishart, 2010), gene set enrichment analysis (GSEA) was conducted using c2.cp.reactome.v2022.1.Hs.symbols.gmt as the reference gene list. RRA analysis screened 2,763 genes with an adjusted P < 0.05. Ridgeline and GSEA plots were used to display the distribution of gene expression in the top six most significant gene sets. Notably, it was observed that the activated pathways were predominantly enriched in neutrophil degranulation (NES = 3.104, P = 1 × 10⁻¹⁰, q = 2.005 ×10⁻⁸) and innate immune system (NES = 2.604, P = 1 × 10⁻¹⁰, q = 2.005 ×10⁻⁸), while the inhibited pathway was enriched in co-stimulation mediated by the CD28 family (NES = −2.734, P = 2.002 × 10⁻⁶, q = 2.027 ×10⁻⁴; Figs. 3A and 3B). These findings indicated that innate immune responses may be activated and adaptive immunity may be inhibited in sepsis.

Immune landscape of peripheral blood in patients with sepsis

Among the above enrichment results, significantly changed genes were involved in several immune-related biological functions. Therefore, immune cell infiltration was evaluated across the aforementioned six datasets via CIBERSORT analysis. The distribution of 22 immune cell types is illustrated by stacked bar charts in Figs. 4A–4F. Cell types that were not present in all samples were filtered out, and their relative proportions were presented using boxplots (Figs. 4A–4F). Notably, the findings showed that the peripheral blood in patients with sepsis was infiltrated by higher proportions of monocytes, M0 macrophages, and neutrophils in the majority of datasets (≥4 datasets), but lower proportions of CD8 T cells, CD4 T cells, plasma, and resting natural killer cells. These results suggest that sepsis induces a reshaping of the composition and distribution of immune cells in the peripheral blood.

The advanced identification of the key gene set related to sepsis through WGCNA

WGCNA was used to analyze the GSE69528 microarray dataset, which comprises the largest number of samples within the six datasets, namely 83 patients with sepsis and 28 healthy controls. The present study constructed an adjacency matrix that followed a scale-free network by setting β = 7 (R² = 0.85) and maintaining high connectivity based on gene distribution (Fig. 5A). The hierarchical clustering method was used to construct gene clusters, and 11 co-expression modules were screened (Fig. 5B). Heatmap visualization revealed that the blue module exhibited the most significant positive correlation with sepsis, whereas the brown module demonstrated the strongest negative correlation with sepsis (Fig. 5C). Correlation analysis between gene significance (GS) and module membership (MM) revealed a strong association between the selected module genes and the blue MM (R = 0.82, P < 1 ×10⁻²⁰⁰; Fig. 5D), implying that genes within the blue module exhibited a strong association with sepsis. Based on the standard of MM > 0.8 and GS > 0.2, 424 genes were selected from the blue module, and their biological functions were revealed using GO and KEGG analyses (Figs. 5E and 5F). The results demonstrated significant enrichment in immune-related terms. The most significantly enriched GO terms for biological processes included leukocyte migration, activation of leukocytes, cells, and myeloid cells involved in immune response activation. For molecular functions, immune receptor activity and pattern recognition receptor activity were significantly enriched. Additionally, the primarily enriched KEGG pathways were neutrophil extracellular trap formation and NOD-like receptor signaling pathway. Ultimately, a gene set including 72 DEGs was obtained by overlapping 164 DEGs from the RRA analysis and 424 genes selected in the WGCNA (Fig. 5G).

The screening of hub genes in the key gene set via PPI analysis combined with Cytoscape software

Based on proteins encoded by the 72 DEGs in the gene set, a PPI network was constructed using the STRING online website (https://cn.string-db.org/). Subsequently, Cytoscape software was used to visualize the interaction network, which consisted of 43 nodes and 62 edges (Fig. 6A). The gene situated in the key node was regarded as a pivotal gene that potentially played a vital physiological function. Important gene modules were identified via the plug-in MCODE. By using default parameters, two gene clusters were screened—Cluster 1 and Cluster 2—and the genes in the two modules were labeled as Set 1 (Fig. 6B). Additionally, the plug-in cytoHubba-MCC was used to identify the top 10 genes (which were labeled as Set 2) and their scores (Figs. 6C and 6D).

It is widely acknowledged that gene expression is regulated in a spatiotemporal manner and forms intricate networks involving the interaction between transcription factors (TFs) and their direct target genes (Xu et al., 2022). Regulatory processes have important implications in various biological processes such as development, homeostasis, and pathogenesis (Xu et al., 2022). The present study employed the plug-in iRegulon of Cytoscape to investigate TFs that regulate Set 2. The findings revealed that TF ELF1 was responsible for regulating the eight genes within Set 2, excluding HP and MCEMP1 (Fig. 6E). By intersecting these eight genes with Set 1, five hub genes were screened: GPR84, hexokinase 3 (HK3), S100A12, CLEC4D and CLEC5A. The expression of these five hub genes was validated by conducting differential expression analysis in the validation dataset GSE28750 and by performing RT-qPCR analysis using blood samples from individuals with sepsis and healthy controls (Figs. 6F and 6G). The results were consistent with the prediction. Additionally, the expression of the predicted TF ELF1 was confirmed through RT-qPCR analysis in peripheral blood monocytes from patients with sepsis and mice after CLP (Table S2A and Fig. 2B).

The prognostic value evaluation of the hub genes by ROC analysis

If the DEGs identified in sepsis were of clinical significance, they should be indicative of disease severity and thus capable of predicting clinical outcomes. To evaluate the prognostic significance of the five hub genes, ROC analysis was performed in GSE54514, which includes survivor and non-survivor data. Unexpectedly, the prognosis of sepsis was significantly associated with four of the five hub genes: HK3, GPR84, CLEC4D, and S100A12. Among these genes, CLEC4D exhibited the highest prognostic value (AUC = 0.898) in distinguishing between sepsis survivors and non-survivors. The prognostic values of the other genes were as follows: GPR84 (AUC = 0.830), S100A12 (AUC = 0.830), and HK3 (AUC = 0.823), and CLEC5A (AUC = 0.741, a moderate prognostic value; Fig. 6H). These four genes were designated as the hub genes with prognostic significance.

The correlation analysis of the hub genes in patients with sepsis

Correlation analysis of the expression data from the GSE28750 validation cohort was conducted to investigate the expression patterns of the four significant hub genes. Prior to the analysis, batch effect removal, normalization, and log2 transformation were applied to process the dataset (Fig. S3A and Fig. 3B). A robust positive correlation was observed among the four hub genes (Fig. 7A). This result indicated that the four significant hub genes were potentially regulated by a shared TF. Additionally, the associations between the hub genes and immune cells were explored. The hub genes exhibited a strong positive association with innate immune cells including neutrophils, activated DCs, and M0 macrophages. Conversely, the hub genes demonstrated a negative association with acquired immune cells including B, CD4 T cells, and CD8 T cells (Fig. 7B). Moreover, most immune genes showed a negative correlation with the hub genes (Fig. 7C). The findings indicated that the four hub genes may be closely associated with immune dysfunction in septic patients.

Construction of a ceRNA network regulating the expression of hub genes

Extensive research on miRNAs is driven by their potential as therapeutic targets for a wide range of diseases. miRNAs can downregulate gene expression and induce gene silencing by binding mRNA molecules. However, lncRNAs, acting as upstream molecules, can influence miRNA function by interacting with miRNA response elements (MREs), leading to the upregulation of gene expression. This interaction between RNAs is called ceRNA (Li et al., 2018). Four online miRNA databases with default parameters were used to identify miRNAs that regulate ELF1 activity. Subsequently, miRNAs that were found to overlap in three databases were selected. In total, 98 miRNAs targeting ELF1 were screened, and a network was constructed using Cytoscape software (Fig. 8A).

The StarBase 2.0 database was used to investigate the interactions between the 98 miRNAs and lncRNAs, and a comprehensive literature search was also conducted. Based on the ceRNA hypothesis, one downregulated miRNA and one upregulated lncRNA previously reported in sepsis were selected (Meng et al., 2023), resulting in the proposal of the following ceRNA network: NEAT1-Homo sapiens (hsa)-miR-495-3p-ELF1. This network was visualized using Cytoscape software and is presented in Fig. 8B. It is speculated that this ceRNA network may serve as a critical regulatory pathway in sepsis.