This study leveraged biospecimens from the Sepsis Genomics Collaborative –a prospective multi-center observational cohort, which has been previously detailed extensively (12–15). Briefly, critically ill children between the ages of 1 week and 18 years who met consensus criteria for septic shock (16) were included and enrolled from 13 pediatric intensive care units (PICUs) across the United States from May 2015, through February 2019. The study protocol received approval from the Institutional Review Boards (IRBs) of the primary site (Cincinnati Children’s Hospital IRBs, Genomics of Septic Shock, IRB ID: 2008–0558 and 2022-0721) and all participating institutions. Informed consent was obtained from parents or legal guardians. All procedures involving human participants adhered to the ethical standards of the participating institutions’ IRBs, the 1964 Helsinki Declaration and its subsequent amendments. Whole blood collected in PAXgene RNA tubes and sera collected within 24 hours of the onset of septic shock were used. No study-related interventions occurred beyond these blood draws. De-identified clinical data were collected daily from days 1 to 7 of PICU admission, with mortality data tracked up to 28 days post-enrollment.

In addition to PICU admission platelet count, the 5 PERSEVERE biomarkers interleukin 8 (IL-8), heat shock protein (HSPA1B), granzyme B (GZMB), matrix metalloprotein 8 (MMP8), and C-C motif chemokine ligand 3 (CCL3) were previously measured in day 1 sera, permitting patient assignment to one of nine PERSEVERE-II Terminal Nodes (TNs) based on the published classification and regression tree (CART) (9, 10). Patients classified to TN 1, 2, 5, and 8 were predicted to be survivors and designated as low risk (<1.9% risk of death). Patients classified to TN 3, 4, 6, 7, and 9 were predicted non-survivors and considered high risk for mortality (16.7%-44.4% risk of death) (10). Demographic and clinical characteristics were compared between children stratified to low and high PERSEVERE-II risk. Pediatric Risk of Mortality (PRISM)-III scores were evaluated as an estimate of baseline illness severity (17). Immunocompromised status, vasoactive and corticosteroid use, prevalence of mechanical ventilation and renal replacement therapy were recorded. Outcome variables included PICU length of stay, PICU-free days, hospital length of stay, 7- and 28-day mortality, and the prevalence of complicated course. Complicated course was defined as the persistence of at least two organ failures at 7 days or mortality by 28 days. Dichotomous variables were compared with the Fisher exact test or chi-squared test. Nonparametric continuous variables were characterized as medians with interquartile ranges (IQRs) and evaluated with the Wilcoxon rank-sum test.

Whole blood was collected in PAXgene Blood RNA tubes and stored at −80 °C. For processing, tubes were thawed at room temperature for 2 h, inverted to homogenize, and 3 mL aliquots were transferred. RNA was isolated using a modified RNeasy Mini protocol on a QIAcube (QIAGEN). Briefly, PAXgene blood/stabilizer was diluted with PBS and centrifuged at 3,000 × g to pellet nucleic acids; pellets were washed with nuclease-free water, re-pelleted (3,000 × g), resuspended in Buffer RLT, treated with Proteinase K, and passed through gDNA-elimination columns. Flow-through was combined with isopropanol, bound to a MinElute column, washed with 80% ethanol, and eluted in RNase-free water. Eluates were heat-denatured (55°C, 5 min) and snap-cooled. RNA quantity was measured by Qubit RNA assays and integrity by BioAnalyzer; samples with RIN < 7 were excluded. Globin RNA was removed with GLOBINclear (Invitrogen) per manufacturer’s instructions. Globin-depleted RNA was quantified (Qubit RNA HS), and 10 ng was used for rRNA depletion and library construction with the SMARTer Stranded Total RNA-seq Kit v2—Pico Input Mammalian (Takara). Libraries were quantified (Qubit dsDNA HS), sized (Fragment Analyzer High Sensitivity Small Fragment kit), pooled, and sequenced on an Illumina NovaSeq 6000 (paired-end, 2 × 100 bp). Per sample, 40–120 million read pairs were generated. FASTQ files were used for downstream processing. Library prep and sequencing were performed at TB-SEQ (Palo Alto, CA).

Raw mRNA counts were mapped to 60,846 Ensembl Gene IDs across samples. The Gene IDs which did not correspond to a known Human Genome Gene Symbol and Entrez ID were removed, leaving 20,239 genes. In addition, 171 Ensembl Gene IDs were true duplicates with identical mRNA count data and were removed. There were nine duplicate Human Genome Gene Symbol pairs, each of which had one member that was less expressed than the other and removed. Finally, there were 26 unique Entrez IDs corresponding to two or three Gene Symbols, 55 Gene Symbols in all. These were manually evaluated, and the Gene Symbol with the lowest average expression was removed. Altogether, this left a gene expression matrix with 20,030 genes across 81 samples.

We conducted Weighted Correlation Network Analysis (WGCNA) (18) to identify gene-modules, representing co-expressed genes, associated with biomarker mortality risk strata and to explore relationships among genes. First, the raw gene expression matrix was normalized using trimmed-mean of M-values (TMM) (19) and low-expression values were removed using edgeR (20) using the default min.count = 10, leaving 13,515 genes. This normalized, log2-transformed matrix was used for network construction. A WGCNA soft threshold of 12 was selected with R² = 0.817 and mean connectivity = 70.200, meeting our goals of scale free topology model fit coefficient > 0.80 and mean connectivity < 100 to achieve maximum correlation strength in addition to appropriate hub connectivity for analysis. (Supplementary eFigure S1). Pearson’s correlation coefficients were calculated to assess the strength of correlation. Signed network construction was utilized to emphasize directional correlation relationships between genes. For each gene co-expression module of interest, its gene expression between conditions was evaluated by comparing module eigengene expression for high- and low- mortality risk patients. WGCNA calculates an eigengene expression value for every patient across every module which represents the first principal component of the gene expression for that patient across all genes in that module. The eigengene expression of high-risk and low-risk patients were compared using a two-tailed heteroscedastic t-test. WGCNA also computes each gene’s correlation with its module and the significance of that correlation; the genes in each module whose correlation with that module has the most significant p-value were identified as driver genes. Similarly, WGCNA calculates each gene’s correlation with the trait of interest and the significance of that correlation, which revealed the genes whose expression was most associated with the high-risk strata. The gene lists in modules of interest were submitted to the Database for Annotation, Visualization and Integrated Discovery (DAVID) to identify functional annotations (21); with only those with a Benjamini-Hochberg adjusted p-value < 0.05 being considered significant.

We identified differentially expressed genes (DEGs) comparing high- and low-mortality risk patients using R package DESeq2 (22). We used a Benjamini-Hochberg adjusted false discovery rate (FDR) threshold of 0.05 to identify DEGs. Heatmap and volcano plots were used to visualize DEGs. The biological relevance of pathways were determined based on Gene Ontology (GO) annotations and Kyoto Encyclopedia of Genes and Genomes (KEGG) using clusterProfiler (22).

We inferred relative differences in cell type abundance comparing risk-strata using CIBERSORTx (23) based on differentially expressed genes. However, this computational tool was originally designed for in silico tissue deconvolution rather than blood and lacks a reference for cell types specific to critically ill patients. To address this limitation, we used single-cell RNA sequencing dataset comprised of critically ill adults with sepsis published by Kwok et al. (11) We calculated a composite gene score as the geometric mean of top 20 overexpressed genes minus the geometric mean of top 20 repressed genes using published methods (24), identified through DEG analyses comparing risk-strata and also available in the single-cell dataset. We mapped this scaled composite score against the Uniform Manifold Approximation and Projection (UMAP) of the Kwok dataset to infer cell types contributing to biological differences between risk-strata.

We used the Kwok dataset to generate pseudobulk gene-expression data. For each individual patient represented in this dataset, gene-expression in each cell type was aggregated using the AggregateExpression function in Seurat. The pseudobulk matrix was batch-corrected using the ComBatseq function from the sva package (v.3.46.0), a negative binomial regression method. DEGs comparing risk strata, identified previously, were selected as features for downstream analyses. To assign risk strata in the reference single cell dataset, we developed a Support Vector Machine (SVM) classification model. The normalized matrix and corresponding labels were randomly split into training and validation sets (80:20 ratio) to train and fine tune the SVM model. This model was applied to the corrected pseudobulk data to assign high- or low-risk labels to all single-cell sepsis samples. Finally, CellChat analysis (v2.1.0) was performed on the single-cell matrix, to infer cell-cell communication networks focusing on the high-risk strata.

To identify key transcriptional regulators among high-risk patients, we submitted gene lists in each WGCNA module with statistically significant association with high-risk strata to the Chip Enrichment Analysis (ChEA3) portal (https://maayanlab.cloud/chea3/) to predict transcription factors (TFs) anticipated to regulate gene co-expression module (25). The most notable transcription factor (TF) is denoted by the lowest mean rank, which indicates the TF predicted by ChEA3 to interact most with the submitted gene lists after searching across multiple libraries including ENCODE, GTEx, ARCHS4, and ReMap. Additionally, we submitted DEGs distinguishing patient risk-strata for Ingenuity Pathway Analysis (IPA, QIAGEN) (26) (QIAGEN) to identify upstream regulators and mechanistic networks that could influence gene expression patterns, focusing on direct interactions between regulators and selecting TFs with highest activation z-scores and p<0.001.

We identified 11 gene co-expression modules using Weighted Correlation Network Analysis (WGCNA). The correlation between these modules and the high-risk strata was calculated, as was the significance level corresponding to each correlation as shown in Figure 1. Fifty-one genes (designated the Pink module; Supplementary Table S1a) had the highest correlation (R² = 0.436) and the most significant p-value for the correlation with the high-risk mortality strata (p = 4.795e^-5). Notably, the gene most strongly associated with this module was MMP8 (p = 1.53e^-36), which encodes for one of the PERSEVERE-II biomarkers. The 5th- and 6th-most strongly associated genes with this module, LCN2 and RETN, respectively, were two of the original 12 candidate PERSEVERE biomarker genes that were pruned from the final model. The gene OLFM4, previously associated with neutrophil subpopulations in pediatric septic shock (27), also has a strong association (p = 7.91e^-15). Three Reactome pathways were significantly associated with this module: “Neutrophil Degranulation,” “Innate Immune System,” and “Immune System.” Three other gene-modules associated with mortality-risk strata are detailed in Supplementary Tables S1b-d. Comparison of gene module eigengene expression between mortality-risk strata is shown in Supplementary eFigure S2, with this 51-gene (“Pink”) module once again showing the greatest difference between groups.

We identified 260 genes that were significantly associated with the high-risk strata at p < 0.001 independent of their gene co-expression module membership (Supplementary Table S2). There was a high degree of overlap between these genes and DEGs identified by DESeq2, detailed subsequently. For example, the top seven genes identified by WGCNA associated with the high-risk cohort were identified by DESeq2 as genes with significantly lower expression in high-risk patients (NHSL2: p < 1.9e-11, LRMP: p < 4.14e-10, IL16: p < 6.16e-9, GLIPR1: p < 7.65e-9, SRPK2: p< 7.73e-9, AOAH: p< 1.33e-9, CALB1: p< 9.55e-9, Supplementary Table S4a). In addition, examining co-expression modules associated with high-risk strata, there were 9 Pink Module genes, 6 Yellow Module genes, and 24 Brown Module genes identified by WGCNA as significantly associated with the high-risk strata independent of module membership. All of these 39 genes were also shown to be over-expressed among high-risk patients with a p < 0.001 by DESeq2 (Supplementary Table S3).

Heatmap visualizing differentially expressed genes comparing high vs low mortality-risk patients is shown in Figure 2a. A total of 2,654 genes (13.3% of all sequenced genes) were differentially expressed at an adjusted p-value < 0.05, of which 1,602 genes were over-expressed and 1,052 were under expressed among high-risk patients relative to those at low-risk of mortality (Supplementary Table S4a). Of note, genes coding for 4 out of the 5 PERSEVERE-II biomarkers - GZMB, CXCL8 (IL-8), HSPA1B, and MMP8 were overexpressed DEGs among high-risk relative to low-risk patients. Figure 2b shows the volcano plot highlighting the most differentially expressed genes (DEGs) based on a log2FC threshold of > ± 1 and adjusted p value of <0.001.

As shown in Figure 3, top panel), genes over-expressed in high-risk samples were involved in the cell cycle, with the most-enriched GO terms being “mitotic cell cycle phase transition,” “chromosome segregation,” “organelle fission,” “nuclear division,” and “regulation of cell cycle phase transition.” The genes repressed among high-risk patients were involved in “positive regulation of cytokine production,” “activation of immune response,” “immune-response regulating signaling pathway,” “immune-response activating signaling pathway,” and “leukocyte-mediated immunity.” Alternatively, GO annotations “T cell proliferation” and “regulation of T cell activation” (Figure 3, bottom panel) and KEGG pathways “T cell receptor signaling pathway” and “B cell receptor signaling pathway” (Supplementary eFigure S3), reflective of the adaptive immune response, were repressed among high-risk patients.

We identified that high mortality-risk patients had a lower average fraction of mature neutrophils (35.5%; 95% CI [29.0%, 42.0%]) compared with low mortality-risk patients (46.3%; 95% CI [42.6%, 50.0%]) based on CIBERSORTx analyses shown in Figure 4. There was no appreciable difference in the proportion of any of the other 21 cell types imputed by CIBERSORTx. As shown in Figure 5, the Kwok et al. dataset had 10 cell types from critically ill adult patients with sepsis. Genes upregulated among high-risk patients were expressed primarily by a small population of developing neutrophils. Further, downregulated genes among patients with high mortality risk were expressed primarily by mature neutrophils. We further identified that developing neutrophils exhibited greater number and strength of inter-cellular interactions, specifically among high-mortality risk patients.

Gene lists of WGCNA submitted to ChEA3, identified LTF, itself a gene in the Pink Module, as the most enriched TF with a mean rank of 1.0 and 10 overlapping genes including MMP8, LCN2, and RETN. Other TFs of interest identified included KLF1 (Brown Module, mean rank of 1.6 and 231 overlapping genes), FOXM1 (Yellow Module, mean rank of 2.6 and 192 overlapping genes), and ZNF12 (Purple Module, mean rank of 8.0 and 13 overlapping genes). Among the top regulator effect networks associated with high mortality-risk patients, identified through IPA analyses of DEGs, was CEBPB with an activation z-score of 5.08 and p-value of 1.11 e^-12, indicating the degree of overlap in genes from the dataset and those modulated by the particular TF. Other key TFs identified included TFEB, MYC, and TBX3. Consistent with results using the former approach, FOXM1 was identified to be activated with a z-score of 3.3 and p-value of overlap of 5.2 e^-11 and KLF1 had an activation z-score of 1.7 and p-value of 1.3 e^-16. In contrast, LTF and ZNF12 were not enriched when using causal network analyses in IPA (Supplementary Table S5).