Transcriptomic profiles of multiple organ dysfunction syndrome phenotypes in pediatric critical influenza

Background Influenza virus is responsible for a large global burden of disease, especially in children. Multiple Organ Dysfunction Syndrome (MODS) is a life-threatening and fatal complication of severe influenza infection. Methods We measured RNA expression of 469 biologically plausible candidate genes in children admitted to North American pediatric intensive care units with severe influenza virus infection with and without MODS. Whole blood samples from 191 influenza-infected children (median age 6.4 years, IQR: 2.2, 11) were collected a median of 27 hours following admission; for 45 children a second blood sample was collected approximately seven days later. Extracted RNA was hybridized to NanoString mRNA probes, counts normalized, and analyzed using linear models controlling for age and bacterial co-infections (FDR q<0.05). Results Comparing pediatric samples collected near admission, children with Prolonged MODS for ≥7 days (n=38; 9 deaths) had significant upregulation of nine mRNA transcripts associated with neutrophil degranulation (RETN, TCN1, OLFM4, MMP8, LCN2, BPI, LTF, S100A12, GUSB) compared to those who recovered more rapidly from MODS (n=27). These neutrophil transcripts present in early samples predicted Prolonged MODS or death when compared to patients who recovered, however in paired longitudinal samples, they were not differentially expressed over time. Instead, five genes involved in protein metabolism and/or adaptive immunity signaling pathways (RPL3, MRPL3, HLA-DMB, EEF1G, CD8A) were associated with MODS recovery within a week. Conclusion Thus, early increased expression of neutrophil degranulation genes indicated worse clinical outcomes in children with influenza infection, consistent with reports in adult cohorts with influenza, sepsis, and acute respiratory distress syndrome.


METHODS
Patients from the original PICFLU Study (1) were selected to be included in this study using the following inclusion/exclusion criteria: 5. Clinical data was available for pSOFA scoring 6. Parent or legal guardian able and willing to provide permission.

II. Patient Exclusion Criteria
1. Inability to collect respiratory specimen for RT-PCR testing within 7 days of illness onset.
2. Nosocomial-acquired infection as determined at the study site by infection control group.
3. Neuromuscular disease requiring chronic mechanical ventilator support through a mask or tracheostomy for neuromuscular weakness. 4. Chronic mechanical ventilator support through a tracheostomy for chronic respiratory failure.
5. End-stage lung disease being evaluated or awaiting lung transplant.
6. Evidence of current pregnancy from clinical management or other documentation. 7. Any of the following pre-existing conditions: immunodeficiencies, chronic lung disease, symptomatic cardiac disease, neuromuscular disease, malignancy, metabolic or mitochondrial disease, or patients who received systemic immunosuppressive medications within the six weeks prior to admission for this acute illness.

pSOFA Scoring
For the Pediatric Sequential Organ Failure Assessment (pSOFA) respiratory scoring in this study, patients on extracorporeal membrane oxygenation (ECMO) support were assigned the maximum respiratory score of 4. Patients with a PaO2/FiO2 (P/F ratio) less than 200 without respiratory support were given a score of respiratory score 2. For patients not receiving vasopressor support (e.g., mean arterial pressure was not available), systolic blood pressure was used to assign a cardiovascular score of 0 or 1 by applying the 2005 international pediatric sepsis and organ dysfunction guidelines (2). A score of 2 or greater qualified as organ dysfunction for the pSOFA score (range 0-4 per organ system).

Multi-Cohort Meta-Analysis Discovery
To expand on the previously published Influenza MetaSignature (IMS) (3), we performed a systematic search for publicly available human gene expression data in peripheral blood and incorporated newer influenza studies. Using the metaIntegrator package in R (v. 2.0.0) (4, 5) in 7 datasets containing 623 samples, we performed multicohort meta-analysis as described previously (3,(6)(7)(8) resulting in a signature of 152 genes (unpublished data). These are listed as "Influenzaassociated genes" under Gene Source column in Table S1.

RNA Extraction
For RNA extraction, the PAXgene Blood RNA Kit (Qiagen, Waltham, MA) was used with the QIAcube instrument as per manufacturer's guidelines with minor modifications based on blood volume and white blood count (WBC). If WBC count (x10 6 WBC/ml) was obtained close to the blood collection time, it was multiplied by the volume of blood (ml) to estimate the total number of cells. If the estimated total WBC count was between 12x10 6 -27.5x10 6 , then the original protocol was followed as per instructions with 80µl final elution. If the total WBC count was <12x10 6 , then BR5 buffer was omitted from QIAcube setup; instead, after the run was complete, 40 µl of BR5 buffer was added manually to the spin column and centrifuged for 1 min at top speed to elute RNA. The RNA sample was then placed back into the QIAcube to incubate at 65 o C for 5 min as per instructions. If the WBC count was ≥27.5x10 6 , the sample was split evenly into two tubes after the second centrifugation step by resuspending the pellet with 700 µl of BR1 (instead of 350 µl). The QIAcube was set up to use two separate spin columns. If no WBC count was available, the resuspended cell pellet was assessed for opacity (increased opacity = increased cellular material) to determine whether to split the sample or follow the original protocol.
All RNA extracts were tested for quality and quantity using the Agilent 2100 Bioanalyzer and Agilent RNA 6000 Nano Kit (Agilent Technologies, Inc., Santa Clara, CA) per the manufacturer's instructions. RNA samples were then aliquoted and stored at -80 o C until use. The PAXgene Blood RNA Kit incorporated an on-column DNase digestion, with additional DNase treatment applied as needed following quality assessment using the DNA-free™ kit (Ambion, Life Technologies, Carlsbad, CA). RNA analysis was again validated on the Bioanalyzer before use.

NanoString
Reporter and Capture ProbeSets were manufactured within the same lot to reduce variability.
Incorporated into these reagents were the NanoString manufactured 6 positive controls and 8 negative controls so that each sample reported control values. To prepare the RNA for the NanoString run, each sample was diluted to 17 ng/µl in RNase-free water (Qiagen) and 5 µl added into small PCR tubes with 8 µl of a mastermix (containing 3 ul Reporter CodeSet and 5 µl hybridization buffer) followed by 2 µl of Capture ProbeSet. After capping the tubes to mix by inversion, they were pulse centrifuged and incubated for16 hours overnight at 65°C (70°C heated lid) on a thermal cycler (details) then stabilized at 4°C and run within 24 hours. To obtain a loading volume of 30 µl, 15 µl of RNase-free water was added to each sample.

Initial Gene Expression Analysis and Quality Control
The NanoString custom panel included a total of 500 mRNA targets: 469 designed for human mRNA detection and 31 for viral/bacterial detection. The latter were included in the normalization but not reported for analysis as they were outside the scope of this study. During quality control, the housekeeping gene, FPGS, was excluded due to low expression and high coefficient of variation. GUSB was also excluded due to its expression being statistically increased in the Prolonged MODS patient group. Samples were excluded from analysis if they produced values outside of the accepted quality control range for NanoString defined technical cartridge characteristics including binding density, field of view, negative and positive control marker counts, or positive normalization factor, and/or if they were below the limit of detection.
A final 214 influenza positive patients with a single, early sample were included in the normalization.

Statistical analysis for paired longitudinal samples
A method was developed to analyze the paired longitudinal data between two groups by Hahn, et al (9). The 469 genes were assessed for 45 longitudinal paired samples using linear regression with an intercept of the MODS outcome on the normalized mRNA counts per gene for all subjects, where log10 age was included as a single covariate. Contrasts (defined as the difference in measurements between two timepoints) were computed for all subjects falling into the categories "MODS recovery" (group 1) and "Prolonged MODS" (group 2). The null distribution for the intercept of the contrast in the linear model is a t-distribution with n-2 degrees of freedom, where n is the number of subjects to be tested (that is, the size of either group 1 or group 2). We aimed to test the intersection hypothesis consisting of "MODS recovery" and the complement of "Prolonged MODS". To this end, we define two hypotheses. First, for recovered MODS (group 1) we tested the null hypothesis that the intercept is zero against the complementary alternative that the intercept is non-zero. For "Prolonged MODS" (group 2), we tested the alternative hypothesis that the intercept is larger (more extreme) than some level L. For this, we fit a linear regression with intercept to the contrasts with log10 age as covariate and recorded the R² of the resulting model. We then chose L by scaling the intercept for group 2 until we explained 25% of the initial R².
Both the tests for group 1 (under the null) and group 2 (under alternative) yielded one p-value per gene. We evaluate all 469 p-values per group separately. That is, we evaluated the 469 p-values in group 1 using the Benjamini-Hochberg procedure to correct for multiple testing using the FDR criterion (10). Likewise, we applied Benjamini-Hochberg to the 469 p-values for group 2. Each time, the Benjamini-Hochberg procedure was applied using a threshold of 0.05/2. We report those genes that are significant in both groups after the multiple testing correction.

Upregulation of LGALS1, GNA15, DUSP4, and TWISTNB in Prolonged MODS
Although LGALS1, GNA15, DUSP4, and TWISTNB were upregulated in Prolonged MODS patients, their significance was the least robust and therefore not as definitive as the other genes.
Still, they were able to differentiate the patients that recovered from MODS within 7 days as opposed to those that had extended organ dysfunction.
LGALS1 and GNA15 are known as for their roles in immune regulation and homeostasis. The protein Galectin-1 produced by LGALS1 has been shown to bind various influenza A subtypes in vitro and improve survival in infected mice (11) and was suggested as a serum cytokine marker for severe COVID-19 (12). Upregulated GNA15 has been proven to be prognostic (along with 10 other biomarkers) for sepsis in neonates (13). DUSP4 is an immune regulator highly expressed by T regs and also controls T helper cell proliferation and differentiation (14). Overexpression of DUSP4 mRNA has been shown to cause defective signaling of T cell-dependent B-cell responses such as decreased antibody production following influenza vaccination in the elderly (15). In the context of sepsis, DUSP4-deficient mice show improved survival (14). Therefore, the increase of DUSP4 in our Prolonged MODS group in which 24% of the children died reflects its detriment in this cohort. Lastly, upregulation of TWISTNB, which encodes a component of RNA polymerase I that helps catalyze the transcription of DNA into RNA, was expected to reflect lower clinical severity and decreased mortality rate as described in Sweeney, et al's, analysis on bacterial sepsis (16) however this was not the case in our cohort. As a fundamental cellular process, its upregulation reflects increased leukocyte ribosomal RNA production which ultimately helped to differentiate the Prolonged MODS group.