The HTA2.0 microarray data of 156 patients with sepsis and 82 healthy subjects was downloaded from GEO with primary data accession number GSE134347 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE134347). Patient characteristics were described in (Scicluna et al., 2020).

Screening of correlated gene pairs was performed by an R package WGCNA (Langfelder and Horvath, 2008). And step-by-step calculation was as follows:

Step 1 The Pearson correlation coefficient between gene X _i and X _j was calculated, and matrix X was converted into the correlation matrix S between genes: S i j = | c o r ( X i , X j ) | (1)

Step 2 Weighted adjacency matrix aij was constructed and suitable soft-threshold power β was selected: a i j = | S i j | β (2)

Step 3 The degree of separation of nodes was calculated. The adjacency matrix is converted into unsigned topological overlap matrix (TOM) to calculate degree of intergenic dissimilarity (Zhang and Horvath, 2005), based on which genes were distributed in different modules: T O M i j = l i j + a i j min ( k i , k j ) + 1 − a i j ; l i j = ∑ u a i u a j u (3)

In the above formula, l _ij represented the sum of product of adjacency coefficients of all common adjacent genes of gene i and j. a _ij represented the adjacency coefficient between gene i and j. k _i represented synthesis of adjacency coefficients of gene i with all neighboring nodes. If TOM _ij was 0, it meant that gene i and j were isolated and not connected to all other genes. If TOM _ij was 1, it meant that these two genes were adjacent to all surrounding genes and were also connected to each other. In other words, TOM represented the similarity of genes, so the dissimilarity between genes could be calculated: D i s t T O M i j = 1 − T O M i j (4)

Step 4 DistTOM _ij was used for hierarchical clustering, and genes were divided into different co-expression modules (Li and Horvath, 2007). WGCNA adopted Dynamic Tree Cut Method to construct the cluster tree, which was a top-down merging algorithm. Through iteration and decomposition of gene clusters, stable gene clusters were eventually achieved (Langfelder et al., 2008). Here, the minimum module size was set as 30 to identify modules and draw dendrogram.

Step 5 Identify trait-related modules. Herein, trait was defined as disease. We defined gene significance (GS) as relationship between gene expression levels and disease. Moreover, module membership (MM) represented the degree of relationship between module feature genes and disease. The higher MM represented the higher correlation between modules and disease. At last, we identified two modules most relevant to disease. One module exhibited positive correlation and another negative correlation.

R package DiffCorr was utilized for the visualization and identification of differential correlations in biological networks. This package was based on Fisher’s z-test and details were explained in (Fukushima, 2013; Zhou et al., 2021).

R package clusterProfiler was implemented to conduct enrichment analysis of clustered genes in blue and darkred modules. We used a hypergeometric distribution test for the classification of enrichment terms. And p values were adjusted by false discovery rate (FDR) method, the cutoff of which was set to be 0.05 (Yu et al., 2012).

Cytoscape (3.9.0) was used to realize visualization of networks (Doncheva et al., 2019).

We applied Student’s t-test to identify genes differentially expressed between sepsis and healthy samples. p values were adjusted by the Benjamini–Hochberg method (Hochberg and Benjamini, 1990). Differentially expressed genes were defined as adjusted p value less than 0.05. We employed Fisher’s z-test to evaluate differential correlations of gene pairs between sepsis and healthy patients. And lfdr less than 0.05 was regarded as significant differential correlations.

Pearson’s correlation coefficient was applied to cluster samples from GSE134347. After removing outliers, a sample clustering tree was drawn (Figure 1B). Co-expression network was constructed from 25,245 coding and non-coding genes through WGCNA approach. We set soft-thresholding power five for satisfying scale-free topology of network, in which the corresponding R² was 0.81 (Figure 1C). And we detected twenty-six modules in the network, as shown in a cluster dendrogram (Supplementary Material S1). The members in each module were listed in Supplementary Material S2. Apart from the grey module consisted of many un-classified members, orange module contained the minimum 33 genes, while the maximum 10,821 genes were included in turquoise module.

Next, we quantified module-trait associations (Figure 2), in which the blue and darkred modules exhibited the most significant associations with sepsis. The corresponding correlation coefficients of blue and darkred modules were 0.88 (p = 2 × 10–78) and −0.77 (p = 2 × 10–48), respectively. In addition, GS and MM analysis demonstrated that genes highly significantly associated with sepsis were also the most crucial factors of modules associated with sepsis (Supplementary Material S3).

Then, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of these two modules were performed. As shown in Figure 3A, genes in the blue module were significantly enriched in ncRNA metabolic process, Herpes simplex virus 1 infection, and Th17 cell differentiation. Notably, as a less well studied subset of CD4⁺ Th cells, reduced Th17 cell responses have been observed in sepsis on account of low expression levels of involved transcription factors, bringing about increased susceptibility of patients to secondary fungal infections (Monneret et al., 2011; Hotchkiss et al., 2013). Our study provided additional evidence on the pivotal roles of Th17 cell function in sepsis. As shown in Figure 3B, genes in the darkred module were enriched in neutrophil activation involved in immune response and neutrophil degranulation. Consistent with previous studies, abnormal behaviors of neutrophil including delayed apoptosis have been witnessed during the early stage of sepsis (Drifte et al., 2013). Meanwhile, the most severely reduced neutrophil function has been investigated in patients with the highest risk of acquiring nosocomial infections (Nedeva, 2021), bestowing the importance of neutrophil responses on sepsis occurrence and progression.

We further chose genes in the blue and darkred modules to estimate differential correlations, which were grouped based on expression patterns in each subtype (sepsis or healthy) under the cluster. molecule function of DiffCorr package. We utilized (1-correlation coefficient) as a distance measure based on cutree function. Two functions named get. eigen.molecule and get. eigen.molecule.graph were applied for module networks visualization (Figure 4). The comp.2. cc.fdr function offered the resultant pair-wise differential correlations among blue and darkred modules.

R package DiffCorr also identified oppositely correlated pairs. For example, two genes positively correlated in sepsis tissues and negatively correlated in healthy tissues, or vice versa. These switched gene pairs were worthy noticed for their crucial roles in understanding molecular mechanisms in the progression of sepsis. Totally one hundred and seventy-four oppositely correlated gene pairs from the blue module and forty-nine from the darkred module were obtained (Supplementary Material S4), whose interaction networks were presented in Figure 5. The top ten significant differentially correlated gene pairs between sepsis and healthy tissues from blue and darkred modules were shown in Table 1.

Lastly, we focused on four hub genes acting as master regulators due to possessing most links with other genes. Zinc finger protein 366 (ZNF366), also known as DC-SCRIPT, belongs to the zinc ring finger protein family and has recently been reported to regulate dendritic cell development (Søndergaard et al., 2015; Wang et al., 2018). ZNF366 has been shown to suppress toll like receptor-mediated expression of IL-10 through modulating NF-κBp65 activation (Søndergaard et al., 2015). Decreased numbers and disabilities of dendritic cells have been widely observed in sepsis followed by immune responses alterations (Roquilly and Villadangos, 2015; Venet and Monneret, 2017; Wang et al., 2020). Additionally, sepsis-induced dendritic cell blockade has been reported to prevent mice from sepsis-induced death (Meng et al., 2017). We also observed its elevated levels in sepsis compared to healthy samples (Figure 6A). Given the great importance of ZNF366 in coordinating dendritic cell function, our study expanded its roles in sepsis progression which needed further validations.

Another zinc finger protein that takes the central position in the differentially co-expressed gene network is MYND-type containing 11 (ZMYND11), which is the negative regulator of NF-κB signaling and vastly impacts the replication of various viruses including Hendra virus and Epstein-Barr virus (Ikeda et al., 2010; Stewart et al., 2013). As previously described, numerous patients with sepsis witnessed the activation of NF-κB pathway initiated by pathogen-associated molecular pattern or danger-associated molecular pattern, emphasizing its great contributions during sepsis (Hayden et al., 2006; Chen et al., 2019; Wang et al., 2019). In line with it, significant lower levels of ZMYND11were detected in sepsis than healthy samples (Figure 6B). Also, impaired NF-kB activation has been correlated with sepsis-induced acute lung injury (Chen et al., 2021). Thus, it could be inferred that ZMYND11 participated in sepsis progression through NF-κB signaling and further mechanisms needed to be explored.

Small p97/VCP-interacting protein (SVIP), localized to the ER membrane by myristoylation, is highly expressed in central nervous system and related to autophagy modulation (Wu et al., 2013; Jia et al., 2019). Further mechanisms have demonstrated that overexpression of SVIP protected hepatocytes from the toxicity of CCL₄ through enhancing LC3 lipidation and activating autophagy (Jia et al., 2019). Prior studies have also found its essential roles in lysosomal dynamic stability and autophagosomal-lysosomal fusion (Johnson et al., 2021). Autophagy dysregulation has been observed in organ injury induced by sepsis, in which autophagy exerted vital effects on programmed cell death pathway and inflammation (Lo et al., 2013; Qiu et al., 2018; Zhao et al., 2019). Meanwhile, animal experiments have endowed autophagy with protective roles in septic brain injury (Su et al., 2015). Also, higher levels of SVIP were observed in sepsis compared to healthy samples (Figure 6C). Based on the above evidence, our study firstly illustrated the underlying roles of SVIP in sepsis through autophagy-related pathways.

Firstly identified in yeast and human placenta, ubiquitin conjugating enzyme E2 H (UBE2H) belongs to the structurally and functionally conserved family of E2s and is involved in ubiquitination and proteasome-mediated protein degradation and regulated by TNF-α signaling (Kaiser et al., 1994; Li et al., 2003; Clague et al., 2015). Increasing evidence has indicated the potential roles of UBE2H in human brain diseases such as amyotrophic lateral sclerosis, Alzheimer’s disease and autistic disorder (Martin et al., 2009; Sokolowski et al., 2018; Lim and Joo, 2020). Nevertheless, there is scant report on the links between UBE2H and sepsis. Previous data has shown that body protein loss during sepsis was caused by upregulation of ubiquitin genes and ubiquitin-proteasome pathway (Hasselgren and Fischer, 1997). Moreover, as a hypoxia-mediated gene, UBE2H may be speculated to participate in tissue hypoxia in septic shock (Rello et al., 2017). Here, we also found elevated levels of UBE2H in sepsis compared to healthy samples (Figure 6D). Therefore, our study firstly revealed the function of UBE2H in sepsis, offering experimental clues for further investigations.