Transcriptomic Profiles in Children With Septic Shock With or Without Immunoparalysis

Background Severe innate immune suppression, termed immunoparalysis, is associated with increased risks of nosocomial infection and mortality in children with septic shock. Currently, immunoparalysis cannot be clinically diagnosed in children, and mechanisms remain unclear. Transcriptomic studies identify subsets of septic children with downregulation of genes within adaptive immune pathways, but assays of immune function have not been performed as part of these studies, and little is known about transcriptomic profiles of children with immunoparalysis. Methods We performed a nested case-control study to identify differences in RNA expression patterns between children with septic shock with immunoparalysis (defined as lipopolysaccharide (LPS)-induced tumor necrosis factor (TNF)α response < 200 pg/ml) vs those with normal LPS-induced TNFα response. Children were enrolled within 48 hours of the onset of septic shock and divided into two groups based on LPS-induced TNFα response. RNA was extracted from whole blood for RNAseq, differential expression analyses using DESeq2 software, and pathway analyses using Ingenuity Pathway Analysis. Results 32 children were included in analyses. Comparing those with immunoparalysis (n =19) to those with normal TNFα response (n = 13), 2,303 transcripts were differentially expressed with absolute value fold change ≥ 1.5 and false discovery rate ≤ 0.05. The majority of downregulated pathways in children with immunoparalysis were pathways that involved interactions between innate and adaptive immune cells necessary for cell-mediated immunity, crosstalk between dendritic cells and natural killer cells, and natural killer cell signaling pathways. Upregulated pathways included those involved in humoral immunity (T helper cell type 2), corticotropin signaling, platelet activation (GP6 signaling), and leukocyte migration and extravasation. Conclusions Our study suggests that gene expression data might be useful to identify children with immunoparalysis and identifies several key differentially regulated pathways involved in both innate and adaptive immunity. Our ongoing work in this area aims to dissect interactions between innate and adaptive immunity in septic children and to more fully elucidate patient-specific immunologic pathophysiology to guide individualized immunotherapeutic targets.


INTRODUCTION
Septic shock is an important source of morbidity and mortality in children worldwide (1,2). While the signs and symptoms of sepsis arise from an overwhelming inflammatory response, work from our group and others has demonstrated significant associations between suppressed innate immune response and adverse outcomes in critically ill septic children (3)(4)(5). Sepsisinduced innate immune suppression is typically occult but can be identified by quantifying production of tumor necrosis factor (TNF)a, an inflammatory cytokine, in whole blood after stimulation with lipopolysaccharide (LPS). Using this standardized assay, we have previously identified a threshold ex vivo LPS-induced TNFa response of < 200 pg/ml, termed immunoparalysis, that is associated with the development of nosocomial infection, prolonged organ dysfunction, and/or increased mortality in multiple studies of children with severe sepsis and/or septic shock (3)(4)(5). While the ex vivo LPS-induced TNFa assay is highly reproducible, its widespread use is limited by the need for a centralized laboratory for reagent preparation and cytokine quantitation, and on-site sample processing at the time of the blood draw. Alternate methods to identify children with clinically relevant ICU-related immune suppression are needed to facilitate multicenter clinical trial design and translation to the bedside.
Previous transcriptomic studies of children with septic shock have identified associations between downregulation of genes involved in adaptive immunity and adverse clinical outcomes, suggesting that gene expression patterns might be useful to identify children with septic shock with clinically relevant immune suppression (6)(7)(8). However, assays of immune cell function were not performed in these studies, and little is known about gene expression profiles in children with sepsis-induced immunoparalysis. We therefore performed a pilot nested casecontrol study to identify differences in RNA expression patterns between children with septic shock with immunoparalysis versus those with normal LPS-induced TNFa response.

Setting
The pediatric intensive care unit (PICU) at Nationwide Children's Hospital is a 54-bed, quaternary care, medical/ surgical intensive care unit that sees approximately 3,000 annual admissions. Patients admitted to a separate cardiothoracic intensive care unit were not included in this study. The study was approved by the Institutional Review Board at Nationwide Children's Hospital.

Patients
After informed consent, children who were within 48 hours of septic shock onset were enrolled. Septic shock was defined by the presence of sepsis (per consensus criteria) and the use of continuous vasoactive medication for shock (9). Children were excluded from the parent study if they had a limitation of care order in place or likely progression to brain death (as determined by the treating attending physician). Children who required extracorporeal membrane oxygenation support within 48 hours of septic shock onset were excluded from these analyses because of the likelihood that their transcriptomic profiles would be very different and an insufficient number of patients to permit stratified analyses. To ensure that enough RNA would be available for analysis, children were also excluded if their total white blood cell count was < 1,000 cells/mm 3 .
LPS-induced TNFa production capacity was assessed within 48 hours of septic shock onset as previously described (5). Briefly, patient samples were drawn into heparinized blood tubes (Becton Dickson, Franklin Lakes, NJ) within 48 hours of the onset of septic shock. 50 microliters of whole blood were added to pre-prepared stimulation solution containing RPMI and 0.5 ng/ml LPS (phenolextracted from Salmonella abortus equi [Sigma, St. Louis, MO]) and incubated for four hours at 37°C. After four hours, the supernatant was collected and stored at -80°C for batch analysis of TNFa, which was quantified by chemiluminescence using the Immulite 1000 automated chemiluminometer (Siemens Healthcare Diagnostics, Deerfield, Il). Stimulation assays were performed in duplicate and values reported are the average values from each set of duplicates.
To ensure adequate separation between groups, children were selected and divided into two groups based on their LPS-induced TNFa response as follows: Immunoparalysis (LPS-induced TNFa response < 200 pg/ml) and normal innate immune function (LPS-induced TNFa response > 830 pg/ml). The definition of immunoparalysis is derived from our previous publications (3-5, 10, 11). The threshold of 830 pg/ml represents the 25%ile of LPS-induced TNFa response in healthy children (5). Children with an initial LPS-induced TNFa response between 200 and 830 pg/ml were not included in these analyses. As we were interested in stable immune phenotypes, children whose LPS-induced TNFa response significantly changed over the course of 48 hours (i.e. transient immunoparalysis or transient normal) were also not included.

Clinical Data Collection and Definitions
Clinical data was abstracted from the electronic medical record into a standardized, electronic case report form by trained research coordinators. Complex chronic conditions were defined as previously published (12). Severity of illness was quantified by the pediatric risk of mortality (PRISMIII) score on the day of septic shock onset (13). Vasoactive inotrope score was used to define intensity of treatment for shock, as previously published (14). Nosocomial sepsis was defined as the onset of septic shock after 48 hours of hospital admission.

RNA Extraction
Blood samples for RNA extraction were drawn at the same time as samples drawn for LPS-induced TNFa production capacity assays. Whole blood (0.6 to 1ml) was drawn into tubes containing 1 ml of Tempus ™ RNA stabilization solution (Applied Biosystems). The tubes were shaken vigorously then immediately frozen and stored at -20°C for batch analysis. RNA was isolated using the MagMax for Stabilized Blood Tubes RNA Isolation kit (Life Technologies) via a modification of the manufacture's protocol to permit small volume blood sampling as follows (15). 100 ul of Tempus 1x PBS was added to each blood sample and vortexed for 30 seconds. Tubes were centrifuged at 4,000 x g for 20 minutes at 4°C to pellet the crude RNA. Pellets were washed with 4ml Tempus Pre-Digestion Wash, centrifuged at 4,000 x g, resuspended in 120 ul of Tempus Resuspension Solution and Proteinase and vortexed for 10 seconds to resuspend the pellet. 10 ul of TURBO ™ DNase was added to the tubes and vortexed for 10 minutes. 50 ul of Binding Solution Concentrate and 20 ul of RNA binding beads were then added and the tube was vortexed for 1 minute. 200 ul of 100% isopropanol was added for another 3-minute vortex. A magnetic stand was used to capture the RNA binding beads with the supernatant discarded. 150 ul of Wash Solution 1 was added to the tube and vortexed for 1 minute. The tube was then placed on the magnetic stand to capture the RNA Binding beads, discarding the supernatant. This wash was repeated once. 150 ul of Wash Solution 2 was added to the tube and vortexed for 1 minute then put back on the magnetic stand. The supernatant was carefully discarded, and the wash was repeated. The tubes were inverted on absorbent paper for 2 minutes to dry the beads then 60 ul of Elution Buffer was added. The tubes were vortexed for 4 minutes and put back on the magnetic stand. The RNAcontaining supernatant was extracted from the tubes and stored at -20°C.

RNA Sequencing
Strand-specific RNA-seq libraries were prepared using NEBNext Ultra II Directional RNA Library Prep Kit for Illumina, following the manufacturer's recommendations. In summary, total RNA quality was assessed using RNA 6000 Nano kit on Agilent 2100 Bioanalyzer (Agilent Biotechnologies) and concentration measured using Qubit RNA HS assay kit (Life Technologies). An aliquot of total RNA was rRNA depleted using NEB's Human/Mouse/Rat RNAse-H based Depletion kit (New England BioLabs). Following rRNA removal, mRNA was fragmented and then used for firstand second-strand cDNA synthesis with random hexamer primers. ds cDNA fragments underwent end-repair and a-tailing and ligated to dual-unique adapters (Integrated DNA Technologies). Adaptor-ligated cDNA was amplified by limit-cycle PCR. Library quality was analyzed on Tapestation High-Sensitivity D1000 ScreenTape (Agilent Biotechnologies) and quantified by KAPA qPCR (KAPA BioSystems). Libraries were pooled and sequenced at 2 x 150 bp read lengths to generate approximately 60-80 million paired-end reads per sample.

Statistical Analyses
Clinical characteristics were compared using Fisher's exact test (categorical variables) or Wilcoxon rank sum test (continuous variables), as appropriate, using Prism 9 (Graph Pad, San Diego, CA). A p-value of < 0.05 was considered significant for these analyses.

RNAseq Data Processing and Differential Expression Analysis
Low-quality reads (q<10) and adaptor sequences were eliminated from raw reads using bbduk version 37.64 (16). Each sample was aligned to the GRCh38.p9 assembly of the Homo Sapiens reference from NCBI (17) using version 2.6.0c of the RNA-Seq aligner STAR (18). Features were identified from Gencode (V28). Feature coverage counts were calculated using feature Counts (19). Normalization and differential expression analysis was performed with DESeq2 as previously described (20). No biases were apparent on visual inspection of PCA plots and therefore no batch correction was applied. Significant differentially expressed genes between the two groups are those that have an absolute value of fold change ≥ 1.5 and a false discovery rate (adjusted p-value) of < 0.05 (5% FDR). To visualize the data, a volcano plot was constructed by plotting the log2 fold change by the negative log adjusted p-value for each transcript. Focusing on the top 50 differentially expressed transcripts (sorted by absolute value of fold change) a heatmap was generated by hierarchical clustering using the pheatmap R package version 1.0.12.

Pathway Analysis
For pathway analyses, the dataset was uploaded to Ingenuity Pathway Analysis (QIAGEN Inc., https://www.quiagenbioinformatics.com/ products/ingenuity-pathway-analysis) and analyzed using core analysis within the software with the log ratio cutoffs set to -1.5 to 1.5 and the false discovery rate set to 0.05. These thresholds resulted in 254 differentially expressed transcripts analyzed (123 upregulated and 131 downregulated).

Patient Characteristics
A total of 68 patients were enrolled in the parent study. For the purposes of these analyses, 19 patients with immunoparalysis and 13 patients with normal LPS-induced TNFa production capacity were included (Supplementary Figure 1). Patient demographics are in Table 1. All patients were treated with vasoactive support, with a median vasoactive-inotrope score of 11.5. Overall, 84% required invasive mechanical ventilation. Mortality was 3%, and ICU length of stay was a median of 7.6 [4.2, 12.8] days. Half of patients had a complex chronic condition, most commonly developmental delay and/or seizure disorder (n = 5), chronic lung disease and/or tracheostomy (n = 2), chronic kidney disease (n = 2), and inflammatory bowel disease (n = 2). 16% of patients were immune compromised at baseline due to immunosuppressive medications, including two patients who were status post hematopoietic stem cell transplant (both in the non-immunoparalysis group). Patients with immunoparalysis had higher severities of illness and were less likely to have a complex chronic condition or gram-negative infection.

Differentially Expressed Transcripts
Comparing patients with immunoparalysis to those with normal LPS-induced TNFa production capacity, 2,303 transcripts were differentially expressed ( Figure 1). A heat map generated by hierarchical cluster analysis of the top 50 differentially expressed transcripts between groups is displayed in Figure 2. A complete list of differentially expressed transcripts is in Supplementary Table 1.

Pathway Analyses
Canonical pathways identified by Ingenuity Pathway Analysis are in Table 2 and

DISCUSSION
Our results support the hypothesis that transcriptome profiling might identify critically ill septic children with immunoparalysis and provide hypothesis-generating data to inform potential mechanisms of sepsis-induced immune suppression in critically ill children. In our cohort, immunoparalysis was associated with downregulation of pathways involved in antigen presentation by innate immune cells, cell-mediated effector T-cell responses, and natural killer cell function. Pathways involved in platelet activation, leukocyte migration, and corticotropin releasing hormone signaling were upregulated. These data suggest that mechanisms of sepsisinduced immunoparalysis likely involve complex interactions among multiple arms of the immune response. Innate immune cells play a critical role in propagating and regulating adaptive immune responses. Thus, it is not surprising that children with immunoparalysis display downregulation of pathways involved in interactions between innate and adaptive immune cells and cell-mediated immunity. Skewing of helper T cell activity away from cell mediated (Th1) immunity and toward more immune suppressive (Th2) responses has long been recognized in adult sepsis (21)(22)(23). Our previous studies similarly implicate adaptive immune suppression in critically ill septic children (5,24,25). In previous transcriptomic analyses of children with septic shock, pathways involved in antigen presentation and T cell receptor signaling have been associated with adverse clinical outcomes, though assays of immune cell function were not performed in these studies (6,8). Ours is among the first studies to demonstrate an association between innate immune suppression, as assessed by stimulated cytokine production, and downregulation of genes involved in antigen presentation and cell mediated immunity. Major histocompatibility complex (MHC) class II molecules expressed on innate immune cells play a key role in antigen presentation and activation of effector T cells. Multiple MHC class II molecules were downregulated in children with immunoparalysis in our cohort and contributed to pathway analyses that identified T cell signaling pathways. Such downregulation of MHC class II molecules would be consistent with suppression of innate immune cell function in sepsis. For instance, low surface expression of the MHC class II molecule, HLA-DR, on the surface of circulating monocytes has been shown to predict mortality and/or the development of nosocomial infection in septic adults, though its utility in defining clinically relevant immune suppression in critically ill children has been less well studied (26,27).
Specific to T cell function, multiple T-cell receptor subunits and the transcription factor, T-bet, were downregulated in children with immunoparalysis. The selective expression of T-bet is essential in Th1 development, and previous studies have shown that T-bet deficient CD4+ T cells fail to trigger an effective Th1 response during bacterial infections (28)(29)(30)(31). CD8+ T-cells lacking T-bet have also been shown previously to have defects in their effector functions (32). CD8+ T-cells lacking both T-bet and EOMES, which is also downregulated in children with immunoparalysis in our dataset, lose their cytotoxicity (33). EOMES can substitute for T-bet under some circumstances, and together they act through STAT5 to regulate perforin expression in CD8+ T-cells (34,35). The under expression of both T-bet and EOMES may provide some explanation for the downregulation of perforin that was also observed in children with immunoparalysis in our cohort. T-bet and EOMES also cooperate to regulate the expression of CD122 to promote cell survival and proliferation of memory CD8+ T-cells (36). Loss of this cooperative regulation may also help explain the observed downregulation of CD122 in our cohort.
Two pathways involved in natural killer cell function were differentially regulated in our cohort. One of the most downregulated genes in children with immunoparalysis was GZMH, a cytotoxic granule protein expressed by natural killer cells (37)(38)(39)(40)(41). Natural killer cells also require proper expression of EOMES and T-bet for normal terminal differentiation, function, and survival (42)(43)(44). Alterations in natural killer cell function have been implicated in hyperinflammatory sepsis syndromes that resemble macrophage activation syndrome (MAS) in both adults   and children (3,10,45). In a recent multicenter observational study of critically ill children with sepsis-induced organ dysfunction, 18% of children with immunoparalysis also met diagnostic criteria for MAS-like sepsis (3). Though not evaluated in the current study, it is possible that alterations in natural killer cell gene expression and/or function may underly a subset of children with both immunoparalysis and hyperinflammatory MAS-like sepsis induced organ dysfunction. We view this as an important area of ongoing research. Two pathways, corticotropin releasing hormone signaling and GPCR-mediated nutrient sensing in enteroendocrine cells, imply disruptions of endocrine signaling in children with immunoparalysis. The most up-regulated protein-coding gene in our dataset was corticotropin releasing hormone receptor 2 (CRHR2). Corticotropin releasing hormone is chiefly responsible for regulating the hypothalamic-pituitary-adrenal axis. Immune cells have been shown to express CRHR2, with differential effects (46,47). For instance, CRHR2 inhibition was shown to decrease neutrophil-mediated inflammation in a murine stresspneumonia model (48). Conversely, stimulation of CRHR2 on splenic B cells led to an increase in B cell apoptosis (49). These findings mirror our findings of upregulation of neutrophil activation along with downregulation of lymphocyte signaling pathways and supports the hypothesis that immunoparalysis likely represents a mixed state of immune dysregulation characterized by simultaneous systemic inflammation and functional immune suppression. Gai signaling, anther G protein signaling pathway, was also upregulated in children with immunoparalysis. G-protein singling is involved in many aspects of cellular function, including immune cell function. Specifically, Gai has previously been implicated in acute inflammatory responses, again suggesting concomitant patterns of acute inflammation and immune suppression (50).
Both leukocyte extravasation and platelet activation pathways were upregulated in children with immunoparalysis, highlighting the potential importance of non-immune cell types, e.g. endothelial cells and platelets, to inflammation in children with sepsis-induced immunoparalysis. Claudin proteins play important roles in endothelial and epithelial cell tight junction function, with claudin regulation both aiding leukocyte extravasation to sites of infection and contributing to capillary leak in septic shock (51). Their roles in immunoparalysis are uncertain. As discussed above, immunoparalysis is often seen concurrently with severe inflammation. It is possible that disruption in tight junction function drives further inflammation by increasing capillary leak leading to worsened shock and/or by disrupting intestinal barrier function. We view this as an important hypothesis for future research. Along similar lines, matrix metalloproteinases have been shown to contribute to neutrophil chemotaxis and migration and to neutrophil antimicrobial activity (phagocytosis and neutrophil extracellular trap formation) (52). While previous studies identified elevated serum concentrations and/or gene expression of MMP-8 in septic non-survivors compared to survivors (53), juvenile animal studies suggest that MMP-8 serves important roles in bacterial clearance, with increased mortality in MMP-8 null mice (52). These studies suggest that MMP-8 elevations seen in septic patients may be reflective of a necessary inflammatory response to clear infection. Again, our findings underscore the importance of understanding mechanisms of concurrent inflammation and immune suppression and of understanding interactions between endothelial cell and neutrophil regulation in pediatric sepsis. Integrins play important roles in both leukocyte and platelet trafficking, adhesion, and activation. ITGB3 binds fibrinogen, serving as a nidus for platelet/ platelet interactions, and may be an important link between inflammation and coagulation. Similarly, upregulation of the platelet glycoprotein VI pathway among children with immunoparalysis points to a potential role for platelets in dysregulated immune responses to sepsis. This is in keeping with a growing body of literature suggesting that platelets play important roles in the innate immune response, though little is known about relationships between platelet and immune function in septic children (54).
Our study has limitations. Because of our sample size, we are unable to account for differences in baseline characteristics, including the presence of complex chronic conditions, between groups or to determine the extent to which these differences may have influenced our results. Future, larger studies are needed to identify changes in gene expression associated with immunoparalysis independent of other clinical characteristics. Second, differences in absolute numbers of circulating leukocytes could also influence our results. Specifically, children with immunoparalysis had lower absolute monocyte counts compared to children without immunoparalysis, which could explain downregulation in gene expression pathways associated with antigen presentation. Our ongoing work in this area aims to determine the relative contribution of differential cell counts to the diagnosis of immunoparalysis and its related outcomes and to identify cell type-specific changes in gene expression associated with sepsis-induced inflammation and immune suppression in children. Third, our cohort contains 16% of children with baseline immune compromise and it is possible that RNA expression patterns were influenced by underlying rather than sepsisinduced immune suppression. However, these subjects were evenly divided between immunoparalysis and nonimmunoparalysis groups and subjects with underlying immune compromise diagnoses have historically been included in studies of sepsis-induced immunoparalysis and outcomes (3)(4)(5). Future studies are needed to ascribe differences in transcriptomic signatures between endogenous and exogenously produced immunoparalysis in children. It should be noted that because we designed our study to evaluate children with severe innate immune suppression/immunoparalysis versus those with normal LPS-induced TNFa production capacity, we are unable to comment on potential changes in gene expression along the spectrum of TNFa production capacity values. Lastly, our current study was focused on identifying transcriptomic signatures associated with immunoparalysis in the first 48 hours of septic shock, thus we are not able to comment on changes over time. Immune dysregulation in septic shock can by highly dynamic and understanding the trajectory of immune function and associated transcriptomic changes remains an important future research goal. In conclusion, our study provides proof of principle that gene expression data might be useful to identify children with immunoparalysis and identifies several key differentially regulated pathways involved in both innate and adaptive immunity. Our ongoing work in this area aims at dissecting interactions between innate and adaptive immunity in children with septic shock, establishing patterns of both immune function and transcriptome expression over time, and more fulling elucidating patient-specific immunologic pathophysiology, including the use of machine learning algorithms to predict immunologic phenotype and clinical outcomes and to guide individualized immunotherapeutic targets.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: ArrayExpress (https://www.ebi.ac.uk/arrayexpress/), E-MTAB-10938.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Nationwide Children's Hospital Institutional Review Board. Written informed consent to participate in this study was provided by the participants' legal guardian/next of kin.

AUTHOR CONTRIBUTIONS
AS and JM conceived of and planned the study, supervised all of the analyses, interpreted the data, and wrote the manuscript. KJ, JH, and MF processed patient samples and collected clinical data. JF, SW, and AW performed and supervised RNA sequencing and data analysis. KB and MH helped plan the study and interpreted the data. All authors contributed to the planning of the study, data acquisition and interpretation. All authors reviewed the manuscript for critical content, and all authors approved the final version of the manuscript submitted for publication. All authors contributed to the article and approved the submitted version. molecules that are predicted to be inhibited based on the IPA dataset. Orange color represents molecules that are predicted to be activated based on the IPA dataset.
Supplementary Figure 3 | Figure depicting molecules in the Th2 Pathway. Green color represents molecules that are downregulated in the IPA dataset. Red color represents molecules that are upregulated in the IPA dataset. Blue color represents molecules that are predicted to be inhibited based on the IPA dataset. Orange color represents molecules that are predicted to be activated based on the IPA dataset.
Supplementary Figure 4 | Figure depicting molecules in the T cell Exhaustion Pathway. Green color represents molecules that are downregulated in the IPA dataset. Red color represents molecules that are upregulated in the IPA dataset. Blue color represents molecules that are predicted to be inhibited based on the IPA dataset. Orange color represents molecules that are predicted to be activated based on the IPA dataset.
Supplementary Figure 5 | Figure depicting molecules in the Crosstalk between Dendritic Cells and Natural Killer Cells Pathway. Green color represents molecules that are downregulated in the IPA dataset. Red color represents molecules that are upregulated in the IPA dataset. Blue color represents molecules that are predicted to be inhibited based on the IPA dataset. Orange color represents molecules that are predicted to be activated based on the IPA dataset.