Dysregulated STAT3 signaling and T cell immunometabolic dysfunction define a targetable, high mortality subphenotype of critically ill children =============================================================================================================================================== * Robert B. Lindell * Samir Sayed * Jose S. Campos * Montana Knight * Andrea A. Mauracher * Ceire A. Hay * Peyton E. Conrey * Julie C. Fitzgerald * Nadir Yehya * Stephen T. Famularo III * Teresa Arroyo * Richard Tustin III * Hossein Fazelinia * Edward M. Behrens * David T. Teachey * Alexandra F. Freeman * Jenna R. E. Bergerson * Steven M. Holland * Jennifer W. Leiding * Scott L. Weiss * Mark W. Hall * Athena F. Zuppa * Deanne M. Taylor * Rui Feng * E. John Wherry * Nuala J. Meyer * Sarah E. Henrickson ## ABSTRACT Sepsis is the leading cause of death of hospitalized children worldwide. Despite the established link between immune dysregulation and mortality in pediatric sepsis, it remains unclear which host immune factors contribute causally to adverse sepsis outcomes. Identifying modifiable pathobiology is an essential first step to successful translation of biologic insights into precision therapeutics. We designed a prospective, longitudinal cohort study of 88 critically ill pediatric patients with multiple organ dysfunction syndrome (MODS), including patients with and without sepsis, to define subphenotypes associated with targetable mechanisms of immune dysregulation. We first assessed plasma proteomic profiles and identified shared features of immune dysregulation in MODS patients with and without sepsis. We then employed consensus clustering to define three subphenotypes based on protein expression at disease onset and identified a strong association between subphenotype and clinical outcome. We next identified differences in immune cell frequency and activation state by MODS subphenotype and determined the association between hyperinflammatory pathway activation and cellular immunophenotype. Using single cell transcriptomics, we demonstrated STAT3 hyperactivation in lymphocytes from the sickest MODS subgroup and then identified an association between STAT3 hyperactivation and T cell immunometabolic dysregulation. Finally, we compared proteomics findings between patients with MODS and patients with inborn errors of immunity that amplify cytokine signaling pathways to further assess the impact of STAT3 hyperactivation in the most severe patients with MODS. Overall, these results identify a potentially pathologic and targetable role for STAT3 hyperactivation in a subset of pediatric patients with MODS who have high severity of illness and poor prognosis. ## INTRODUCTION Sepsis is defined as life-threatening organ dysfunction that develops in the setting of immune dysregulation (*1*) and represents a leading cause of adult and pediatric mortality worldwide (*2*). Pediatric sepsis is associated with distinct epidemiology and outcomes compared to adult sepsis (*3*), with the majority of pediatric sepsis deaths occurring in immunocompromised children (*4–6*) and in the setting of persistent multiple organ dysfunction syndrome (MODS) (*7, 8*). MODS represents a “final common pathway” by which severe, systemic inflammation arising from diverse clinical insults (e.g. sepsis, trauma, shock) leads to progressive organ failure due to a combination of endothelial, epithelial, mitochondrial, and immunologic dysfunction (*9–11*); in children, sepsis is the most common cause of MODS (*12–14*). Identifying modifiable pathobiology is an essential first step to successful translation of core biologic insights into precision therapeutics for critical illness (*15*). In pediatric sepsis, both innate and adaptive immune dysfunction (*16–18*) are associated with secondary infection (*19*), persistent organ dysfunction (*20*) and mortality (*21*). Mitochondrial dysfunction, a hallmark of this sepsis-associated immune suppression (*22, 23*), is also associated with organ failure in pediatric sepsis (*24–26*). Yet despite the established link between immune dysregulation and both morbidity and mortality in pediatric sepsis, successful immune modulation in sepsis has been elusive (*27*) and it remains unclear which host immune factors contribute causally to adverse sepsis outcomes (*28*). A successful precision medicine approach to pediatric sepsis will require an understanding of the molecular events which underlie the development of organ failure, a knowledge of which events are reversible, and an ability to identify high-risk patients in real time (*29*). Existing approaches to identify pediatric sepsis subphenotypes based on clinical and laboratory features have identified several high-risk patient subgroups (*30–33*), but these subgroups lack defined mechanisms of immune dysregulation and thus cannot inform treatment decisions. A wide range of pathogenic inborn errors of immunity (IEI) have been identified in cohorts of pediatric patients with sepsis (*34*), suggesting that diverse genetic factors may contribute to the dysregulated host immune phenotypes that develop in the setting of critical illness. Recent attempts to define this heterogeneous host response in adults using molecular approaches (*35–41*) have yielded new mechanistic insights into sepsis pathobiology, suggesting that deep immune phenotyping may be able to overcome the immunologic heterogeneity that has previously hampered precision sepsis efforts. Given the role of immune dysregulation in sepsis pathobiology and the lack of novel targeted therapies for these high-risk patients, we designed a longitudinal cohort study of pediatric MODS patients with and without sepsis with the goal of identifying severity-associated MODS subphenotypes. We hypothesized that severity-associated MODS subphenotypes would be associated with targetable mechanisms of immune dysregulation. Using longitudinal proteomics, cytometry, and transcriptional analysis, we defined three prognostic subphenotypes in critically ill children with and without sepsis and identified an association between pathologic STAT3 hyperactivation and T cell immunometabolic dysregulation in the subset of pediatric MODS patients with the highest severity of illness and worst prognosis. ## RESULTS ### Proteomic heterogeneity of pediatric MODS To identify MODS subphenotypes in patients with and without sepsis, we collected peripheral whole blood samples in a prospective, observational cohort of 88 pediatric patients with MODS. Blood was collected at MODS onset and then twice weekly through death or resolution of MODS. At each timepoint, peripheral blood mononuclear cells (PBMC) and heparin plasma biospecimens were cryopreserved for later analysis, as shown in **Fig. 1A**. We compared these patients to a separate cohort of 25 pediatric healthy control (HC) participants. Our initial analyses included a 1536-marker proteomics panel from Olink Proteomics selected for patient stratification and identification of candidate therapeutic targets (*42*) (**Table S1**) and a 35-marker spectral flow cytometry-based immune phenotyping panel optimized to measure PBMC abundance and activation state (**Table S2**). Demographics and clinical outcomes were similar between our 35 MODS patients with sepsis and 53 MODS patients without sepsis (**Fig. 1B**), and the only organ dysfunction category that differed between groups was hematologic dysfunction, which was driven by increased thrombocytopenia among patients with sepsis at MODS onset (platelet count median [IQR]: 119 [99-187] vs 199 [131-278], *p*<0.001; normal range 150-450 x103/μl). Length of stay, cumulative PELOD-2 (Pediatric Logistic Organ Dysfunction-2) (*13*) organ dysfunction score, and survival to PICU (pediatric intensive care unit) discharge did not differ between MODS subgroups. Sepsis is characterized by a surge in pro-inflammatory cytokines, and while we noted the expected increased proinflammatory cytokines (interleukin (IL)-1β, IL-6, tumor necrosis factor (TNF), IL-18, monocyte chemoattractant protein (MCP)-1, interferon (IFN)-γ, all p<0.0001) in patients with sepsis compared to HC participants, many of these cytokines were also markedly elevated in MODS patients without sepsis. We found that mean IL-1β, TNF, and IL-6 levels did not differ between patients with and without sepsis, while IFN-γ was higher in patients with sepsis (*p*<0.0001) (**Fig. 1C**). We then used principal component analysis to visualize the heterogeneity within the MODS cohort and overlap between patients with and without sepsis. We noted separation between HC participants and MODS patients but substantial overlap between MODS patients with and without sepsis (**Fig. 1D**). Patients with MODS also exhibited substantial heterogeneity along principal component 1 and principal component 2, dimensions which are strongly associated with proteins corresponding to shock, growth regulation, and inflammatory signaling (**Table S3**). This biologic overlap suggests that MODS develops in the setting of overwhelming proinflammatory signals and may represent a final common pathway of critical illness, independent of the inciting diagnosis. To better understand the shared proteomic landscape of MODS patients with and without sepsis, we constructed a heatmap of row-normalized protein expression with the full proteomics dataset and used hierarchical clustering to visualize heterogeneity within the cohort at MODS onset. As shown in **Fig. 1E**, hierarchical clustering separates HC participants from patients with MODS but does not discriminate between MODS patients with sepsis and without sepsis, findings which are complementary to the principal component analysis. ![Fig 1.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2024/06/11/2024.06.11.24308709/F1.medium.gif) [Fig 1.](http://medrxiv.org/content/early/2024/06/11/2024.06.11.24308709/F1) Fig 1. Shared proteomic features in pediatric MODS patients with and without sepsis. **(A)** Schematic of participant cohorts and sample processing pipeline using plasma proteomics and spectral flow cytometry. **(B)** Cohort demographics, exposures, and clinical outcomes, stratified by MODS etiology (i.e. sepsis vs non-sepsis). Comparisons between groups using chi-squared and Wilcoxon rank sum test, as appropriate. **(C)** Proinflammatory cytokines are increased in patients with MODS but do not differ between patients with and without sepsis, except for IFN-γ which is increased in patients with sepsis compared to patients without sepsis. Comparisons were made between groups as shown and determined by Wilcoxon rank sum test. Box-and-whiskers plots include a box indicating median and interquartile range, and whiskers extending to the data point not further than the 1.5x interquartile range. **(D)** Principal component analysis based on expression of 1,448 proteins analyzed by proximity extension assay in 186 samples from 88 patients with MODS and 25 pediatric healthy control (HC) participants reveals substantial heterogeneity and overlap between MODS patients with and without sepsis. **(E)** Clustered heatmap showing row-normalized protein expression among 88 patients at MODS onset and 25 HC participants. Colored annotation bar indicates MODS subgroup. ### Identification of three severity-associated MODS subphenotypes Given the overlapping proteomic profile between MODS patients with and without sepsis, we sought to identify MODS subphenotypes based on plasma protein expression at MODS onset. First, we used linear mixed-effects models to identify plasma proteins associated with severity of illness, defined by the PELOD-2 organ dysfunction score, after adjustment for age and sex as fixed effects and day from MODS onset as a random effect. This yielded a set of 214 plasma proteins which were significantly associated with illness severity (**Table S4**). We then employed consensus clustering – an unsupervised machine learning methodology for subclass discovery (*43*) – to identify distinct clusters of patients with MODS based on protein expression (**Fig. 2A**). Optimal cluster number (*k*=3) was identified via Monte Carlo bootstrapping and confirmed by elbow method and gap statistic (**Fig. S1**). To understand the association between cluster membership and clinical outcomes, we tested the effect of MODS subphenotype on the cumulative incidence of both mortality and survival to PICU discharge with the Fine-Gray subdistribution hazard model (*44*). In this competing risk survival analysis, patients in Group C have higher cumulative incidence of death (*p*=0.03) and lower cumulative incidence of survival to PICU discharge (*p*=0.04) compared to patients in Clusters A and B, with separation occurring in the first week and persisting to day +28 (**Fig. 2B**). Though defined by protein expression alone, the identified MODS subphenotypes also differ by clinical features and outcomes (**Fig. 2C**). MODS etiology varies by Group, with most sepsis patients in Group B and C and all trauma patients in Group A. Group C patients have higher severity of illness, more non-cardiopulmonary organ failures, and increased cumulative organ dysfunction scores and mortality compared to Groups A/B. Etiology of MODS and computed subphenotype for each patient is detailed in **Table S5**. Immunocompromised diagnoses were identified in 13% of MODS patients (11/88), and patients with immunocompromised status were not imbalanced across MODS subphenotypes (p=0.74), as detailed in **Table S6**. Noting that these protein-derived MODS subphenotypes are associated with an ordinal increase in number of organ failures and mortality across subphenotypes, we hypothesized a reduced set of plasma proteins could successfully identify MODS subphenotypes and could thus be more suitable for translation to the clinical setting. Elastic net regularization is a common approach to generate a high-performing sparse model with good prediction accuracy (*45*). To define a parsimonious protein signature, we trained an ordinal elastic net model (*46*) and generated a 24-protein signature which successfully discriminates the three subphenotypes (**Fig. S2**). We next sought to evaluate the performance of this parsimonious protein model through two complementary approaches. First, we tested the association between the elastic net proteins and severity of illness. Using a linear mixed-effects model, we estimated the fold change in protein expression associated with one standard deviation change in PELOD-2 organ dysfunction score (**Fig. 2D**). This model demonstrates that modest changes in expression of these 24 proteins are associated with meaningful differences in illness severity. Second, we tested the discrimination of our 24-protein signature for MODS subphenotype. As shown in the heatmap of protein expression at MODS onset in **Fig. 2E**, hierarchical clustering of elastic net proteins effectively separates MODS subphenotypes but does not discriminate between MODS patients with and without sepsis. To quantify the discrimination of the parsimonious elastic net protein set for MODS subphenotypes, we used linear discriminant analysis to calculate the polytomous discrimination index (PDI), a measure of rank-based discrimination performance, for each MODS subphenotype. Category-specific PDI indicated excellent discrimination for each subphenotype (Group A 0.98, Group B 0.96, Group C 0.99), and overall PDI for the model was 0.98. Taken together, these results suggest that our reduced 24-protein signature reflects severity of illness and successfully discriminates subphenotypes at MODS onset. ![Fig. 2.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2024/06/11/2024.06.11.24308709/F2.medium.gif) [Fig. 2.](http://medrxiv.org/content/early/2024/06/11/2024.06.11.24308709/F2) Fig. 2. Three MODS subphenotypes based on protein expression at MODS onset. **(A)** Principal component analysis based on expression of 214 severity-associated proteins in 88 patients at MODS onset reveals separation among three MODS subphenotypes defined through Monte Carlo consensus clustering (*43*). **(B)** Cumulative incidence of survival to PICU discharge and PICU mortality by MODS day among three MODS subphenotypes. **(C)** MODS etiology, cohort demographics, exposures, and clinical outcomes, stratified by MODS subphenotype. Comparisons between groups by chi-squared and Wilcoxon rank sum test, as appropriate. **(D)** Forest plot showing estimated fold-change in normalized protein expression associated with a +1 standard deviation increase in PELOD-2 score for 24 proteins identified through ordinal elastic net regression (*46*) which differ across MODS subphenotypes. Forest plot includes a point estimate of fold-change and whiskers which represent the 95% confidence interval around this estimate. **(E)** Heatmap showing normalized protein expression of the 24 proteins identified by ordinal elastic net regression among 88 patients at MODS onset. Colored annotation bars indicate MODS subphenotype and sepsis status. ### Immune cell frequency and activation vary by MODS subphenotype Having identified three MODS subphenotypes (Groups A, B, and C) based on severity-associated protein expression and developed a parsimonious model to classify patients with MODS, we hypothesized that protein-derived subphenotypes would be associated with differences in cellular immunophenotype. For this analysis, we performed high dimensional spectral flow cytometry on cryopreserved PBMCs obtained at MODS onset using a custom-designed 35 marker panel (**Table S2**) which includes both phenotypic and functional markers. After arcsinh scaling (*47*) and quality control with flowAI (*48*), we performed FlowSOM metaclustering (*49*) and identified 14 PBMC populations by surface and intracellular marker expression. **Fig. S3** presents a representative example of our manual gating strategy for this immune phenotyping panel, which we used to confirm the identity of the FlowSOM metaclusters. To visualize immunophenotypic differences between HC participants and the three MODS subphenotypes (**Fig. 3A**), we subsampled the data to 100,000 cells per MODS subphenotype and applied t-distributed Stochastic Neighbor Embedding with Compute Unified Device Architecture (tSNE-CUDA) dimensionality reduction (*50*). Differences in cell populations by MODS subphenotype are quantified in stacked bar plots (**Fig. 3B**), showing ordinal reduction of lymphoid and myeloid effector cells across MODS groups. The proportional abundance of non-naïve (central memory, effector memory, and terminally differentiated) CD8+ T cells were markedly reduced in Group C patients (**Fig. 3C**), and this loss of non-naïve CD8+ T cells was strongly associated with severity of illness by linear regression (*p*<0.001). A similar ordinal trajectory was noted across many cell types (**Fig. 3D**), including T cells, B cells, NK cells, and dendritic cells (Cuzik test of trend *p*<0.001 for each cell type shown). As seen with the reduced frequency of non-naïve T cells, reduction in frequency of each of these cell types was strongly associated with severity of illness by linear regression (all *p*<0.001). In addition to shifts in cellular frequency, we hypothesized that protein-derived subphenotypes would be associated with differences in immune cell activation, as measured by expression of markers of proliferation (Ki67) and activation (CD38 and HLA-DR), across MODS subphenotypes. Representative bivariate plots of Ki67 expression in non-naïve CD4+ and CD8+ T cells and cytotoxic NK cells demonstrate increased proliferation in Group B and Group C patients compared to Group A patients (**Fig 3E**). Similarly, CD38 and HLA-DR co-expression in non-naïve CD8+ T cells was markedly increased in Group B and Group C patients compared to Group A, indicative of CD8+ T cell activation in these groups (**Fig. 3F**). Corresponding boxplots in **Fig. 3G** quantify these differences in proliferation and activation by MODS subphenotype. Taken together, these data suggest that CD8+ T cells have markedly reduced abundance but concurrently have increased markers of proliferation and activation in Group C patients, and that these cellular phenotypes are associated with worsening severity of illness. ![Fig. 3.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2024/06/11/2024.06.11.24308709/F3.medium.gif) [Fig. 3.](http://medrxiv.org/content/early/2024/06/11/2024.06.11.24308709/F3) Fig. 3. PBMC effector cell abundance and activation state vary by plasma protein-derived MODS subphenotype. **(A)** PBMC immune phenotype tSNE-CUDA projection demonstrates marked differences in frequency of lymphoid and myeloid lineages between HC participants and patients with MODS, and among MODS subphenotypes. Differential abundance across groups by cell subtype was determined by Kruskal-Wallis test. **(B)** Stacked bar plots show differences in CD3+ and CD3-cell abundance by MODS subphenotype corresponding to the FlowSOM metaclusters in the tSNE-CUDA projections at left. **(C)** Compared to HC participants, the abundance of CD8+ T cell subsets are reduced in patients with MODS, with an ordinal trajectory across severity subgroups noted in CD8+ central memory cells, CD8+ effector memory cells, and CD8+ Temra cells. Effector CD8+ cells in Group C patients are significantly reduced in frequency compared to other MODS subgroups. Comparisons were made between groups as shown and determined by Wilcoxon rank sum test. **(D)** Shifts in other lymphoid and myeloid subsets across MODS subgroups demonstrate a loss of peripheral non-naïve CD4+ and CD8+ T cells, γδ-T cells, memory B cells, cytotoxic NK cells, and plasmacytoid dendritic cells, with the most marked differences in the Group C patients with the highest disease severity. Comparisons were made between groups as shown and determined by Wilcoxon rank sum test. **(E)** In non-naïve CD4+, non-naïve CD8+, and cytotoxic NK cells, Ki67 expression is increased Group B and Group C patients compared to Group A (all *p*<0.05), indicative of increased proliferation of these effector cell subgroups. **(F)** In non-naïve CD8+ T cells, CD38 and HLA-DR co-expression is markedly increased Group B and Group C patients compared to Group A (both *p*<0.01), indicative of CD8+ T cell activation in these groups. **(G)** Boxplots quantifying differences in proliferation and activation by MODS subphenotype. Comparisons were made between groups as shown and determined by Wilcoxon rank sum test. ### Multiple concurrent mechanisms of immune dysregulation define Group C patients Because Group C patients have distinct clinical, proteomic, and cellular features compared to other patients with MODS, we hypothesized that expression of key canonical inflammatory pathways would differ by MODS subphenotype. Using the proteomics dataset, we first examined differential protein expression after adjustment for patient age, sex, severity of illness, and days since MODS onset using a linear mixed-effects model. In unadjusted analysis, 1061/1448 measured proteins were differentially expressed in Group C (**Fig. 4A**), and in our adjusted model 1003/1061 proteins remained differentially expressed. For 98% (980/1003) of these differentially expressed proteins, expression was upregulated in Group C compared to Group A/B. Using the adjusted protein expression dataset, we then performed pathway enrichment analysis of our plasma proteomics data using complementary group-and patient-based strategies. First, we used Ingenuity Pathway Analysis (IPA, Qiagen) (*51*) to identify enrichment of 22 canonical pathways in the Group C proteome in comparison to Group A/B (**Fig. 4B**). Noting that 8 of the top 10 differentially expressed pathways identified by IPA were related to hyperinflammatory signaling, we then applied Gene Set Variation Analysis (GSVA) (*52*) to study enrichment of five canonical proinflammatory pathways on the individual patient level using Human Molecular Signatures Database (MSigDB) hallmark gene sets (*53*). Using GSVA, we assigned patient-level protein module enrichment scores for each pathway of interest based on expression of proteins in the corresponding Hallmark gene set, as detailed in Methods. Module enrichment scores for each pathway are visualized in Fig. 4C and demonstrate markedly increased IL-6/JAK/STAT3 and IL-2/JAK/STAT5 module enrichment in Group B and Group C patients (each *p*<0.001 vs Group A), modest increase in TNF/NFkB signaling in Group B and Group C patients (each *p*<0.05 vs Group A), and an ordinal reduction in PI3K/AKT/MTOR signaling across MODS subphenotypes (Cuzik test of trend *p*<0.001). Conversely, IFN-γ response module enrichment scores were minimally increased in patients with MODS and did not vary by MODS subphenotype. We then assessed the correlation between Hallmark module enrichment and PBMC immune cell subset abundance and activation in patients with MODS (**Fig. 4D**). We noted that IL-6/JAK/STAT3 module enrichment had the strongest correlation with proliferation and activation of PBMC subsets previously identified in **Fig. 3E-F**, including non-naïve CD4+ and CD8+ proliferation (Ki67 expression) and activation (HLA-DR/CD38 coexpression and PD-1/CD39 coexpression). Classical and non-classical monocyte HLA-DR expression was inversely correlated with IL-6/JAK/STAT3 module enrichment score, while immune regulatory cell populations (regulatory T cells, myeloid-derived suppressor cells) were positively correlated with IL-6/JAK/STAT3 module enrichment score. Finally, we studied the longitudinal expression of canonical cytokine storm markers (*54*) in patients with MODS to understand the duration of immune dysregulation in Group C patients. Normalized protein expression by day since MODS onset by subphenotype is shown in **Fig. 4E**. We noted that IL-6, MCP-1, and IL-18 expression are significantly higher in Group C patients for the first 7 days after MODS onset (each *p*<0.001 at day +7), while IL-1β and TNF expression remain different for only the first 4 days (*p*=0.05 and *p*=0.009 respectively at day +4) and IFN-γ does not differentiate MODS subphenotypes. ![Fig. 4.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2024/06/11/2024.06.11.24308709/F4.medium.gif) [Fig. 4.](http://medrxiv.org/content/early/2024/06/11/2024.06.11.24308709/F4) Fig. 4. Group C patients have a hyperinflammatory plasma protein phenotype characterized by multiple concurrent mechanisms of immune dysregulation. **(A)** Volcano plot comparing proteomic profile of Group C patients to Group A/B patients, showing that 1061/1448 measured proteins were differentially expressed in Group C. **(B)** After correction for patient age, sex, severity of illness and days since MODS onset with a linear mixed-effects model, 1003 proteins remained differentially expressed in group C. Based on these differentially expressed proteins, Ingenuity Pathway Analysis was used to identify enrichment of 22 canonical pathways in the Group C proteome after Benjamini-Hochberg correction. **(C)** Individual patient module enrichment scores calculated using GSVA and MSigDB Hallmark pathways demonstrate markedly increased IL-6/JAK/STAT3 and IL-2/JAK/STAT5 signaling in Group B and Group C patients, modest increase in TNF/NFkB signaling in Group B and Group C patients, and an ordinal reduction in PI3K/AKT/MTOR signaling across MODS subphenotypes. Comparisons were made between groups as shown and determined by Wilcoxon rank sum test. **(D)** Spearman’s correlation between GSVA module enrichment score and PBMC abundance, activation, and proliferation within MODS patients. Correlations with Benjamini-Hochberg corrected *p*-value <0.05 are shown. **(E)** Longitudinal expression of canonical cytokine storm markers demonstrates persistence of immune dysregulation beyond MODS onset. For each panel, longitudinal normalized protein expression is presented by days from MODS onset through day 14. Loess regression lines with 95% confidence intervals are presented for each MODS subphenotype. ### Cell-specific signatures of immunometabolic dysregulation in patients with STAT3 hyperactivation Having identified a link between STAT3 pathway signaling, severity of illness, relative decrease in frequency of effector cells populations, and increased T cell proliferation and activation by flow cytometry, we next hypothesized that STAT3 signaling would have cell-specific effects in Group C patients. To assess the impact of STAT3 hyperactivation on immune cell phenotype and function in patients with MODS, we conducted a single cell RNA sequencing experiment (10X Genomics, Pleasanton, CA) on cryopreserved PBMCs obtained at MODS onset for Group C patients (n=9) and pediatric HC participants (n=3). IL-6/JAK/STAT3 module enrichment scores were significantly different between the 9 patients with MODS and 3 HC participants (0.331 vs. −0.494, *p*<0.0001) selected for the transcriptomics experiment. **Fig. 5A** shows the study schematic for this analysis. We sorted live CD45+ PBMCs using a Cytek Aurora CS cytometer (Cytek Biosciences, Fremont, CA) and then profiled 10,000 cells per patient using a 5’ RNA tag single cell sequencing approach (scRNA-seq; 10X Genomics, Pleasanton, CA). Libraries were sequenced to a depth of ∼30,000 reads per cell on a Novaseq S2 (Illumina, San Diego, CA). Transcripts were aligned to the GRCh38 reference genome using Cell Ranger v8.0 (10X Genomics, Pleasanton, CA). After quality control and integration, cell identities were inferred using the Azimuth reference-based mapping pipeline (*55*) and used to define 14 PBMC populations by transcriptional profile. To visualize immunophenotypic differences between HC participants and Group C patients in the scRNA-seq experiment (**Fig. 5B**), we subsampled the data to 30,000 cells per group and applied Uniform Manifold Approximation and Projection (UMAP) dimensionality reduction (*56*). Differences in cell population abundance are quantified in stacked bar plots (**Fig. 5C**), and demonstrate shifts in CD4+ lymphocytes, CD8+ lymphocytes, B cells, NK cells, and monocytes abundance between HC participants and Group C (all *p*<0.001). Compared to the flow cytometry immune profiling assay results, our transcriptomics analysis revealed a higher proportion of non-naïve CD4+ T cells, but other cell subset frequencies are similar between modalities. These shifts in differential abundance between flow cytometry and transcriptional data may reflect inherent differences in cell subset identification using transcriptional and cell surface markers (*57*). Consistent the flow cytometry data set, the transcriptomic data set reflects the loss of effector cell subsets in Group C patients compared to HC participants, including loss of non-naïve CD8+ T cells and cytotoxic NK cells. To understand cell-specific differences in STAT3 signaling in patients with MODS, we performed pseudobulk analysis of differential pathway expression using UCell (*58*). Pseudobulk analysis aggregates transcriptional profiles by cell subtype prior to statistical testing and has been shown to outperform other methods of differential expression analysis in single cell datasets (*59*). By pseudobulk analysis, we noted that CD4+, CD8+, and TCRgd+ T cells from patients with MODS had increased IL-6/JAK/STAT3 module enrichment scores compared to controls (all *p*<0.001) while monocytes, NK cells, and B cell IL-6/JAK/STAT3 module enrichment did not differ between groups, as shown in the split enrichment plots in **Fig. 5D**. Transcriptomic immune dysfunction scores have recently been applied to adult and pediatric sepsis, critical influenza, and COVID-19 for risk stratification (*36–40, 60, 61*). Among these, the quantitative sepsis response signature score (SRSq) is a prognostic metric summarizing the host transcriptional response to sepsis based on expression of 19 genes by bulk RNA sequencing (*40*). Because higher SRSq is associated with proinflammatory signaling, lymphocyte dysfunction, severity of illness, and higher risk of poor outcomes, we hypothesized that SRS gene enrichment would be associated with STAT3 signaling in Group C PBMC subsets. For each PBMC subset of interest, we generated a linear regression model to compare proinflammatory module enrichment to enrichment of the extended SRS gene set. The heatmap in **Fig. 5E** depicts the correlation between SRS and inflammatory module enrichment scores by cell type, and the adjusted R2 is displayed for all significant associations. While the STAT3 pathway had the highest degree of correlation to SRSq, the positive correlation was moderate and only applied to gene expression within the T cell compartment. STAT3 pathway activation has been shown to downregulate innate and adaptive immune responses in the tumor microenvironment (*62*), attenuating Th1 lymphocyte responses (*63*) and enhancing MDSC (*64*) and regulatory T cell function (*65*). STAT3 signaling also alters immune cell mitochondrial respiration, enhancing both oxidative phosphorylation (*66*) and glycolysis (*67*). Having identified loss of effector T cells and NK cells in our experimental data and increased STAT3 signaling in our scRNA-seq data, we next hypothesized that STAT3 hyperactivation would be associated with altered PBMC immunometabolic state. To investigate the link between STAT3 signaling and PBMC immunometabolism, we first calculated Hallmark glycolysis and oxidative phosphorylation pathway module enrichment scores using UCell (*58*) as detailed in Methods. We then constructed bivariate plots of glycolysis and oxidative phosphorylation module enrichment scores for both Group C and HC PBMCs (**Fig. 5F**). Compared to HC PBMCs, we identified two distinct populations of PBMCs in Group C – one with high module enrichment scores for both glycolysis and oxidative phosphorylation, and another with low module enrichment scores for glycolysis and oxidative phosphorylation. These populations are denoted on the bivariate plot with labeled boxes. As shown in the ridge plot in **Fig. 5F**, among cells with increased metabolic activity, we noted that IL-6/JAK/STAT3 module enrichment scores were lower than median pseudobulk values for each PBMC subset, suggesting that increased STAT3 signaling is associated with suppressed immunometabolic function. Based on this insight, we hypothesized that hyperactive STAT3 signaling would be associated with reduced effector cell metabolic activity. Because multiple inflammatory signaling pathways are concurrently active in Group C patients, we then constructed a mixed effects logistic regression model to determine the association between module enrichment score and immunometabolic profile of non-naïve CD4+, non-naïve CD8+, and cytotoxic NK cells in patients with MODS at the single cell level. As shown in **Fig. 5G**, IL-6/JAK/STAT3 and IL-2/JAK/STAT5 module scores at the single-cell level are associated with low glycolysis and oxidative phosphorylation module scores, while IFN-γ response and PI3K/AKT/MTOR module scores are associated with high glycolysis and oxidative phosphorylation module scores at the single-cell level (all *p*<0.01). These results suggest that STAT3 and STAT5 hyperactivation in T cells in patients with MODS are associated with poor immunometabolic function, which may partially explain the link between hypercytokinemia and T cell dysfunction in sepsis. Bulk PBMC mitochondrial function has been previously associated with adverse outcomes in pediatric sepsis (*24–26*), and these results suggest a mechanism by which this immunometabolic dysregulation may occur. ![Fig. 5.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2024/06/11/2024.06.11.24308709/F5.medium.gif) [Fig. 5.](http://medrxiv.org/content/early/2024/06/11/2024.06.11.24308709/F5) Fig. 5. Cell-specific signatures of immunometabolic dysregulation in MODS. **(A)** Study schematic for the transcriptomics experiment, which included 9 patients with MODS and 3 pediatric HC participants. Cells were sorted, sequenced, and analyzed using a bioinformatics pipeline to compare pathway enrichment scores by condition and cell type. **(B)** PBMC immune phenotype UMAP projection demonstrates differences in lymphoid and myeloid lineages between healthy and Group C patients. Differential abundance by cell subtype was determined by Wilcoxon rank sum test. **(C)** Stacked bar plots show differences in CD3+ and CD3-cell abundance by group corresponding to the metaclusters in the UMAP projections to the left. **(D)** Pseudobulk analysis of module enrichment scores across PBMC subsets showed that CD4+, CD8+, and TCRgd+ T cells from patients with MODS had increased IL-6/JAK/STAT3 module enrichment scores compared to controls while monocytes, NK cells, and B cell IL-6/JAK/STAT3 module enrichment scores did not differ between groups. **(E)** Heatmap depicting the correlation between SRSq and inflammatory pathway by cell type. Box color depicts the adjusted R2, and the R2 value is displayed for all significant associations. **(F)** Bivariate plots of single-cell glycolysis and oxidative phosphorylation enrichment scores for both Group C and HC PBMCs. Labeled boxes denote populations of increased and decreased immunometabolic pathway enrichment scores. At right, ridge plot denotes IL-6/JAK/STAT3 module enrichment by PBMC subset. **(G)** The estimated effect of inflammatory pathway signaling on immunometabolic activity of effector T and NK cells at the single cell level. Pathway beta coefficients and 95% confidence intervals are derived from mixed effects logistic regression model of the association between module enrichment score and immunometabolic profile of non-naïve CD4+, non-naïve CD8+, and cytotoxic NK cells. Positive coefficients are associated with increased glycolysis/oxidative phosphorylation activity while negative coefficients are associated with decreased glycolysis/oxidative phosphorylation activity. ### STAT3 as a candidate target for precision immunomodulation Having established increased plasma expression of STAT3 target proteins in Group C patients and that STAT3 module enrichment was associated with dysregulated cellular immunity and transcriptional evidence of immunometabolic dysregulation, we next sought to assess whether STAT3 activation represented a physiologic or pathologic response in the setting of critical illness. The STAT3 signaling pathway is a canonical inflammatory pathway associated with capillary leak, endothelial dysfunction, emergency granulopoiesis, and lymphocyte dysregulation (*68*). A recent human study has identified an association between immature CD66b+ neutrophils and STAT3-mediated emergency granulopoiesis in adult patients with sepsis (*69*), and two recent preclinical studies have demonstrated a protective effect of selective STAT3 inhibition in mice with CLP-induced sepsis (*70, 71*). Informed by these studies, we hypothesized that the degree and duration of STAT3 activation in Group C patients is pathologic and therefore a candidate target for precision immunomodulation. To test this hypothesis, we compared plasma protein expression between patients with MODS and patients with rare, monogenic inborn errors of immunity (IEI) involving constitutive activation or inactivation of the STAT1 and STAT3 signaling pathways. IEIs can provide context to help decode complex phenotypes by defining *in vivo* human immune signatures which correspond to signaling pathway perturbation. For this analysis, IEI patients were categorized as STAT1 gain-of-function (GOF) (n=9), STAT3 GOF (n=5), STAT1 autosomal dominant and dominant-negative (DN) (n=1), or STAT3 DN (n=3). Protein expression was measured using a 384-marker inflammatory proteomics panel from Olink Proteomics (**Table S7**), and we analyzed pathway expression in bulk by GSEA (*72*) and at the individual patient level using GSVA (*52*). **Fig. 6A** shows the study schematic for this analysis, which includes 88 patients with MODS, 18 patients with IEIs, and 25 HC participants. Orthogonal to our Ingenuity Pathway Analysis findings, we first confirmed population-level proteomic enrichment (normalized enrichment score 1.53, *p*-value 0.004) of the Hallmark IL-6/JAK/STAT3 signaling pathway in Group C patients at MODS onset compared to HC participants using GSEA (**Fig. 6B**). Leading edge proteins and other enriched Hallmark pathways from this analysis are shown in **Table S8**. To assess expression of STAT target proteins across MODS and IEI patients, we identified 31 measured proteins from the KEGG JAK/STAT gene set and performed hierarchical clustering based on normalized protein expression. As shown in **Fig. 6C**, we noted that STAT1 GOF and STAT3 GOF co-localize with Group C patients while STAT1 DN and STAT3 DN patients co-localize with Group A patients in hierarchical clustering. This protein-level analysis also highlights the residual biologic heterogeneity within our MODS severity-defined subphenotypes. Having demonstrated relative enrichment of STAT target proteins in MODS Group C patients, STAT1 GOF patients, and STAT3 GOF patients, we next used GSVA to assess patient-level enrichment of five inflammatory pathways in MODS and IEI patients. **Fig. 6D** shows a clustered heatmap of module enrichment scores for each patient in the dataset. While the cohort of STAT3 GOF and some STAT1 GOF patients show increased IL-6/JAK/STAT3 module enrichment scores, only STAT3 GOF patients display concurrent decreased PI3K/AKT/MTOR module scores. This pattern of concurrent increased STAT3 signaling and decreased mTOR activation was associated with PBMC immunometabolic dysregulation in our scRNA-seq analysis (**Fig. 5G**). Noting the similarities between Group C patients and STAT3 GOF patients, we next compared IL-6/JAK/STAT3 module enrichment scores across MODS and IEI patients (**Fig. 6E**). Compared to Group A patients, IL-6/JAK/STAT3 module enrichment scores were increased in Group C patients (*p*<0.001) and STAT3 GOF patients (*p*=0.001), and many patients with MODS had enrichment scores which exceeded the mean score of STAT3 GOF patients. Because patients with STAT3 GOF have constitutive pathway activation leading to immune dysregulation, this analysis suggests that Group C patients exhibit pathologic dysregulation of the STAT3 signaling pathway at MODS onset. Finally, we sought to determine the trajectory of STAT3 dysregulation in Group C patients. For this analysis, we computed longitudinal IL-6/JAK/STAT3 module enrichment scores among 26 patients with MODS in whom ≥3 longitudinal samples were obtained. We noted that Group C patients had the highest IL-6/JAK/STAT3 module enrichment scores, and that they remained persistently elevated through the first two weeks of MODS (**Fig. 6F**). In contrast, IL-6/JAK/STAT3 module enrichment scores gradually decreased in Group B patients and rapidly decreased in Group A patients. In this longitudinal analysis, it was uncommon for patient IL-6/JAK/STAT3 module enrichment scores to increase over time, with two Group C patients increasing from a low IL-6/JAK/STAT3 module enrichment scores to a higher score around one week after MODS onset. ![Fig. 6.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2024/06/11/2024.06.11.24308709/F6.medium.gif) [Fig. 6.](http://medrxiv.org/content/early/2024/06/11/2024.06.11.24308709/F6) Fig. 6. STAT3 signaling as a candidate target for precision immunomodulation in a subset of patients with MODS based on plasma proteomics. **(A)** Study schematic for this analysis, which includes 88 patients with MODS, 18 patients with IEIs, and 25 HC participants. Normalized protein expression was analyzed using GSEA and GSVA to determine similarities in protein expression between MODS subphenotypes and monogenic IEIs. **(B)** GSEA analysis of population-level enrichment of the Hallmark IL-6/JAK/STAT3 signaling pathway in Group C patients at MODS onset compared to HC participants. **(C)** Heatmap showing normalized protein expression of 31 KEGG JAK/STAT target proteins among 88 patients at MODS onset and 18 patients with inborn errors of immunity (IEI) characterized by STAT1 and STAT3 gain-of-function (GOF) and dominant-negative (DN) mutations. Colored annotation bars indicate MODS subphenotype and IEI subphenotype. **(D)** Clustered heatmap of patient-level enrichment of five inflammatory pathways in MODS patients and IEI patients using GSVA. Colored annotation bars indicate MODS subphenotype and IEI subphenotype. **(E)** Comparison of Hallmark IL-6/JAK/STAT3 signaling module enrichment scores by MODS subphenotype and IEI subphenotype. Comparisons were made between groups as shown and determined by Wilcoxon rank sum test. **(F)** Longitudinal IL-6/JAK/STAT3 module enrichment scores among 26 MODS patients with ≥3 sample timepoints. Solid lines represent individual patient trajectories; dashed lines represent loess regression lines for each MODS subphenotype. ## DISCUSSION In a prospective cohort of children admitted to a quaternary care PICU with MODS from diverse etiologies, we have identified an association between pathologic STAT3 hyperactivation and T cell immunometabolic dysregulation in a subset of pediatric patients with MODS with the highest severity of illness and poor prognosis. These patients can be identified by a 24-protein signature at MODS onset regardless of MODS etiology, and they continue to exhibit pathologic STAT3 hyperactivation for up to two weeks. Pathologic STAT3 signaling is targetable through approved anti-cytokine antibodies, JAK inhibitors, and experimental STAT inhibitors (*73*), and thus represents a viable target for precision medicine in pediatric patients with MODS. During the COVID-19 pandemic, some patients requiring supplemental oxygen were noted to benefit from a combination of dexamethasone and a JAK inhibitor (*74*). In randomized controlled trials, use of the JAK inhibitors baricitinib and tofacitinib were associated with reduced mortality (*75–78*), faster recovery (*79*), and lower rates of disease progression (*75*) in multiple adult trials focusing on varying levels of illness severity. Several available JAK inhibitor agents have good bioavailability, rapid onset of action, and short half-lives, features that are favorable for use in the critical care setting. If STAT3 hyperactivation is in fact a causal contributor to organ failure in pediatric MODS, the safety and efficacy of baricitinib in critical COVID-19 suggests a strong rationale for targeting JAK signaling in future pediatric clinical trials. Our highly granular approach yields insights complementary to many previous translational adult and pediatric sepsis studies and identifies a novel targetable mechanism of immune dysregulation in critically ill children with and without sepsis. As other investigators have shown, the highest severity patients in our cohort have innate (*17, 21*) and adaptive (*18, 20*) immune dysfunction and mitochondrial dysfunction (*25, 26, 80*). Our findings redemonstrate the association we previously reported (*23*) between T cell dysfunction, mitochondrial dysfunction, and clinical outcomes in pediatric sepsis, and suggest a potential causative molecular mechanism for sepsis-associated T cell immunometabolic dysregulation. Sepsis-3 criteria define sepsis as organ dysfunction that develops in the setting of a dysregulated response to infection (*1*), and nearly all patients with sepsis in this study had evidence of immune dysregulation. To our surprise, we identified features of immune dysregulation in many patients with MODS who did not have suspected infection, and our subphenotypes include patients with and without sepsis. These findings supplement the findings of previous studies (*60, 81*) which have focused on identifying clinical and transcriptional differences between patients with sepsis compared to sterile inflammation. Our findings demonstrating overlapping signatures of septic and sterile MODS may be explained by the high severity of illness in our cohort, as these patients may exhibit a final common pathway of critical illness which develops in the setting of overwhelming proinflammatory signals. Conversely, it is possible that populations of patients with MODS are enriched with children who have undiagnosed primary immune regulatory disorders. The PHENOMS study estimated that inborn errors of immunity could be present in up to 60% of pediatric patients with sepsis (*34*); increased prevalence of IEIs in pediatric MODS cohorts may therefore influence the association between immune dysregulation and clinical outcomes. Future multiomic immune profiling studies should also incorporate genomic data to identify novel genotype-phenotype associations which will further our mechanistic understanding of immune dysregulation in the setting of pediatric critical illness. Leveraging the proteomics dataset, we used penalized regression to identify a prognostic protein signature detectable at MODS onset and demonstrated excellent discrimination using linear discriminant analysis. Elastic net regularization is a common approach applied to feature selection in high-dimensional data which generates a high-performing sparse model with good prediction accuracy (*45*). Prior to potential use for prognostication or selection of targeted therapeutics, the proposed parsimonious 24-protein MODS severity model will require prospective validation. It will also benefit from validation across multiple proteomics platforms, as aptamer-based and antibody-based approaches are susceptible to both technical and genetic variation among samples (*82*). A validated severity model derived from the plasma proteome would present an opportunity for reverse translation to the bedside, as clinical labs in many centers can perform in-house multiplex protein assays. The longitudinal analysis of protein and pathway expression in this study provides evidence that immune dysregulation present at MODS onset persists for more than a week in many patients. Persistent immune dysregulation is amenable to consideration of a trial of precision therapeutics in which patients could be identified at MODS onset and then assigned to a treatment arm based on biomarker-based prognostic and predictive enrichment strategies. Equally relevant, prognostic subphenotypes defined at MODS onset demonstrated relatively stable proinflammatory pathway enrichment scores through time, suggesting that subphenotype membership is patient/episode-specific and minimally impacted by time from MODS onset. Cellular features of immune dysregulation associated with STAT3 hyperactivation were identified through both spectral flow cytometry and single cell transcriptomics and could serve as a surrogate marker for evaluating *ex vivo* response to candidate targeted therapeutics, including cytokine blockade and JAK inhibition. Assessment of PBMC mitochondrial function (e.g. via Seahorse or SCENITH (*83*) assays) could also serve as a candidate preclinical readout for response to immunomodulatory therapy. In our single cell analysis, IL-6/JAK/STAT3 module scores were moderately correlated with enrichment of sepsis response signature (SRS) genes in lymphoid subsets and not associated with SRS genes in myeloid subsets. This stands in contrast to a recent report of adult ICU patients with and without sepsis in which proportion of immature CD14+ MS1 monocytes was positively correlated with expression of SRS genes and TNF signaling (*84*). An alternative pathway of immune dysregulation has recently been described in which immature CD66b+ neutrophils drive emergency granulopoiesis via STAT3, which represents an immunocompromised disease endotype in adult patients with sepsis (*69*). Taken together, these observations suggest that immune dysregulation in pediatric MODS is somewhat distinct from that captured by the SRS gene signature, and that immature neutrophil driven emergency granulopoiesis may contribute to STAT3 hyperactivation and resulting immune dysregulation in our sickest pediatric patients with MODS. Further studies are required to test these hypotheses. Recent human and murine data identify a central role for STAT3 signaling in the pathogenesis of sepsis (*69–71*), and our results concordantly demonstrate associations between STAT3 signaling, clinical outcomes, and features of immune dysregulation implicated in pediatric sepsis pathobiology. However, the JAK/STAT signaling pathway is complex, and while we have focused on STAT3 as a candidate target for precision immunomodulation, it is likely that additional STAT pathways (including STAT5) are also dysregulated and may be targetable in pediatric patients with MODS. Future *ex vivo* studies assessing the impact of JAK inhibition on cellular immune phenotype and transcriptional profile are necessary to dissect the molecular mechanisms which underly these associations. Rare, monogenic IEIs offer unique insights into the mechanisms of human immune dysregulation and can provide context to help decode complex sepsis subphenotypes, and is a strength of our study. Typically identified in children with severe presentations of specific infectious diseases (*85*), IEIs are known to be enriched in pediatric patients with COVID-19 (*86–89*), influenza (*90–92*), and sepsis (*34, 93, 94*). Genotype-phenotype associations between monogenic IEIs and disease susceptibility may offer insights into causal pathobiology of illness which subsequently inform precision therapeutics (*95*). Our comparison between patients with MODS and patients with gain-of-function and dominant-negative STAT disorders allows us to learn from human disease the impact of chronic amplification of these key pathways and suggests that the extent of STAT3 hyperactivation in the most severe MODS subphenotype acutely exceeds that of patients with constitutive STAT3 activation, and that this hyperactivation is slow to resolve among survivors. Going forward, it will be important to replicate our prognostic MODS subphenotypes in a validation cohort using rapid, quantitative diagnostic test, as the highly multiplexed Olink assays are performed in batch in a reference lab. We must also understand the determinants of STAT3 hyperactivation in our defined subset of pediatric patients with MODS, which we hypothesize may reflect a host-specific susceptibility to excess inflammatory response in the setting of infectious and sterile insults. Integration of genomics, proteomics, and functional cellular assays into a single informatics pipeline may uncover new genotype-phenotype associations and highlight the key markers of immune dysregulation in patients with MODS. Limitations of our study include an *a priori* sample size that may have only been powered to resolve 3 subphenotypes of a heterogeneous syndrome. This reflects a pragmatic approach to study design and does not preclude the existence of other pediatric MODS subphenotypes. Additionally, our flow cytometry and transcriptomics experiments relied on cryopreserved PBMC samples. While this approach improved study feasibility and minimized batch effect associated with longitudinal sample collection, it is possible that rare and fragile cell populations may have been impacted by a single freeze/thaw cycle. The use of PBMCs also precludes cellular analysis of granulocyte populations, particularly neutrophils, which may be of particular relevance in the setting of STAT3 hyperactivation. Finally, we conducted our analysis on circulating blood cells and plasma protein expression based on the assumption that circulating immune cells reflect what happens in the injured tissues of patients with MODS, but the functional similarities between circulating immune cells and tissue resident immune cells are not known in the setting of hyperinflammation. Collectively, our longitudinal multiomics analysis defines a prognostic plasma proteomic signature present at MODS onset which is shared between patients with and without sepsis, demonstrates concordance across proteomic, cellular, and transcriptional analysis, and identifies a pathologic and potentially reversible role for STAT3 hyperactivation in a subset of pediatric patients with MODS who have high severity of illness and poor prognosis. These findings advance our understanding of MODS immunobiology and highlight potential opportunities for a precision medicine approach to the treatment of hyperinflammation in critically ill children. ## MATERIALS AND METHODS ### Study Design and Participants #### Experimental design After obtaining IRB approval (IRB #19-017032), we enrolled patients with multiple organ dysfunction syndrome (MODS) into a prospective observational cohort study in the Pediatric Intensive Care Unit and Cardiac Intensive Care Unit at Children’s Hospital of Philadelphia. Patients with new dysfunction of ≥2 organs defined by modified Proulx criteria (*96*) (**Table S9**) within the last three calendar days and age >40 weeks post-conceptual age and <18 years were eligible for enrollment. Exclusion criteria included limitations of care orders at the time of eligibility, clinical suspicion for brain death, and prior enrollment. This study co-enrolled with the multicenter PediAtric ReseArch of Drugs, Immunoparalysis and Genetics during MODS (PARADIGM) study, with shared inclusion/exclusion criteria and case report form but independent biospecimen collection, processing, and analysis. Patient (or legal guardian) consent (and assent, if appropriate) were obtained prior to study enrollment in accordance with our IRB approved protocol. We enrolled June 2020 to December 2022, when we reached our prespecified enrollment target of 88 patients, which we expected to achieve 90% power to detect moderate differences (Cohen’s d = 0.5) between three MODS subgroups, based on pilot flow cytometry and proteomics data. Longitudinal whole blood samples were obtained twice weekly in sodium heparin tubes from onset of MODS through death or resolution of all organ failure. #### Clinical metadata Clinical metadata were abstracted by a nurse research coordinator using a standardized case report form into REDCap (*97*) for the parent PARADIGM study. Data quality was verified through manual and automated queries. MODS onset was defined as the calendar day that a patient developed new dysfunction of ≥2 organs defined by modified Proulx criteria **(Table S9)** (*9*). MODS inciting diagnosis was abstracted from attending physician daily progress notes and coded as “sepsis,” “trauma,” “cardiopulmonary bypass,” or “non-infected.” Sepsis was identified as the MODS inciting diagnosis based on clinical suspicion for infection by the primary team and positive microbiologic or virologic testing. Patients with “culture negative sepsis” were classified as sepsis if they developed organ dysfunction in the setting of a clinical suspicion for infection and received sepsis therapies, in concordance with pediatric sepsis guidelines (*98, 99*). Immunocompromised status was defined by diagnosis of active malignancy, prior hematopoietic cell transplant, or primary immune deficiency syndrome. Daily data were extracted from MODS onset through day +28 to allow for daily calculation of the PELOD-2 organ dysfunction score (*13*). We defined cumulative PELOD-2 score through day +28 as our primary outcome because it incorporates both degree and duration of organ failure into a composite outcome variable which is associated with meaningful differences in long-term mortality (*100*) and health related quality of life (*101*) in pediatric patients with sepsis. For each patient and healthy control participant, a one-way hash function was used to generate a “hashed patient ID” number which cannot be linked to patient records or used to identify patients by anyone outside of the study team. These hashed patient IDs are included in supplemental tables where required. #### Healthy control participant samples Cryopreserved PBMC and heparin plasma samples from 25 pediatric participants without immunological disease were identified from an existing pediatric biorepository (IRB #18-015920) to serve as a healthy control group. Participant (or legal guardian) consent (and assent, if appropriate) were obtained prior to study enrollment in accordance with our IRB approved protocol. Age range and sex of healthy control participants are shown in **Table S10**. #### Inborn errors of immunity patient samples Heparin plasma samples from patients with STAT1 gain-of-function (9), STAT1 dominant-negative (1), STAT3 gain-of-function (5), and STAT3 dominant-negative (3) mutation were identified from collaborators and used as comparators for MODS patients with dysregulated STAT3 signaling. These participants were consented in accordance with local IRB protocols and samples were shared through collaborative research agreements. ### Study Procedures #### Sample processing Biospecimens were collected starting within 3 calendar days of MODS onset and continued twice weekly through death, recovery from all acute organ failures, or until six samples were obtained. Blood samples were obtained in heparin plasma tubes by clinical research coordinators and processed on site within one hour of acquisition. Platelet poor plasma was separated by centrifugation and snap frozen on dry ice. Peripheral blood mononuclear cells (PBMCs) were isolated using SepMate (STEMCELL Technologies, Vancouver, Canada) density gradient centrifugation using Lymphoprep media. Whole blood was mixed 1:1 with PBS and layered onto Lymphoprep gradient. SepMate tubes were centrifuged and the buffy coat suspension was spun down for cell isolation. ACK lysis was completed and cells were resuspended in complete RPMI (cRPMI; RPMI 1640 supplemented with 10% FBS, 1% L-Glutamine, 1% Pen-Strep) for counting prior to cryopreservation in 500μl of freezing media (90% FBS, 10% DMSO). #### Plasma proteomics by proximity extension assay For each patient and timepoint, 100μl of heparin plasma was submitted to Olink Proteomics (Uppsala, Sweden) for analysis on the Explore 1536 proximity extension assay (PEA) platform. PEA is a dual-recognition assay which uses paired antibodies labeled with DNA oligonucelotides to identify and quantify protein expression through Next Generation Sequencing (*42*). After normalization and quality control, protein concentrations are reported in log2-transformed Normalized Protein eXpression (NPX) units. In total, we analyzed 229 samples from 131 patients across three experiments using 8 shared bridging samples for between-plate-normalization. 1448 proteins met quality control thresholds in all three experiments and were included in downstream analyses (**Table S11**). #### Spectral flow cytometry staining and acquisition Cryopreserved PBMCs were thawed in 10mL cRPMI and 1x106 cells per sample were plated in a 96-well round-bottom plate. Cell pellets were sequentially incubated with live/dead blue with Fc block, surface antibody stain with Brilliant Stain buffer, permeabilization reagent, and intracellular antibody stain with Brilliant Stain buffer. After staining, cell pellets were resuspended in 1.6% PFA and held at 4°C until acquisition on a Cytek Aurora spectral flow cytometer (Cytek Biosciences, Fremont, CA) the following morning. In total, we analyzed 303 samples from 113 patients across eight experiments. Please see **Table S2** for details of flow cytometry antibodies and buffers and Supplemental Methods for complete staining protocol. #### Flow cytometric cell sorting for scRNA-seq Cryopreserved PBMCs were thawed in 10mL cRPMI and 1x106 cells were transferred to a new tube for staining. Cell pellets were incubated with 100μl mixture of live/dead blue and CD45 antibody, then washed and passed through a 35µm nylon mesh cell strainer prior to acquisition. Using a Cytek Aurora CS cell sorter (Cytek Biosciences, Fremont, CA), live singlet CD45+ cells were sorted into 1.5ml Eppendorf tubes prefilled with 500μl of PBS + 20% FBS. Sorted cells were washed with 500μl PBS + 10% FBS twice to remove sheath fluid EDTA prior to 10X processing. Please see Supplemental Methods for complete FACS staining protocol. #### Single cell RNA sequencing (scRNA-seq) We utilized the 10X Chromium Next GEM Single Cell 5’ Kit v2 and Chromium Single Cell Human TCR Amplification Kit (10X Genomics, Pleasanton, CA) for single cell analyses of sorted PBMC samples (see above) from 9 patients in MODS Group C and 3 pediatric HC participants. Single-cell isolation and library preparation were performed in the Center for Applied Genomics at Children’s Hospital of Philadelphia. Sequencing was performed using the Illumina S2 flow cell (Illumina, San Diego, CA). Data were demultiplexed and processed using the Cell Ranger 8.0 analytical pipeline (10X Genomics, Pleasanton, CA), with reads aligned to the GRCh38 reference genome. ### Statistical Analyses #### Identification of MODS subphenotypes We constructed linear mixed effects models for each protein to determine the association between normalized protein expression and cumulative PELOD-2 score, with age and sex modeled as fixed effects and day from MODS onset modeled as a random effect. Proteins which were significantly associated with cumulative PELOD-2 score after Benjamini-Hochberg correction (FDR *p*<0.05) were retained for further analysis. We used consensus clustering to define subphenotypes and identified optimal *k* via the Monte Carlo reference-based consensus clustering algorithm (*43*) which uses the proportion of ambiguous clustering (PAC) score to test the hypothesis that *k*=*n* clusters is more informative than *k*=1. Optimal number of clusters in the MODS data (k=3) was confirmed by Monte Carlo bootstrapping, elbow method, and gap statistic. The clinical and proteomic datasets contained no missing data elements, thus imputation was not necessary for our analyses. #### Competing risk survival analysis We fit a proportional subdistribution hazards regression model to assess the effect of covariates on the subdistribution of death and survival to PICU discharge in a competing risk setting and then estimated the cumulative incidence function from this model for each MODS subphenotype and outcome using the Fine and Gray model (*44*). We then tested the hypothesis that the subdistribution hazard differed between Group C and Group A/B for both death and survival to PICU discharge. #### Model reduction via ordinal elastic net We defined a parsimonious protein signature to classify MODS subphenotypes by training an ordinal elastic net model (*46*) for feature selection using the severity-associated proteins previously identified through linear mixed effects modeling. This approach fit a semi-parallel elementwise link multinomial-ordinal regression model with elastic net penalty which defines coefficients for each protein:subphenotype and then shrinks the model to a single coefficient for each protein via elastic net penalty to minimize Gini index. We assessed the performance of the parsimonious elastic net protein set to discriminate the three subphenotypes using principle component analysis, linear mixed effects post-estimation, and polytomous discrimination index (*102*). #### Identification of spectral flow cytometry metaclusters After arcsinh scaling (*47*) and quality control with flowAI (*48*), we performed FlowSOM metaclustering (*49*) to identify 14 PBMC populations by surface and intracellular marker expression. PBMCs were first clustered into k=60 FlowSOM metaclusters and then combined in a stepwise fashion to generate the 14 canonical populations identified in the figure and used for downstream analysis. Metacluster similarity was determined by surface and intracellular marker expression and tSNE-CUDA (*50*) proximity and then confirmed with manual gating. Proliferation and activation markers were analyzed using bivariate plots. #### Single cell transcriptomics analytic pipeline After sequencing and alignment to the GRCh38 reference genome, transcriptomics data were analyzed using the Seurat 5.0 preprocessing and integration pipeline (*103*). Briefly, after review of standard quality control metrics, we selected cells with 1200-20000 transcripts per cell, 750-5000 genes per cell, and <20% mitochondrial RNA content for downstream analysis. After standard normalization and scaling, we performed anchor-based RPCA integration across samples and conditions. After quality control and integration, cell identities were inferred using the Azimuth reference-based mapping pipeline (*55*) and used to define 14 PBMC populations by transcriptional profile. Results of each quality control and integration step are shown in **Fig. S4**. Pseudobulk analysis was performed by condition, sample, and cell type using Seurat’s AggregateExpression function. Differential expression analysis was performed using DESeq2 (*104*). The association of effector cell pathway enrichment and immunometabolic profile at the single cell level was determined using a generalized linear mixed-effects model incorporating pathway enrichment scores as fixed effects and patient identity as a random effect. #### Pathway enrichment analysis We first performed pathway enrichment analysis by subphenotype by analyzing estimated protein expression in Group C patients (adjusted for age, sex, and PELOD-2 score) with Ingenuity Pathway Analysis (Qiagen) (*51*). We then analyzed pathway expression in the protein dataset in bulk by GSEA (*72*) and at the individual patient level using GSVA (*52*), with a focus on enrichment of five canonical proinflammatory pathways using Human Molecular Signatures Database (MSigDB) Hallmark gene sets (*53*): > “TNF Signaling via NFkB” (HALLMARK\_TNFA\_SIGNALING\_VIA\_NFKB), “IL-6/JAK/STAT3 Signaling” (HALLMARK\_IL6\_JAK\_STAT3_SIGNALING), “IL-2/JAK/STAT5 Signaling” (HALLMARK_IL2_STAT5_SIGNALING), > > “IFN-γ Response” (HALLMARK_INTERFERON_GAMMA_RESPONSE), > > “PI3K/AKT/MTOR Signaling” (HALLMARK\_PI3K_AKT_MTOR_SIGNALING). We also applied GSVA to pseudobulk and single cell transcriptomics data in a similar approach using UCell (*58*) and extended our analysis to include two relevant Hallmark immunometabolic pathways: > “Glycolysis” (HALLMARK_GLYCOLYSIS), > > “Oxidative Phosphorylation” (HALLMARK_OXIDATIVE_PHOSPHORYLATION). A list of the gene sets used in GSEA and GSVA analysis in this manuscript are included in **Table S12**. Throughout the manuscript, we refer to Hallmark pathways as “modules” as opposed to “gene sets” because the analytic approach is applied to both proteomics data (via GSVA) and transcriptomics data (via UCell). #### SRS module concordance To determine the concordance between Hallmark proinflammatory module enrichment scores and the extended sepsis response signature (SRS) gene set, we calculated an extended SRS gene set enrichment score for each cell in the transcriptomics dataset using UCell (*58*). We then generated a linear regression model to compare to enrichment of the extended SRS gene set to each Hallmark proinflammatory pathway. We calculated the adjusted R2 for the correlation between extended SRS and inflammatory pathway enrichment by cell type. #### Analysis software All computational analysis was completed using R 4.3.2 and Bioconductor 3.18. Flow cytometry data was processed using FlowJo 10.9. Figures were compiled in Adobe Illustrator. ## Supporting information Supplementary Materials - Tables [[supplements/308709_file03.xlsx]](pending:yes) Supplementary Materials [[supplements/308709_file04.pdf]](pending:yes) ## Data Availability Code and de-identified study data will be made available at https://github.com/Lindell-Lab/MODS-Subphenotypes at the time of peer-reviewed publication. [https://github.com/Lindell-Lab/MODS-Subphenotypes](https://github.com/Lindell-Lab/MODS-Subphenotypes) ## Funding This study was funded in part by K12HD047349 (RBL), K08AI135091 (SEH), and pilot grants from the Thrasher Research Fund (RBL) and the CHOP Research Institute (RBL). The parent PARADIGM study is funded by R01HD095976 (MWH). Healthy control participant recruiting was funded by a Career Award for Medical Scientists from the Burroughs Wellcome Fund (SEH). Additional funding sources include the Barbara Brodsky Foundation and institutional startup funds from the Division of Critical Care Medicine. ## Author contributions Conceptualization: RBL, AFZ, EJW, NJM, SEH Methodology: RBL, JCF, NY, MK, SLW, MWH, AFZ, EJW, NJM, SEH Investigation: RBL, SS, JSC, AAM, CAH, PC, STF, TA, RT, AFF, JREB, SMH, JWL Data Curation: RBL, STF, TA Formal Analysis: RBL, MK, HF, DMT, RF, NJM, SEH Writing – Original Draft: RBL, NJM, SEH Writing – Reviewing & Editing: all authors Visualization: RBL Funding Acquisition: RBL, JCF, EMB, DTT, EJW, NJM, SEH Supervision: NJM, SEH ## Competing interests Children’s Hospital of Philadelphia and the University of Pennsylvania filed a provisional US patent application on May 22, 2024 related to the prognostic molecular subphenotypes presented in this manuscript (63/650,523; Identification of High Mortality Sub-Phenotypes of Pediatric Multiple Organ Dysfunction Syndrome (MODS) for Management and Treatment Thereof). RBL, NJM, and SEH are named as co-inventors on that patent application. All other authors declare that they have no competing interests. ## Data and materials availability Code and de-identified study data will be made available at [https://github.com/Lindell-Lab/MODS-Subphenotypes](https://github.com/Lindell-Lab/MODS-Subphenotypes) at the time of publication. MINSEQE-complaint transcriptional data will also be made available in the GEO repository at the time of publication. ## ACKNOWLEDGEMENTS We owe a debt of gratitude to all the patients and families who contributed clinical data and blood samples to this study, and to all involved clinical teams. We thank the many groups throughout CHOP and the University of Pennsylvania who contributed to this study, including the Pediatric Sepsis Program, the Dysregulated Immune Response Team, the Institute for Immunology & Immune Health, and the Department of Anesthesia and Critical Care. In particular, we thank the CHOP clinical research coordinator team (Nola Juste, Stephen Famularo, Teresa Arroyo, Barrington Bucknor, Ryan Burnett, Jenny Bush, Noah Buzinkai, Amanda Bwint, Amanda Cornetta, Shawn Dickey, Mary Ann DiLiberto, Rebecca Douglas, Emily Duffey, Glory Edioma, Emily Kirwin, Cindee Levow, Carly Machon, Alanah McKelvey, Erin Nfonoyim, Kathlyn Phengchomphet, Myooran Sivarupan, Taylor Slocumb, Bryn Spaide, Julia Surow, Samantha Switzer, Derrick Tam, Wandave Tizhe, Krithika Tupil, and Jia Yuan) for screening and enrolling eligible families in the study, the CHOP BioRC Specimen Processing Unit (Richard Tustin III, Annemarie Butler, Kate Kearns, Nicole Sundo, and Molly Tancini) for processing longitudinal human biospecimens, and the CHOP Center for Applied Genomics (Diana Slater, Cuiping Hou, and Mikaela Hufnell) for their assistance with single cell processing and sequencing. We are grateful to the PediAtric ReseArch of Drugs, Immunoparalysis and Genetics during MODS (PARADIGM) study team for allowing this study to proceed as a single-center ancillary to their multicenter study, and for generously sharing curated metadata from the parent study. Special thanks to members of the Collaborative Pediatric Critical Care Research Network (CPCCRN), Pediatric Acute Lung Injury and Sepsis Investigators (PALISI) Network, and Pediatric Critical Care and Trauma Scientist Development Program (PCCTSDP) for their advice during the development and conduct of this study. * Received June 11, 2024. * Revision received June 11, 2024. * Accepted June 11, 2024. * © 2024, Posted by Cold Spring Harbor Laboratory The copyright holder for this pre-print is the author. All rights reserved. The material may not be redistributed, re-used or adapted without the author's permission. ## REFERENCES AND NOTES 1. 1. M. Singer et al., The Third International Consensus Definitions for Sepsis and Septic Shock (Sepsis-3). JAMA 315, 801–810 (2016). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1001/jama.2016.0287&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=26903338&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F11%2F2024.06.11.24308709.atom) 2. 2. K. E. Rudd et al., Global, regional, and national sepsis incidence and mortality, 1990– 2017: analysis for the Global Burden of Disease Study. The Lancet 395, 200–211 (2020). 3. 3. S. L. Weiss et al., Global epidemiology of pediatric severe sepsis: the sepsis prevalence, outcomes, and therapies study. Am J Respir Crit Care Med 191, 1147–1157 (2015). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1164/rccm.201412-2323OC&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=25734408&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F11%2F2024.06.11.24308709.atom) 4. 4. M. S. Zinter, C. C. Dvorak, A. Spicer, M. J. Cowan, A. Sapru, New Insights Into Multicenter PICU Mortality Among Pediatric Hematopoietic Stem Cell Transplant Patients. Crit Care Med 43, 1986–1994 (2015). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1097/CCM.0000000000001085&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=26035280&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F11%2F2024.06.11.24308709.atom) 5. 5. R. B. Lindell et al., High Levels of Morbidity and Mortality Among Pediatric Hematopoietic Cell Transplant Recipients With Severe Sepsis: Insights From the Sepsis PRevalence, OUtcomes, and Therapies International Point Prevalence Study. Pediatr Crit Care Med 18, 1114–1125 (2017). 6. 6. R. B. Lindell, A. Nishisaki, S. L. Weiss, D. M. Traynor, J. C. Fitzgerald, Risk of Mortality in Immunocompromised Children With Severe Sepsis and Septic Shock. Crit Care Med 48, 1026–1033 (2020). 7. 7. S. L. Weiss et al., The Epidemiology of Hospital Death Following Pediatric Severe Sepsis: When, Why, and How Children With Sepsis Die. Pediatr Crit Care Med 18, 823–830 (2017). [PubMed](http://medrxiv.org/lookup/external-ref?access_num=http://www.n&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F11%2F2024.06.11.24308709.atom) 8. 8. J. C. Lin et al., New or Progressive Multiple Organ Dysfunction Syndrome in Pediatric Severe Sepsis: A Sepsis Phenotype With Higher Morbidity and Mortality. Pediatr Crit Care Med 18, 8–16 (2017). [PubMed](http://medrxiv.org/lookup/external-ref?access_num=28060151&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F11%2F2024.06.11.24308709.atom) 9. 9. F. Proulx et al., The pediatric multiple organ dysfunction syndrome. Pediatr Crit Care Med 10, 12–22 (2009). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1097/PCC.0b013e31819370a9&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=19057438&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F11%2F2024.06.11.24308709.atom) 10. 10. J. A. Carcillo et al., Pathophysiology of Pediatric Multiple Organ Dysfunction Syndrome. Pediatr Crit Care Med 18, S32–S45 (2017). [PubMed](http://medrxiv.org/lookup/external-ref?access_num=26652251&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F11%2F2024.06.11.24308709.atom) 11. 11. S. L. Weiss et al., Refining the Pediatric Multiple Organ Dysfunction Syndrome. Pediatrics 149, S13–S22 (2022). 12. 12. S. Leteurtre et al., Validation of the paediatric logistic organ dysfunction (PELOD) score: prospective, observational, multicentre study. Lancet 362, 192–197 (2003). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/S0140-6736(03)13908-6&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=12885479&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F11%2F2024.06.11.24308709.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000184232000008&link_type=ISI) 13. 13. S. Leteurtre et al., PELOD-2: an update of the PEdiatric logistic organ dysfunction score. Crit Care Med 41, 1761–1773 (2013). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1097/CCM.0b013e31828a2bbd&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=http://www.n&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F11%2F2024.06.11.24308709.atom) 14. 14. L. J. Schlapbach et al., Scoring Systems for Organ Dysfunction and Multiple Organ Dysfunction: The PODIUM Consensus Conference. Pediatrics 149, S23–S31 (2022). 15. 15. D. M. Maslove et al., Redefining critical illness. Nat Med 28, 1141–1148 (2022). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/s41591-022-01843-x&link_type=DOI) 16. 16. H. R. Wong, R. J. Freishtat, M. Monaco, K. Odoms, T. P. Shanley, Leukocyte subset-derived genomewide expression profiles in pediatric septic shock. Pediatr Crit Care Med 11, 349–355 (2010). [PubMed](http://medrxiv.org/lookup/external-ref?access_num=20009785&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F11%2F2024.06.11.24308709.atom) 17. 17. M. W. Hall et al., Immunoparalysis and nosocomial infection in children with multiple organ dysfunction syndrome. Intensive Care Med 37, 525–532 (2011). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1007/s00134-010-2088-x&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=21153402&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F11%2F2024.06.11.24308709.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000287512900021&link_type=ISI) 18. 18. J. A. Muszynski et al., Early adaptive immune suppression in children with septic shock: a prospective observational study. Crit Care 18, R145 (2014). 19. 19. J. A. Muszynski et al., Innate immune function predicts the development of nosocomial infection in critically injured children. Shock 42, 313–321 (2014). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1097/SHK.0000000000000217&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=24978895&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F11%2F2024.06.11.24308709.atom) 20. 20. J. A. Muszynski et al., Early Immune Function and Duration of Organ Dysfunction in Critically III Children with Sepsis. Am J Respir Crit Care Med 198, 361–369 (2018). 21. 21. M. W. Hall et al., Innate immune function and mortality in critically ill children with influenza: a multicenter study. Crit Care Med 41, 224–236 (2013). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1097/CCM.0b013e318267633c&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=23222256&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F11%2F2024.06.11.24308709.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000313150300025&link_type=ISI) 22. 22. R. S. Hotchkiss et al., Accelerated lymphocyte death in sepsis occurs by both the death receptor and mitochondrial pathways. J Immunol 174, 5110–5118 (2005). [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6ODoiamltbXVub2wiO3M6NToicmVzaWQiO3M6MTA6IjE3NC84LzUxMTAiO3M6NDoiYXRvbSI7czo1MDoiL21lZHJ4aXYvZWFybHkvMjAyNC8wNi8xMS8yMDI0LjA2LjExLjI0MzA4NzA5LmF0b20iO31zOjg6ImZyYWdtZW50IjtzOjA6IiI7fQ==) 23. 23. R. B. Lindell et al., Impaired Lymphocyte Responses in Pediatric Sepsis Vary by Pathogen Type and are Associated with Features of Immunometabolic Dysregulation. Shock 57, 191–199 (2022). 24. 24. S. L. Weiss et al., Mitochondrial dysfunction in peripheral blood mononuclear cells in pediatric septic shock. Pediatr Crit Care Med 16, e4–e12 (2015). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1097/PCC.0000000000000277&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=25251517&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F11%2F2024.06.11.24308709.atom) 25. 25. S. L. Weiss et al., Persistent Mitochondrial Dysfunction Linked to Prolonged Organ Dysfunction in Pediatric Sepsis. Crit Care Med, (2019). 26. 26. S. L. Weiss et al., Mitochondrial Dysfunction is Associated With an Immune Paralysis Phenotype in Pediatric Sepsis. Shock 54, 285–293 (2020). 27. 27. J. C. Marshall, Why have clinical trials in sepsis failed? Trends Mol Med 20, 195–203 (2014). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.molmed.2014.01.007&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=24581450&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F11%2F2024.06.11.24308709.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000335105700004&link_type=ISI) 28. 28. R. B. Lindell, N. J. Meyer, Interrogating the sepsis host immune response using cytomics. Crit Care 27, 93 (2023). 29. 29. R. B. Lindell, N. J. Meyer, Charting a course for precision therapy trials in sepsis. Lancet Respir Med 12, 265–267 (2024). 30. 30. H. R. Wong et al., Developing a clinically feasible personalized medicine approach to pediatric septic shock. Am J Respir Crit Care Med 191, 309–315 (2015). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1164/rccm.201410-1864OC&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=25489881&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F11%2F2024.06.11.24308709.atom) 31. 31. J. A. Carcillo et al., A Multicenter Network Assessment of Three Inflammation Phenotypes in Pediatric Sepsis-Induced Multiple Organ Failure. Pediatr Crit Care Med 20, 1137–1146 (2019). 32. 32. N. Yehya et al., Temperature Trajectory Sub-phenotypes and the Immuno-Inflammatory Response in Pediatric Sepsis. Shock 57, 645–651 (2022). 33. 33. L. N. Sanchez-Pinto et al., Patterns of Organ Dysfunction in Critically Ill Children Based on PODIUM Criteria. Pediatrics 149, S103–S110 (2022). 34. 34. K. F. Kernan et al., Prevalence of Pathogenic and Potentially Pathogenic Inborn Error of Immunity Associated Variants in Children with Severe Sepsis. J Clin Immunol 42, 350–364 (2022). 35. 35. B. P. Scicluna et al., Classification of patients with sepsis according to blood genomic endotype: a prospective cohort study. Lancet Respir Med 5, 816–826 (2017). 36. 36. E. E. Davenport et al., Genomic landscape of the individual host response and outcomes in sepsis: a prospective cohort study. Lancet Respir Med 4, 259–271 (2016). 37. 37. T. E. Sweeney et al., Unsupervised Analysis of Transcriptomics in Bacterial Sepsis Across Multiple Datasets Reveals Three Robust Clusters. Crit Care Med 46, 915–925 (2018). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1097/CCM.0000000000003084&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=29537985&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F11%2F2024.06.11.24308709.atom) 38. 38. C. W. Seymour et al., Derivation, Validation, and Potential Treatment Implications of Novel Clinical Phenotypes for Sepsis. JAMA 321, 2003–2017 (2019). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1001/jama.2019.5791&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=http://www.n&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F11%2F2024.06.11.24308709.atom) 39. 39. H. R. Wong et al., Prospective clinical testing and experimental validation of the Pediatric Sepsis Biomarker Risk Model. Sci Transl Med 11, (2019). 40. 40. E. Cano-Gamez et al., An immune dysfunction score for stratification of patients with acute infection based on whole-blood gene expression. Sci Transl Med 14, eabq4433 (2022). 41. 41. P. Sinha et al., Identifying molecular phenotypes in sepsis: an analysis of two prospective observational cohorts and secondary analysis of two randomised controlled trials. Lancet Respir Med 11, 965–974 (2023). 42. 42. L. Wik et al., Proximity Extension Assay in Combination with Next-Generation Sequencing for High-throughput Proteome-wide Analysis. Mol Cell Proteomics 20, 100168 (2021). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.mcpro.2021.100168&link_type=DOI) 43. 43. C. R. John et al., M3C: Monte Carlo reference-based consensus clustering. Sci Rep 10, 1816 (2020). 44. 44. G. R. Fine JP, A Proportional Hazards Model for the Subdistribution of a Competing Risk. Journal of the American Statistical Association 94, 496–509 (1999). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.2307/2670170&link_type=DOI) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000081058500019&link_type=ISI) 45. 45. H. Zou, T. Hastie, Regularization and Variable Selection Via the Elastic Net. Journal of the Royal Statistical Society Series B: Statistical Methodology 67, 301–320 (2005). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1111/j.1467-9868.2005.00503.x&link_type=DOI) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000227498200007&link_type=ISI) 46. 46. M. J. Wurm, P. J. Rathouz, B. M. Hanlon, Regularized Ordinal Regression and the ordinalNet R Package. J Stat Softw 99, (2021). 47. 47. H. Chen et al., Cytofkit: A Bioconductor Package for an Integrated Mass Cytometry Data Analysis Pipeline. PLoS Comput Biol 12, e1005112 (2016). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1371/journal.pcbi.1005112&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=27662185&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F11%2F2024.06.11.24308709.atom) 48. 48. G. Monaco et al., flowAI: automatic and interactive anomaly discerning tools for flow cytometry data. Bioinformatics 32, 2473–2480 (2016). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1093/bioinformatics/btw191&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=27153628&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F11%2F2024.06.11.24308709.atom) 49. 49. S. Van Gassen et al., FlowSOM: Using self-organizing maps for visualization and interpretation of cytometry data. Cytometry A 87, 636–645 (2015). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1002/cyto.a.22625&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=25573116&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F11%2F2024.06.11.24308709.atom) 50. 50. R. R. D. M. Chan, F. Huang and J. F. Canny,, paper presented at the 2018 30th International Symposium on Computer Architecture and High Performance Computing, Lyon, France, 2018. 51. 51. A. Kramer, J. Green, J. Pollard, Jr., S. Tugendreich, Causal analysis approaches in Ingenuity Pathway Analysis. Bioinformatics 30, 523–530 (2014). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1093/bioinformatics/btt703&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=24336805&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F11%2F2024.06.11.24308709.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000332032100010&link_type=ISI) 52. 52. S. Hanzelmann, R. Castelo, J. Guinney, GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics 14, 7 (2013). 53. 53. A. Liberzon et al., The Molecular Signatures Database (MSigDB) hallmark gene set collection. Cell Syst 1, 417–425 (2015). 54. 54. D. C. Fajgenbaum, C. H. June, Cytokine Storm. N Engl J Med 383, 2255–2273 (2020). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1056/NEJMra2026131&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=http://www.n&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F11%2F2024.06.11.24308709.atom) 55. 55. Y. Hao et al., Integrated analysis of multimodal single-cell data. Cell 184, 3573–3587 e3529 (2021). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/J.CELL.2021.04.048&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=http://www.n&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F11%2F2024.06.11.24308709.atom) 56. 56. L. McInnes, J. Healy, N. Saul, L. Grobberger, UMAP: Uniform Manifold Approximation and Projection. Journal of Open Source Software 3, 861 (2018). 57. 57. Y. Kashima et al., Potentiality of multiple modalities for single-cell analyses to evaluate the tumor microenvironment in clinical specimens. Sci Rep 11, 341 (2021). 58. 58. M. Andreatta, S. J. Carmona, UCell: Robust and scalable single-cell gene signature scoring. Comput Struct Biotechnol J 19, 3796–3798 (2021). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.csbj.2021.06.043&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=34285779&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F11%2F2024.06.11.24308709.atom) 59. 59. J. W. Squair et al., Confronting false discoveries in single-cell differential expression. Nat Commun 12, 5692 (2021). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/s41467-021-25960-2&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=34584091&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F11%2F2024.06.11.24308709.atom) 60. 60. T. E. Sweeney, A. Shidham, H. R. Wong, P. Khatri, A comprehensive time-course-based multicohort analysis of sepsis and sterile inflammation reveals a robust diagnostic gene set. Sci Transl Med 7, 287ra271 (2015). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1126/scitranslmed.aaa8038&link_type=DOI) 61. 61. T. E. Sweeney et al., A community approach to mortality prediction in sepsis via gene expression analysis. Nat Commun 9, 694 (2018). 62. 62. T. Wang et al., Regulation of the innate and adaptive immune responses by Stat-3 signaling in tumor cells. Nat Med 10, 48–54 (2004). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/nm976&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=14702634&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F11%2F2024.06.11.24308709.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000187743600038&link_type=ISI) 63. 63. K. Takeda et al., Enhanced Th1 activity and development of chronic enterocolitis in mice devoid of Stat3 in macrophages and neutrophils. Immunity 10, 39–49 (1999). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/S1074-7613(00)80005-9&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=10023769&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F11%2F2024.06.11.24308709.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000078428700005&link_type=ISI) 64. 64. Y. Nefedova et al., Regulation of dendritic cell differentiation and antitumor immune response in cancer by pharmacologic-selective inhibition of the janus-activated kinase 2/signal transducers and activators of transcription 3 pathway. Cancer Res 65, 9525–9535 (2005). [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6NjoiY2FucmVzIjtzOjU6InJlc2lkIjtzOjEwOiI2NS8yMC85NTI1IjtzOjQ6ImF0b20iO3M6NTA6Ii9tZWRyeGl2L2Vhcmx5LzIwMjQvMDYvMTEvMjAyNC4wNi4xMS4yNDMwODcwOS5hdG9tIjt9czo4OiJmcmFnbWVudCI7czowOiIiO30=) 65. 65. E. Zorn et al., IL-2 regulates FOXP3 expression in human CD4+CD25+ regulatory T cells through a STAT-dependent mechanism and induces the expansion of these cells in vivo. Blood 108, 1571–1579 (2006). [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6MTI6ImJsb29kam91cm5hbCI7czo1OiJyZXNpZCI7czoxMDoiMTA4LzUvMTU3MSI7czo0OiJhdG9tIjtzOjUwOiIvbWVkcnhpdi9lYXJseS8yMDI0LzA2LzExLzIwMjQuMDYuMTEuMjQzMDg3MDkuYXRvbSI7fXM6ODoiZnJhZ21lbnQiO3M6MDoiIjt9) 66. 66. J. Wegrzyn et al., Function of mitochondrial Stat3 in cellular respiration. Science 323, 793–797 (2009). [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6Mzoic2NpIjtzOjU6InJlc2lkIjtzOjEyOiIzMjMvNTkxNS83OTMiO3M6NDoiYXRvbSI7czo1MDoiL21lZHJ4aXYvZWFybHkvMjAyNC8wNi8xMS8yMDI0LjA2LjExLjI0MzA4NzA5LmF0b20iO31zOjg6ImZyYWdtZW50IjtzOjA6IiI7fQ==) 67. 67. J. J. Balic et al., STAT3 serine phosphorylation is required for TLR4 metabolic reprogramming and IL-1beta expression. Nat Commun 11, 3816 (2020). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/s41467-020-17669-5&link_type=DOI) 68. 68. R. Clere-Jehl et al., JAK-STAT Targeting Offers Novel Therapeutic Opportunities in Sepsis. Trends Mol Med 26, 987–1002 (2020). 69. 69. A. J. Kwok et al., Neutrophils and emergency granulopoiesis drive immune suppression and an extreme response endotype during sepsis. Nat Immunol 24, 767–779 (2023). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/s41590-023-01490-5&link_type=DOI) 70. 70. S. Xu et al., Phospho-Tyr705 of STAT3 is a therapeutic target for sepsis through regulating inflammation and coagulation. Cell Commun Signal 18, 104 (2020). 71. 71. S. Imbaby et al., Beneficial effect of STAT3 decoy oligodeoxynucleotide transfection on organ injury and mortality in mice with cecal ligation and puncture-induced sepsis. Sci Rep 10, 15316 (2020). 72. 72. A. Subramanian et al., Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A 102, 15545–15550 (2005). [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6NDoicG5hcyI7czo1OiJyZXNpZCI7czoxMjoiMTAyLzQzLzE1NTQ1IjtzOjQ6ImF0b20iO3M6NTA6Ii9tZWRyeGl2L2Vhcmx5LzIwMjQvMDYvMTEvMjAyNC4wNi4xMS4yNDMwODcwOS5hdG9tIjt9czo4OiJmcmFnbWVudCI7czowOiIiO30=) 73. 73. X. Hu, J. Li, M. Fu, X. Zhao, W. Wang, The JAK/STAT signaling pathway: from bench to clinic. Signal Transduct Target Ther 6, 402 (2021). 74. 74.74. COVID-19 Treatment Guidelines Panel, Coronavirus Disease 2019 (COVID-19) Treatment Guidelines. National Institutes of Health. *Available at* [https://www.covid19treatmentguidelines.nih.gov/](https://www.covid19treatmentguidelines.nih.gov/). Accessed April 1, 2024. 75. 75. P. O. Guimaraes et al., Tofacitinib in Patients Hospitalized with Covid-19 Pneumonia. N Engl J Med 385, 406–415 (2021). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1056/nejmoa2101643&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=http://www.n&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F11%2F2024.06.11.24308709.atom) 76. 76. V. C. Marconi et al., Efficacy and safety of baricitinib for the treatment of hospitalised adults with COVID-19 (COV-BARRIER): a randomised, double-blind, parallel-group, placebo-controlled phase 3 trial. Lancet Respir Med 9, 1407–1418 (2021). 77. 77. E. W. Ely et al., Efficacy and safety of baricitinib plus standard of care for the treatment of critically ill hospitalised adults with COVID-19 on invasive mechanical ventilation or extracorporeal membrane oxygenation: an exploratory, randomised, placebo-controlled trial. Lancet Respir Med 10, 327–336 (2022). 78. 78. R. C. Group, Baricitinib in patients admitted to hospital with COVID-19 (RECOVERY): a randomised, controlled, open-label, platform trial and updated meta-analysis. Lancet 400, 359–368 (2022). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/S0140-6736(22)01109-6&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=http://www.n&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F11%2F2024.06.11.24308709.atom) 79. 79. A. C. Kalil et al., Baricitinib plus Remdesivir for Hospitalized Adults with Covid-19. N Engl J Med 384, 795–807 (2021). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1056/nejmoa2031994&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=http://www.n&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F11%2F2024.06.11.24308709.atom) 80. 80. R. B. Lindell et al., Impaired Lymphocyte Responses in Pediatric Sepsis Vary by Pathogen Type. medRxiv, (2021). 81. 81. S. Thair et al., Gene Expression-Based Diagnosis of Infections in Critically Ill Patients-Prospective Validation of the SepsisMetaScore in a Longitudinal Severe Trauma Cohort. Crit Care Med 49, e751–e760 (2021). 82. 82. M. Pietzner et al., Synergistic insights into human health from aptamer- and antibody-based proteomic profiling. Nat Commun 12, 6822 (2021). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/s41467-021-27164-0&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=34819519&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F11%2F2024.06.11.24308709.atom) 83. 83. R. J. Arguello et al., SCENITH: A Flow Cytometry-Based Method to Functionally Profile Energy Metabolism with Single-Cell Resolution. Cell Metab 32, 1063–1075 e1067 (2020). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.cmet.2020.11.007&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=33264598&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F11%2F2024.06.11.24308709.atom) 84. 84. G. G. F. Leite et al., Monocyte state 1 (MS1) cells in critically ill patients with sepsis or non-infectious conditions: association with disease course and host response. Crit Care 28, 88 (2024). 85. 85. J. L. Casanova, Severe infectious diseases of childhood as monogenic inborn errors of immunity. Proc Natl Acad Sci U S A 112, E7128–7137 (2015). [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6NDoicG5hcyI7czo1OiJyZXNpZCI7czoxMjoiMTEyLzUxL0U3MTI4IjtzOjQ6ImF0b20iO3M6NTA6Ii9tZWRyeGl2L2Vhcmx5LzIwMjQvMDYvMTEvMjAyNC4wNi4xMS4yNDMwODcwOS5hdG9tIjt9czo4OiJmcmFnbWVudCI7czowOiIiO30=) 86. 86. C. I. van der Made et al., Presence of Genetic Variants Among Young Men With Severe COVID-19. JAMA 324, 663–673 (2020). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1001/jama.2020.13719&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=http://www.n&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F11%2F2024.06.11.24308709.atom) 87. 87. P. Y. Lee et al., Immune dysregulation and multisystem inflammatory syndrome in children (MIS-C) in individuals with haploinsufficiency of SOCS1. J Allergy Clin Immunol 146, 1194–1200 e1191 (2020). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.jaci.2020.07.033&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=32853638&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F11%2F2024.06.11.24308709.atom) 88. 88. H. Abolhassani et al., Inherited IFNAR1 Deficiency in a Child with Both Critical COVID-19 Pneumonia and Multisystem Inflammatory Syndrome. J Clin Immunol 42, 471–483 (2022). 89. 89. L. Schidlowski, A. P. D. Iwamura, S. U. D. Covid, A. Condino-Neto, C. Prando, Diagnosis of APS-1 in Two Siblings Following Life-Threatening COVID-19 Pneumonia. J Clin Immunol 42, 749–752 (2022). 90. 90. M. J. Ciancanelli et al., Infectious disease. Life-threatening influenza and impaired interferon amplification in human IRF7 deficiency. Science 348, 448–453 (2015). [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6Mzoic2NpIjtzOjU6InJlc2lkIjtzOjEyOiIzNDgvNjIzMy80NDgiO3M6NDoiYXRvbSI7czo1MDoiL21lZHJ4aXYvZWFybHkvMjAyNC8wNi8xMS8yMDI0LjA2LjExLjI0MzA4NzA5LmF0b20iO31zOjg6ImZyYWdtZW50IjtzOjA6IiI7fQ==) 91. 91. N. Hernandez et al., Life-threatening influenza pneumonitis in a child with inherited IRF9 deficiency. J Exp Med 215, 2567–2585 (2018). [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6MzoiamVtIjtzOjU6InJlc2lkIjtzOjExOiIyMTUvMTAvMjU2NyI7czo0OiJhdG9tIjtzOjUwOiIvbWVkcnhpdi9lYXJseS8yMDI0LzA2LzExLzIwMjQuMDYuMTEuMjQzMDg3MDkuYXRvbSI7fXM6ODoiZnJhZ21lbnQiO3M6MDoiIjt9) 92. 92. H. K. Lim et al., Severe influenza pneumonitis in children with inherited TLR3 deficiency. J Exp Med 216, 2038–2056 (2019). [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6MzoiamVtIjtzOjU6InJlc2lkIjtzOjEwOiIyMTYvOS8yMDM4IjtzOjQ6ImF0b20iO3M6NTA6Ii9tZWRyeGl2L2Vhcmx5LzIwMjQvMDYvMTEvMjAyNC4wNi4xMS4yNDMwODcwOS5hdG9tIjt9czo4OiJmcmFnbWVudCI7czowOiIiO30=) 93. 93. S. Asgari et al., Exome Sequencing Reveals Primary Immunodeficiencies in Children with Community-Acquired Pseudomonas aeruginosa Sepsis. Front Immunol 7, 357 (2016). 94. 94. A. Borghesi et al., Whole-exome Sequencing for the Identification of Rare Variants in Primary Immunodeficiency Genes in Children With Sepsis: A Prospective, Population-based Cohort Study. Clin Infect Dis 71, e614–e623 (2020). 95. 95. A. J. Kwok, A. Mentzer, J. C. Knight, Host genetics and infectious disease: new tools, insights and translational opportunities. Nat Rev Genet 22, 137–153 (2021). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/s41576-020-00297-6&link_type=DOI) 96. 96. F. Baudin et al., Modalities and Complications Associated With the Use of High-Flow Nasal Cannula: Experience in a Pediatric ICU. Respir Care 61, 1305–1310 (2016). [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6ODoicmVzcGNhcmUiO3M6NToicmVzaWQiO3M6MTA6IjYxLzEwLzEzMDUiO3M6NDoiYXRvbSI7czo1MDoiL21lZHJ4aXYvZWFybHkvMjAyNC8wNi8xMS8yMDI0LjA2LjExLjI0MzA4NzA5LmF0b20iO31zOjg6ImZyYWdtZW50IjtzOjA6IiI7fQ==) 97. 97. P. A. Harris et al., Research electronic data capture (REDCap)--a metadata-driven methodology and workflow process for providing translational research informatics support. J Biomed Inform 42, 377–381 (2009). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.jbi.2008.08.010&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=18929686&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F11%2F2024.06.11.24308709.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000264958800018&link_type=ISI) 98. 98. B. Goldstein, B. Giroir, A. Randolph, S. International Consensus Conference on Pediatric, International pediatric sepsis consensus conference: definitions for sepsis and organ dysfunction in pediatrics. Pediatr Crit Care Med 6, 2–8 (2005). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1097/01.PCC.0000149131.72248.E6&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=15636651&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F11%2F2024.06.11.24308709.atom) 99. 99. L. J. Schlapbach et al., International Consensus Criteria for Pediatric Sepsis and Septic Shock. JAMA 331, 665–674 (2024). 100.100. J. J. Zimmerman et al., Critical Illness Factors Associated With Long-Term Mortality and Health-Related Quality of Life Morbidity Following Community-Acquired Pediatric Septic Shock. Crit Care Med 48, 319–328 (2020). 101.101. K. L. Meert et al., Trajectories and Risk Factors for Altered Physical and Psychosocial Health-Related Quality of Life After Pediatric Community-Acquired Septic Shock. Pediatr Crit Care Med 21, 869–878 (2020). 102.102. B. Van Calster et al., Extending the c-statistic to nominal polytomous outcomes: the Polytomous Discrimination Index. Stat Med 31, 2610–2626 (2012). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1002/sim.5321&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=22733650&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2024%2F06%2F11%2F2024.06.11.24308709.atom) 103.103. Y. Hao et al., Dictionary learning for integrative, multimodal and scalable single-cell analysis. Nat Biotechnol 42, 293–304 (2024). 104.104. M. I. Love, W. Huber, S. Anders, Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 15, 550 (2014).