Whole blood transcriptomics identifies subclasses of pediatric septic shock

Abstract

Background

Sepsis is a highly heterogeneous syndrome, which has hindered the development of effective therapies. This has prompted investigators to develop a precision medicine approach aimed at identifying biologically homogenous subgroups of patients with septic shock and critical illnesses. Transcriptomic analysis can identify subclasses derived from differences in underlying pathophysiological processes that may provide the basis for new targeted therapies. The goal of this study was to elucidate pathophysiological pathways and identify pediatric septic shock subclasses based on whole blood RNA expression profiles.

Methods

The subjects were critically ill children with cardiopulmonary failure who were a part of a prospective randomized insulin titration trial to treat hyperglycemia. Genome-wide expression profiling was conducted using RNA sequencing from whole blood samples obtained from 46 children with septic shock and 52 mechanically ventilated noninfected controls without shock. Patients with septic shock were allocated to subclasses based on hierarchical clustering of gene expression profiles, and we then compared clinical characteristics, plasma inflammatory markers, cell compositions using GEDIT, and immune repertoires using Imrep between the two subclasses.

Results

Patients with septic shock depicted alterations in innate and adaptive immune pathways. Among patients with septic shock, we identified two subtypes based on gene expression patterns. Compared with Subclass 2, Subclass 1 was characterized by upregulation of innate immunity pathways and downregulation of adaptive immunity pathways. Subclass 1 had significantly worse clinical outcomes despite the two classes having similar illness severity on initial clinical presentation. Subclass 1 had elevated levels of plasma inflammatory cytokines and endothelial injury biomarkers and demonstrated decreased percentages of CD4 T cells and B cells and less diverse T cell receptor repertoires.

Conclusions

Two subclasses of pediatric septic shock patients were discovered through genome-wide expression profiling based on whole blood RNA sequencing with major biological and clinical differences.

Trial Registration This is a secondary analysis of data generated as part of the observational CAF-PINT ancillary of the HALF-PINT study (NCT01565941). Registered March 29, 2012.

Background

Septic shock, a severe form of sepsis, characterized by profound circulatory and metabolic abnormalities that require the need for a vasopressor to maintain mean arterial pressures, is associated with mortality rates greater than 40% [1]. Despite efforts to develop therapies beyond current standard practices over the past thirty years, over 60 large sepsis clinical trials have failed to identify significant positive results [2]. Many experts now believe that the failure is largely due to the heterogeneity of sepsis patients [3]. This has prompted investigators to develop a precision medicine approach aimed at identifying biologically homogenous subgroups of patients with septic shock and critical illnesses [4,5,6]. The hope is to develop targeted therapeutic interventions for sepsis patients based on pathophysiological processes and clinical phenotypes rather than a one-size-fits-all-approach. [7,8,9,10,11].

Wong et al., pioneers in the field, were the first to delineate septic shock subgroups among pediatric patients, based on underlying biological perturbations and establish a potentially prognostic and predictive sepsis subclassification system. Their work was based on microarray-based gene expression data, discovering distinct subclasses characterized by differential gene expression related to adaptive immunity and glucocorticoid receptor signaling. Remarkably, these subclasses, derived independently of clinical data, exhibited variations in mortality rates with a suggestion of differential steroid benefit between subclasses, highlighting their clinical significance. However, differences, if any, in circulating protein biomarkers and immune cell repertoires between these subclasses are not known [5].

Although microarray platforms have traditionally served as established and reliable tools for gene expression profiling, the emergence of RNA sequencing (RNA-seq) introduces a potent technology capable of directly sequencing amplified cDNA to measure transcript abundance. This relatively recent approach enables sensitive transcript detection and encompasses a broader quantitative spectrum of expression level changes compared to microarrays. Consequently, it yields more accurate estimations of absolute transcript levels, identifies a greater number of differentially expressed protein-coding genes, and exhibits enhanced concordance between RNA-Seq and protein expression measurements. Moreover, RNA-Seq holds the potential for comprehensive biological exploration such as cell deconvolution and imputation of B and T Cell repertoires surpassing the capabilities of microarrays [12,13,14].

In this study, we leveraged available RNA-Seq data from a clinical trial of children with circulatory and/or respiratory failure and compared whole genome expression between cases with septic shock and mechanically ventilated controls. Next, we used an unsupervised discovery-based approach to identify two pediatric septic shock subclasses based on their transcriptomic signatures and compared their clinical characteristics, ionotropic use, and circulating biomarkers. Finally, to further expand on the knowledge of biological differences between the subclasses, we imputed immune cell abundance, and T and B cell repertoires from the sequencing data and compared these between the subclasses.

Methods

Aim

The first aim was to elucidate genetic and biological differences in pediatric septic shock patients compared to critically ill controls. The second aim was to identify clinically significant pediatric septic shock subclasses based on whole blood RNA expression profiles and compare clinical characteristics, circulating protein biomarkers, immune cells, and adaptive immune repertoires.

Design, setting, patient characteristics

This was a nested case–control study. Subjects included in this study were a subset of the patients from the coagulation and fibrinolysis in pediatric insulin titration trial (CAF-PINT) [11] an ancillary to the heart and lung failure-pediatric insulin titration trial (HALF-PINT) (ref 23). All were mechanically ventilated, had blood glucose levels greater than 150 mg per deciliter and were treated with continuous insulin infusion randomized to one of two targeted glycemic control ranges. Blood samples were collected prior to initiation of the randomized intervention. Septic shock cases were defined as patients who had a documented clinical diagnosis of sepsis and were started on inotropes within 72 h prior to blood sampling. Controls were defined as patients who had no documentation of infection or sepsis, no positive cultures and did not require inotropic support prior to randomization (Fig. 1, Table 1) (Additional file 1, 2).

Fig. 1

Flowchart of patients in this study after applying exclusion and inclusion criteria. All septic shock cases were diagnosed with sepsis and were started on inotropes < 72 h from blood sampling. All controls had no documented diagnosis of sepsis or pneumonia, no positive cultures, and no inotrope use, and were started on mechanical ventilation < 72 h from blood sampling

Table 1 Baseline Characteristics of Septic Shock Cases and Controls

RNA isolation, sequencing, and expression quantification

Total RNA was extracted from whole blood using the PAXgene Blood RNA kit modified for pediatric use. Next, sequencing libraries were prepared using the Nugen universal plus kit with polyA capture and sequenced with the NovaSeq S4 system (Illumina) to generate 2 × 150 base paired-end reads to a target depth of 50 million read-pairs per sample. After quality control, 20,010 protein-coding genes were left for analysis. Additional details on these methods are found in the Additional file 1: Online Data Supplement.

Differential gene expression analysis

For Aim 1, we compared baseline gene expression between cases and controls using negative binomial generalized linear models (EdgeR). Differentially expressed genes (DEGs) were defined as genes with an absolute log₂-fold change of greater than 0.5 and a false discovery rate (FDR) of less than 0.05. Then, we used EnrichR, a pathway enrichment tool, to identify significant GO terms and KEGG pathways in our list of DEGs [15]. Terms and pathways with adjusted p values less than 0.05 were retained. For Aim 2, patient values for these statistically significant DEGs underwent variance-stabilizing transformation using the DESeq 2 package and were plotted in a heatmap with unsupervised hierarchical clustering. The DESeq2 package implements variance-stabilizing transformation which is a logarithmic transformation that aims to stabilize the variance across the range of expression levels. Spearman correlation was used for both genes and samples, and complete linkage method was used. Unsupervised hierarchical clustering is a technique to group individuals into clusters based on similarities, with the distance or difference between individuals represented by a dendrogram. As was previously done by Wong et al. [16], two septic shock subclasses were ultimately chosen based on first-order dendrogram branching (Fig. 3A), as Subclass 2 patients appeared visually more similar to controls, and Subclass 1 was clustered separately in a different clade. A K means clustering analysis was also performed, and using the silhouette method, we determined two clusters was most optimal which confirmed two septic shock subgroups as most optimal. Overall the K means cluster analysis supported the original analysis with a notable 78% agreement in subclass assignment providing substantial robustness with our original clusters. Clinical features, circulating biomarkers, immune cells and adaptive immune repertoires were then contrasted using subclass membership. We also compared gene expression and performed pathway enrichment to analyze the biological differences between these two subclasses.

Comparing clinical characteristics

We compared the clinical outcomes between cases and controls and then between the two subclasses using nonparametric linear regression analysis with bootstrapping, adjusting for baseline differences between treatment groups including PRISM score and steroid use.

Biomarker quantification

Plasma levels of twelve biomarkers of inflammation and endothelial injury were assayed using the Human Magnetic Luminex Screening Assay with the Luminex 200 (R&D Systems) according to the manufacturer’s instructions at the UCLA Immune Assessment Core facility. Values were compared across septic shock subclasses using the Wilcoxon rank-sum test and the false discovery rate criterion to adjust for multiple testing.

Measuring whole blood cell type composition

We applied the gene expression deconvolution interactive tool (GEDIT), a deconvolution tool that utilizes bulk gene expression data to model the most likely combination of cells that would produce the presented bulk transcriptome to estimate cell type fractional composition within whole blood [17, 18]. The abundance of each cell type in each subclass was then compared using a Wilcoxon rank-sum test, followed by a Bonferroni correction.

Profiling of adaptive immune repertoires

We then sought to quantify adaptive immune responses based on the recombination landscape of genes encoding B and T cell receptors (BCR and TCR). We utilized ImReP, a computational method for profiling the adaptive immune repertoire from RNA-Seq data [19]. Subsequently, a comparison between the two subgroups was performed using a Wilcoxon rank-sum test on the number of reads associated with each receptor type.

Results

Patients

Of 303 CAF-PINT patients, we identified n = 46 with septic shock and n = 52 controls without sepsis or shock (Fig. 1). Septic shock patients were older, whereas controls were younger and more likely to be admitted for asthma, bronchiolitis, postoperative care, or trauma (Table 1).

Identification of septic shock subclasses

We first identified genes that were differentially expressed between septic shock cases and controls. A total of 840 DEGs were identified: 530 were upregulated (Additional File 2) and 310 were downregulated (Additional File 3) in patients with septic shock compared to controls (FDR < 0.1). Pathway enrichment analysis using EnrichR showed that septic shock patients had increased expression of pathways related to neutrophil degranulation and glutathione metabolism and decreased expression of pathways related to antigen processing and presentation and T cell activity (Fig. 2, Additional file 4: Table E1). The 840 genes were subjected to unsupervised hierarchical clustering based on a correlation distance matrix as shown in the heatmap in Fig. 3A. Overall, septic shock cases appear clustered to the right, and critically ill controls are clustered to the left. Interestingly, there was heterogeneity among the septic shock cases. The transcriptomic signature of cases on the right appears significantly different from the septic shock cases on the left, which are more similar to controls. Based on first-order branching of the dendrogram of the columns, two major septic shock subclasses were identified. We designated these subclasses arbitrarily as Subclasses 1 and 2 (Fig. 3B).

Fig. 2

Gene network displaying functional enriched terms with associated DEGs between septic shock cases and controls (A). Upregulated pathways in septic shock cases compared to controls (B). Downregulated pathways in septic shock cases compared to controls

Fig. 3

Heatmaps A Heatmap depicting the expression of all 840 differentially expressed genes (DEGs) between septic shock cases and controls. Each row represents a DEG; blue indicates that a gene is downregulated, while red indicates that a gene is upregulated. Each column is a patient, so each column depicts the transcriptomic signature of each patient. The dendrogram at the top shows unsupervised hierarchical clustering of the patients such that the patients with similar transcriptomic signatures are clustered together. B This heatmap depicts the same patients and genes as in Fig. 2A, but now the controls are sorted on the left, and the cases sorted on the right. Clustering based on each patient’s transcriptomic signature was performed within each group

Differing clinical characteristics between the septic shock subclasses

Given the biological differences between the two subclasses, we investigated differences in clinical outcomes. We first compared baseline characteristics between the two groups, and the demographic and clinical data are shown in Table 1. We found that Subclass 1 patients were on average significantly older and had higher average blood glucose levels at baseline. There was no difference between subgroups in ICU admission illness severity, as calculated by the pediatric risk of mortality (PRISM) score. Despite the absence of a difference in mortality risk, these two subgroups had significantly different clinical outcomes (Table 2, Additional file 4: Table E9). Over the course of the study, Subclass 1 patients had statistically higher maximum pediatric logistic organ dysfunction (PELOD) scores indicating more severe multiple organ dysfunction. Subclass 1 also had higher maximum Vasoactive-Inotropes Scores, indicating that they required more cardiovascular support. Subclass 1 also had fewer hospital-free days and longer ICU stays than Subclass 2. There were no statistically significant differences in ventilator-free days. No differences in mortality were seen, as only three patients total died after 90 days of follow-up.

Table 2 Clinical characteristics and outcomes of Subclass 1 and Subclass 2

Biomarker analysis

Given the finding of greater innate and myeloid-derived inflammation as well greater circulating endothelial cells among patients in Subclass 1, we next tested whether gene expression findings could be recapitulated using plasma protein markers. We examined a panel of twelve candidate inflammatory biomarkers of sepsis to identify differences between the two subclasses. These markers were initially chosen to provide insights into derangements of inflammation and thrombosis pathways based on their association with hyperglycemia and pediatric critical illness [20,21,22]. Subclass 1 had significantly higher levels of PAI1, IL6, IL8, IL10, ANG2, TREM1, and TNFR than Subclass 2 (Table 3). There were no significant differences between subclasses of IL4, thrombomodulin, TFPI, P-selectin, and ICAM1. In comparing each subclass individually to controls, Subclass 1 had significantly higher levels of PAI1, IL6, IL8, IL10, ANG2, TREM1, and TNFR-1, while Subclass 2 and controls did not have many significant differences relative to non-septic critically ill controls (Additional file 4: Table E6).

Table 3 Differences in Inflammatory Biomarker Levels between Subclasses 1 and 2

Differential gene expression across septic shock subclasses

Using these two septic shock subclasses, Subclass 1 (n = 21) and Subclass 2 (n = 25), as comparison groups, we performed differential gene expression analysis and found a total of 2486 DEGs. A total of 1275 of these genes were upregulated (Additional File 5) and 1211 were downregulated (Additional File 6) in Subclass 1 compared to Subclass 2 (FDR < 0.05). Pathway enrichment analysis of the 2486 genes revealed upregulation in neutrophil-mediated immunity and downregulation in adaptive immunity pathways including B and T cell activity, in Subclass 1 compared to Subclass 2 (Fig. 4, Additional file 4: Table E2).

Fig. 4

Gene network displaying functionally enriched terms with associated DEGs between Subclass 1 and Subclass 2. A Upregulated Pathways in Subclass 1 compared to Subclass 2 B Downregulated Pathways in Subclass 1 compared to Subclass 2

Cell deconvolution analyses and profiling of adaptive immune repertoires

To quantify cell type populations in our samples, we performed expression-based cell deconvolution analysis [17]. In concordance with our pathway analysis, we found that compared to Subclass 2, Subclass 1 had significantly lower percentages of CD4 + T cells (0.9 ± 1.5 vs. 5.2 ± 4.3; p < 0.001), B cells (2.7 ± 2.0 vs. 4.6 ± 3.3; p = 0.024), and dendritic cells (3.0 ± 0.9 vs. 3.8 ± 1.0, p = 0.007), and a higher percentage of endothelial cells (4.8 ± 0.7 vs 4.2 ± 0.6, p = 0.004) and adipocyte cells (2.6 ± 0.5 vs. 2.4 ± 0.5; p = 0.005) (Fig. 5A, Additional file 4: Table E3). There were no significant fractional differences in other cell types, including CD8 + T cells or neutrophils, between the two subclasses (Additional file 4: Table E3). We considered that in addition to a decrease in T cell populations, there may also be differences in the T cell repertoire between the subclasses, which can be a surrogate for monitoring the effectiveness of the adaptive immune system. Using ImReP, we found that compared to Subclass 2, Subclass 1 had significantly fewer T cell clonotypes (Fig. 5B, Additional file 4: Table E4) and lower activities and diversity of TCRA, TCRB, and TCRD. (Fig. 5C, D, Additional file 4: Table E5).

Fig. 5

A Box plot depicting differences in whole blood cell type abundances between Subclass 1 and 2. Subclass 1 has significantly lower percentages of CD4 T cells (0.9 ± 1.5 vs. 5.2 ± 4.3; p < 0.001) and B cells (2.7 ± 2.0 vs. 4.6 ± 3.3; p = 0.024), and a higher percentage of endothelial cells compared to Subclass 2 (4.8 ± 0.7 vs. 4.2 ± 0.6). There were no significant differences in percentages of CD8 T cells between the two subclasses. B Number of TCR reads per 1 million RNA-Seq reads in septic shock subclasses. Subclass 1 had significantly fewer TCR reads normalized per RNA sequencing reads than Subclass 2 (3.91 ± 4.30 vs. 10/03 ± 8.46; p = 0.001). C Subclass 1 had significantly fewer mean clonotypes than Subclass 2 (59.70 ± 59.49 vs. 155.96 ± 138.4; p = 0.002). D Inverse Simpson index, a measure of mean diversity of the T cell repertoire, between Subclass 1 versus Subclass 2. Subclass 1 had a less diverse mean T cell repertoire than Subclass 2 on the Inverse Simpson index (43.84 ± 45.11 vs. 114.45 ± 108.8; p = 0.002)

Discussion

We provide important insights into the pathophysiology of septic shock and the heterogeneity of biological perturbations within patients with septic shock. We compared the expression profiles of children with septic shock to mechanically ventilated controls and found upregulated expression of innate immune and neutrophil pathways and downregulation of adaptive immune pathways in children with septic shock. On further analysis, the expression profiles of patients with septic shock clustered around two subclasses with differential upregulation of innate immunity and downregulation of adaptive immunity. While both subclasses were clinically identified as “septic shock,” Subclass 1 was characterized by upregulation of innate immunity and biomarker profiles consistent with a hyperimmune response along with concomitant downregulation of both B and T cells and lesser T cell receptor diversity. This same subclass was associated with worse clinical outcomes, underscoring the importance of addressing sepsis heterogeneity. Subclass 2, while having the same clinical diagnosis, appeared biologically much more similar to our critically ill controls.

Our study suggests that patients with severe initial perturbations of the innate and adaptive immune systems, such as those in Subclass 1, are likely to have worse outcomes including higher PELOD scores, higher maximum inotrope usage, and longer intensive care unit and hospital stays, despite the lack of significant differences in illness severity or the pediatric risk of mortality (PRISM) score on initial presentation, indicating that conventional clinical scores alone are often are not sufficient to capture the pathophysiological complexity or outcomes of septic shock. An interesting clinical difference in our cohort is that Subclass 1 patients tended to be older and had modestly higher average blood glucose levels at baseline. [23]. IL-6 is known to contribute to hyperglycemia through insulin resistance, and multiple retrospective analyses have reported that hyperglycemia is associated with adverse outcomes [24].

Not only did Subclass 1 have upregulated genes related to neutrophil immunity compared to Subclass 2, but these patients also had significantly higher levels of proinflammatory cytokines (PAI1, IL6, IL8, ANG2, TREM1, and TNFR1), and a greater degree of organ dysfunction. Early unchecked innate immune-driven inflammation has been associated with a more profound degree of organ injury. Concomitantly, sepsis has also been shown to cause adaptive immunosuppression with a marked loss of T and B cells [25]. Subclass 1 had significant downregulation of genes related to T and B cell activity, and a lower percentage of CD4 T cells and B cells compared to Subclass 2. Severe T cell dysfunction, leading to decreased T cell cytotoxicity and T cell apoptosis in the setting of a high antigen load and elevated cytokines, is known to occur in septic shock [25].It is a source of debate whether it is the innate immune-driven hyperinflammation, adaptive immunosuppression, or a contribution of both is the driver of morbidity and mortality in sepsis. A recent study supported the latter, as they demonstrated that during sepsis, the proliferation of a large population of immature neutrophils inhibited the proliferation and activation of CD4 + T cells, and a subset of patients with higher frequencies of the immature neutrophils had poorer outcomes [26]. A potential consideration here is that the different sepsis subclasses could represent different chronological stages of the same underlying pathobiology. We attempted to mitigate this confounder by excluding patients who initiated vasoactive support > 72 h prior to enrollment.

Comparison to sepsis subclasses in the literature

In the adult population, at least five independent research groups have identified sepsis subgroups describing an adaptive immunity suppression phenotype with corresponding higher mortality such as Scicluna et al.’s MARS1 endotype [27], Davenport et al.’s SRS1 subgroup [28], and Sweeney et al.’s Inflammopathic subtype [29,30,31].

In the pediatric population, Hector Wong et. al. were the first to identify two pediatric sepsis endotypes, A and B. Endotype A was characterized by upregulation of innate immunity pathways and repression of pathways related to the adaptive immune system and glucocorticoid receptor signaling. However, the gene expression pattern that differentiates adult SRS groups was not enriched in the pediatric endotypes [32, 33]. This suggests that there may be differences between adults and children, highlighting the importance of studying sepsis specifically in the pediatric population.

Our findings and subclassification are consistent with previously published subclasses in the pediatric sepsis population. Our Subclass 1 was analogous to the Wong et al. pediatric Endotype A, which was identified in a separate pediatric cohort. Similar to Endotype A in Dr Wong’s cohort, Subclass 1 was also characterized by repression of the adaptive immune system related to T cell activation and had worse clinical outcomes than Subclass 2. We compared our gene list with the previously validated 100 gene signature that differentiated Endotype A and B, and found 12 of the 100 genes present in our list (GNAI3, PLCG1, CD3E, CD247, NCR3, ARPC5, ZAP70, FYN, SEMA6B, TLR8, CAMK2D, TLR8) [34]. In 2021, Muszynski et al. also identified a subclass of septic shock children with immunoparalysis with worse clinical outcomes [35]. They found that this subclass’s transcriptomic profiles demonstrated upregulated pathways in leukocyte extravasation, and downregulation in adaptive immunity pathways [35]. When comparing our Subclass 1 DEGs to Musznski’s immunoparalysis subclass, shared upregulated DEGs included COL17A1, LAMA2, FFAR3, RAP1GAP, XCR1, MMP27, and MMP8, and shared downregulated DEGs included KLRC1 and IL2RB. Thus, while Wong’s Endotype A, Muszynski’s immunoparalysis subclass, and our Subclass 1 are similar, challenges remain in identifying an appropriate panel of candidate genes out of the often large lists of DEGs that can be generalized across the pediatric population. More research in different pediatric cohorts is needed to develop a consensus subclassification system.

At the molecular level, T cell diversity plays a key role in the ability of the adaptive immune system to effectively mount a response to invading pathogens. Each TCR is made up of alpha (TCRA) and beta chains (TCRB) or delta (TCRD) and gamma (TCRG) chains. TCR diversity is generated through V(D)J recombination in the early stages of T cell maturation in the thymus and is critical to effectively recognizing antigen peptides. In the adult sepsis population, studies have shown that septic patients present with a marked decrease in TCR diversity after the onset of shock, which is associated with mortality and the development of nosocomial infections [36, 37]. However, there is a paucity of literature looking at TCR diversity in pediatric sepsis, and one study demonstrated that the TCR repertoires in adults and children are discrepant and thus difficult to directly compare [38]. Our study demonstrated that Subclass 1 patients demonstrated reduced TCR diversity and clonality, including decreased diversity of TCRA, TCRB, and TCRD, with associated worse clinical outcomes at the onset of septic shock.

In the adult population, there is evidence that in sepsis, a decrease in circulating B cells is associated with a poor prognosis [39]. Additionally, adult studies suggest that the B cell depletion is selective, and IL-10-producing B cells may actually increase and exacerbate immunosuppression [39]. In this pediatric study, Subclass 1 had decreased levels of B cells, but significantly higher levels of IL-10 with worse outcomes, similar to adult findings. The mechanism behind the depletion of B cells in sepsis is not well understood. One theory is that profound sepsis can impair bone marrow production leading to decreases in B cell numbers. Alternatively, some studies suggest sepsis signals can possibly trigger B cell apoptosis [39]. More studies are needed to investigate this phenomenon.

Our study is novel, as for the first time it not only provides a comprehensive evaluation of dysregulated pathways based on genome-wide differential gene expression but also enhances it with cell deconvolution and T and B cell receptor diversity estimations in the context of plasma biomarkers and clinical characteristics in septic shock to better characterize and sub-phenotype the biological perturbations in septic shock that are biologically plausible and clinically relevant in the pediatric population. A potential future direction of this subclassification study would be to create a classifier model that could eventually become useful in the clinical setting to classify pediatric septic shock patients and inform clinical decisions.

Limitations of the study

First, given that this was a retrospective analysis of a prospective pediatric clinical trial, a limitation was that a suitable case and control population were identified post-retrospectively. We assumed accurate and timely documentation of sepsis and initiation of inotropes for cases, and clinician-documented reasons for PICU admission and initiation of mechanical ventilation and negative blood cultures to create a control group. However even in the unlikely scenario that an underlying infection was missed in our controls, our cases required inotropes, indicating that on the clinical spectrum, sepsis in our cases was certainly much more severe. The use of mechanically ventilated controls ensured that our results were specific to septic shock and not just a result of generalized critical illness. Another limitation was that we could not differentiate between the different sources of sepsis each patient may have had. Finally, our sample sizes were relatively small. A larger sample size would have increased statistical power to detect clinically meaningful differences within each subclass; therefore, mortality was not a reported outcome in this study given that only three of the patients died within our cohort. Therefore, while Subclass 1 had worse outcomes, we cannot state that they had higher mortality.

Conclusions

We identified two subclasses of children with septic shock based on differential gene expression using RNA-Seq. These two subclasses have differential regulation of genes related to the immune system that is relevant to the pathophysiology of sepsis and septic shock. One subclass is characterized by upregulation of innate immunity pathways and repression of adaptive immunity pathways, with lower levels of T cells and B cells. This subclass is associated with clinically worse outcomes. Thus, subclassifying patients with septic shock based on genome-wide expression data can aid in identifying both new targets for therapies and which patients would likely benefit from them.