key: cord-0328564-yv57frpw authors: Sacco, K.; Castagnoli, R.; Vakkilainen, S.; Liu, C.; Delmonte, O. M.; Oguz, C.; Kaplan, I. M.; Alehashemi, S.; Burbelo, P. D.; Bhuyan, F.; de Jesus, A. A.; Dobbs, K.; Rosen, L. B.; Cheng, A.; Shaw, E.; Vakkilainen, M. S.; Pala, F.; Lack, J.; Zhang, Y.; Fink, D.; Oikonomou, V.; Snow, A. L.; Dalgard, C. L.; Chen, J.; Sellers, B. A.; Montealegre Sanchez, G. A.; Barron, K.; Rey, E.; Vial, M. C.; Poli, M. C.; Licari, A.; Montagna, D.; Marseglia, G. L.; Licciardi, F.; Ramenghi, U.; Discepolo, V.; Lo Vecchio, A.; Guarino, A.; Eisenstein, E. M.; Imberti, L.; Sottini, A.; Biondi, A.; Mato, S.; Gerstbac, title: Multi-omics approach identifies novel age-, time- and treatment-related immunopathological signatures in MIS-C and pediatric COVID-19 date: 2021-09-27 journal: nan DOI: 10.1101/2021.09.24.21263853 sha: 5c66e55c60c0bccb3fee16abd67bc1ce75841365 doc_id: 328564 cord_uid: yv57frpw Pediatric COVID-19 (pCOVID-19) is rarely severe, however a minority of SARS-CoV-2-infected children may develop MIS-C, a multisystem inflammatory syndrome with significant morbidity. In this longitudinal multi-institutional study, we used multi-omics to identify novel time- and treatment-related immunopathological signatures in children with COVID-19 (n=105) and MIS-C (n=76). pCOVID-19 was characterized by enhanced type I IFN responses, and MIS-C by type II IFN- and NF-{kappa}B dependent responses, matrisome activation, and increased levels of Spike protein. Reduced levels of IL-33 in pCOVID-19, and of CCL22 in MIS-C suggested suppression of Th2 responses. Expansion of TRBV11-2 T-cell clonotypes in MIS-C was associated with inflammation and signatures of T-cell activation, and was reversed by glucocorticoids. The association of MIS-C with the combination of HLA A*02, B*35, C*04 alleles suggests genetic susceptibility. MIS-C B cells showed higher mutation load. Use of IVIG was identified as a confounding factor in the interpretation of autoantibody levels. Following infection with severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), most children develop mild and self-limiting symptoms of coronavirus disease 2019 (COVID-19) [1] although severe cases and fatal outcomes have been also reported [2] . Furthermore, approximately 3-4 weeks after exposure to SARS-CoV-2, some children develop a hyperinflammatory response resembling Kawasaki Disease (KD) and toxic shock syndrome that has been termed multisystem inflammatory syndrome in children (MIS-C) [3, 4] . The mechanisms underlying the phenotypic heterogeneity of pediatric and MIS-C remain ill-defined. Older age, male sex, obesity, co-existing comorbidities genetic defects of TLR3-and TLR7-dependent type I interferon (IFN) pathway, and neutralizing autoantibodies against type I IFNs are associated with more severe clinical outcomes in adults with [5] [6] [7] [8] . More limited information is available on the immune response to SARS-CoV-2 infection in children [9] . Elevated serum levels of several inflammatory biomarkers, an expansion of T cell clonotypes expressing the T-cell receptor TRBV11-2 gene (possibly in response to a SARS-CoV-2 superantigen), and presence of autoantibodies directed against several self-antigens have been reported in MIS-C [10] [11] [12] [13] [14] . 7 To explore the early immune and inflammatory response that characterizes pCOVID-19 and MIS-C, we measured levels of 50 soluble biomarkers in serum or plasma obtained from 55 children with pCOVID-19 within 7 days since onset of symptoms (median, 2 days; IQR, 0.75-3 days), and in 48 children with MIS-C within 7 days from hospitalization (median, 2 days; IQR, 1-4 days). Soluble biomarkers were also measured in 62 MIS-C patients > 7 days after admission (median, 14 .5 days; IQR, 10.75-32 days) and in 53 pHC. Distinctive signatures characterized pCOVID-19 and MIS-C. Higher levels of IFN-α2a were detected in pCOVID-19 (Figure 2A and Figure S1A ). We corroborated these findings using a whole-blood 28-gene NanoString assay capturing expression of type I IFN-stimulated genes in both myeloid and lymphoid cells [17] . High levels of IFN-α2a in pCOVID-19 were associated with a higher type I IFN score ( Figure 2B ). In addition, pCOVID-19 was also characterized by low levels of IL-33 (Figure 2A) , an epithelial and endothelial cell alarmin. To investigate how immune and inflammatory responses evolve over time in pCOVID-19, we analyzed levels of soluble biomarkers in blood samples that were collected in 105 pCOVID-19 patients at various time points after onset of symptoms or positive PCR and observed a rapid decline of IFN-α2a, IFN-γ, IL-10, CXCL10, CCL2 and ferritin ( Figure S2A ). increase in biomarkers related to type II IFN signaling (IFN-γ, CXCL9, CXCL10), macrophage activation (IL-6, sTNFRI, IL-10, sCD25, IL-17, TNF-α, sCD163, CCL2, CCL3, CCL4, ferritin), endothelial injury (VEGF, sVCAM-1/sCD106, sE-Selectin/sCD62E), matrisome-related inflammation [18] (MMP-9, sST2/sIL-33R, CX3CL1) and septic shock (LBP), and low levels of for use under a CC0 license. This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this this version posted September 27, 2021. ; https://doi.org/10.1101/2021.09.24.21263853 doi: medRxiv preprint 8 CCL22. For most of these biomarkers, levels tended to decrease at later time points during hospitalization (Figure 2A and Figure S2B ), likely indicative of inflammation resolution and recovery. Consistent with this broad inflammatory signature, NanoString analysis of a 15-gene type II IFN-dependent and of 11 NF-κB-responsive genes revealed significantly higher scores in paired samples obtained from MIS-C patients at earlier versus later time points during hospitalization (Figure 2C) , and a similar pattern was observed also for type I IFN score ( Figure 2B ). Although some pCOVID-19 patients showed elevation of several inflammatory biomarkers ( Figure 2A and Figure S2B ), their response was generally not as prominent as in MIS-C (Table S1) , and no differences were observed in the NanoString type II IFN and NF-B scores between pCOVID-19 and pHC ( Figure 2C) . Feature importance analysis based on random forest classification (that also included age, sex and ethnicity) identified low levels of IL-33 and increased levels of IL-6, TNF-α, IL-15, ferritin, IL-10, CCL2, IFN-α2a, IFN- and MPO as the most important parameters distinguishing pCOVID-19 from pHC ( Figure 2D ). Using the same approach, elevated levels of ferritin, sST2/sIL-33R, MPO, sTNFRI, IL-6, sVCAM-1, CXCL10, MMP-9, IL-10, IL-15, sCD25, Reg3A, and low levels of CCL22 emerged as the most important parameters distinguishing early MIS-C from pHC ( Figure 2E) . Furthermore, random forest classification identified molecules involved in matrisome (sST2/sIL-33R and MMP-9), intestinal inflammation and myocardial damage (Reg3A) and T cell homeostasis (CCL22) as the most important factors distinguishing MIS-C from pCOVID-19 ( Figure S2C) . Multivariate regression analysis identified IL-33 and CCL22 as the only biomarkers whose levels were significantly different in pCOVID-19 vs. pHC, and in MIS-C vs. pCOVID-19, respectively ( Table S2 ). The prominent inflammatory signature of MIS-C was associated with significantly elevated levels of soluble Spike protein, a pattern that was observed for use under a CC0 license. This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this this version posted September 27, 2021. 9 only in a minority of cases of pCOVID-19 ( Figure 2F) . Finally, anti-Spike (anti-S) and anti-Nucleocapsid (anti-N) antibody levels were significantly higher in MIS-C than in pCOVID-19 ( Figure S2D) , consistent with the limited time interval between onset and symptoms and sample collection in the pCOVID-19 group. To investigate whether age plays an important role in modulating inflammatory responses, we compared levels of soluble biomarkers measured in 8 children and in 26 adults with a similar, moderate degree of COVID-19 severity, whose blood samples were collected within 7 days of symptom onset and within 7 days of admission, respectively. For this purpose, we took advantage of soluble biomarker data that we had previously generated in aCOVID-19 [19] . Groups of pHC (n=53) and adult healthy controls (aHC, n=45) were studied in parallel. For most biomarkers (38/48), blood levels differed between pHC and aHC (Table S1), indicating that age plays an important role in setting baseline immune status. After adjustment for these baseline differences, we identified three groups of biomarkers: a first group in which serum levels were significantly different in pHC and aHC (Figure S1B) , indicating that the difference of unadjusted blood levels observed between pCOVID-19 and aCOVID-19 is probably driven by age, rather than COVID-19 itself. A second group of biomarkers differed significantly in pCOVID-19 versus aCOVID-19, but not between pHC and aHC ( Figure S1C) , suggesting that the nature and severity of inflammatory responses induced by SARS-CoV-2 infection differentially affects patients of different age. 10 To gain additional insights into the broad circulating protein landscape of MIS-C and pCOVID-19 compared to healthy children, we performed proteomic profiling of a subgroup of subjects using SOMAscan® [20] , an aptamer-based proteomics assay, measuring 1305 human protein analytes in plasma samples (Figure 3) . In 10 pCOVID-19 patients, we observed a limited number of up-and down-regulated proteins (26 and 25, respectively) relative to 4 pHC, including increased levels of myeloid activation-associated proteins (MPO, IL18R1, TNFAIP6, ACP5), and SIGLEC7 (an inhibitor of NK cell pyroptosis and inflammasome activation) [21] (Figure 3A 11 with systemic glucocorticoids and IVIG, which has been associated with an improved clinical outcome [16] . However, how these interventions modulate the inflammatory response has not been elucidated. We identified 9 patients for whom we had paired samples available prior to (median, 1 day prior; IQR 1.5 to 0 days) and after (median, +6 days, IQR +2.5 to +10 days) glucocorticoid administration. Biomarkers associated with type II IFN response (IFN-γ, CXCL9), T cell activation (sCD25), cell adhesion (sE-Selectin/sCD62E) and monocyte/macrophage activation (sTNFRII, M-CSF, ferritin, IL-6) decreased following systemic glucocorticoids ( Figure S4A) . We performed correlation analysis among biomarkers, factoring time interval from hospitalization, use of glucocorticoids and/or IVIG, and time interval from treatment, and observed a negative correlation between length of hospitalization and levels of most soluble biomarkers in patients who had received glucocorticoids, irrespective of whether IVIG was administered or not ( Figure S4B -C). Levels of LBP correlated more closely with other inflammatory biomarkers in patients treated with both glucocorticoids and IVIG than in those treated with glucocorticoids only (Figure S4B -C). Furthermore, random forest regression analysis identified LBP as a variable of higher predictive importance in patients who received both glucocorticoids and IVIG (Figure S4D) , suggesting that the addition of IVIG may play a role in controlling signs of septic shock, of which LBP is a biomarker [22] . 12 available for 3 MIS-C patients. We also performed CITE-seq profiling on sorted non-naïve T and B cells to enhance TCR and BCR clonality analysis. We performed unsupervised clustering using the denoised and scaled by background (DSB)-normalized 188 surface protein markers (see Methods) to derive 24 annotated coarser level cell populations ( Figure 4A) . We integrated the CITE-seq data with previously published aCOVID-19 datasets by doing cell-cluster label transfer [24, 25] . The predicted transferred labels are largely concordant with the annotations ( Figure S5A ). As shown in Figure S5B , the frequency of non-classical monocytes was reduced in MIS-C patients; a similar pattern was observed in aCOVID-19, and especially in those with more severe disease [24] . Furthermore, reduced frequencies of plasmacytoid dendritic cells (pDC) were detected in aCOVID-19 (especially in severe disease) and a similar trend was observed in MIS-C. The proportion of pDCs was not altered in pCOVID-19, which could explain the elevated type I IFN responses in pCOVID-19 compared to aCOVID-19 [26] . Another characteristic of pCOVID-19 was the increased frequency of CD8+ memory T cells, which was also noted in adults with less severe disease ( Figure S5B ). We next systemically assessed cell type specific transcriptional changes among pHC, pCOVID-19, and MIS-C ( Figure 4B and Figure S5D ). Linear models incorporating the primary variables of interest, including age and time effects (in MIS-C versus pCOVID-19), were fit for the "pseudobulk" gene expression values for each cell cluster/subset (see Methods). GSEA was then performed based on ranked differential expression genes to search for enriched biological pathways [27] . Strong T and B cell activation signatures were observed in both pCOVID-19 and MIS-C groups compared to pHC, together with increased antigen presentation in both innate and adaptive cell populations ( Figure 4B and Figure S5D ). Consistent with a recent report of increased cytotoxicity of NK cells in MIS-C patients [13] , we also observed enrichment of the for use under a CC0 license. This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this this version posted September 27, 2021. ; 13 gene set "KEGG_Natural_Killer_cell_mediated_cytotoxicity" in CD16 hi NK cells from MIS-C patients ( Figure 4B, left panel) . This enrichment signal was mainly driven by "leading edge"(LE) genes such as KIR2DL3, KIR3DL2, IFNGR1, CD244, LCK, and FASLG. CD16 hi NK cells in pCOVID-19 did not have this signature (Figure 4B , left panel and Figure S5D ). Type I IFN signatures (including gene signatures induced by live viral challenge/vaccination [28, 29] ) were strongly elevated in almost all immune cell subsets in pCOVID-19 but only in a few MIS-C adaptive cell populations; MIS-C exhibited broadly lower type I IFN signatures across cell types compared to pCOVID-19 ( Figure 4B and Figure S5D ). Consistent with our prior CITE-seq analysis in adults [24] , time effect analysis (using CITE-Seq data from samples collected at various time points during hospitalization) hinted that the type I IFN signature in pCOVID-19 decreased over time in most cell types (Figures 4C) . Although classical monocyte cell frequencies were similar, the mRNA based uniform manifold approximation and projection (UMAP) visualization of monocytes showed separation among pHC, pCOVID-19, and MIS-C ( Figure S5C, left panel) . Specifically, MIS-C monocytes showed significantly higher levels of CD163 expression and of several S100A family inflammatory genes; the latter were also increased (although to a lesser degree) in pCOVID-19 monocytes compared to pHC ( Figure S5C, middle and This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this this version posted September 27, 2021. ; 14 [28] , which could partly be due to steroid and IVIG use. Intriguingly, the lymphocytes (CD4+ and CD8+ T cells and NK cells) and DC populations tended to have lower inflammatory signatures instead in pCOVID-19 than both MIS-C and pHC: e.g., pCOVID-19 had lower LE gene signature scores (of the HALLMARK_TNFα_via_NFκB signaling gene set) compared to MIS-C and controls ( Figures 4B, D and E and S5D) . This repressed inflammatory gene signatures in nonmonocyte populations in pCOVID-19 could point to unique response mechanisms in children compared to adults. Bulk high throughput sequencing of TCR (TRB) repertoire was performed to analyze the breadth of the SARS-CoV-2 specific TCR repertoire by identifying exact CDR3 amino acid sequence matches between the SARS-CoV-2-specific sequences previously reported in the ImmuneCODE database and our dataset, as previously described (See Methods). Breadth values represent the fraction of TRB clonotypes that are SARS-CoV-2 specific in each repertoire. By analyzing the TCR repertoire from blood samples obtained at various time intervals from onset of symptoms and hospitalization, we observed a modest increase in breadth of SARS-CoV-2 specific clonotypes in MIS-C and pCOVID-19 compared to pHC (Figure S6A ), indicating that SARS-CoV-2 had not induced an increased proportion of unique SARS-CoV-2-specific CDR3 clonotypes in infected children. High-throughput sequencing analysis of TRB Variable (TRBV) gene usage revealed markedly increased frequency of TRBV11-2 clonotypes in MIS-C (Figure 5A) , confirming previous reports [13, 14, 30, 31] . Interestingly, such increased frequency of TRBV11-2 clonotypes was restricted to MIS-C samples that were collected soon after hospitalization, whereas a rapid decline in the for use under a CC0 license. This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. proportion of TRBV11-2 clonotypes was observed thereafter (Figure 5B ), as also reported by others [31] . Both the increased TRBV11-2 usage and the progressive decline in the frequency of TRBV11-2 clonotypes were confirmed in CITE-seq profiling of CD4+ T cells ( Figure 5C ) and total T cells of MIS-C patients. Computational analysis revealed enrichment of unique SARS-CoV-2 specific CDR3 clonotypes among TRBV11-2 positive clonotypes in all groups (MIS-C, pCOVID-19, pHC); however, such enrichment was significantly lower in MIS-C compared to pCOVID-19 and pHC ( Figure S6B ). In silico structural modeling had previously shown that TRBV11-2 binds through its CDR2 (and not CDR3) region to a superantigen-like motif in the SARS-CoV-2 Spike protein [14, 30, 32] . High-throughput sequencing of TRB repertoire demonstrated that TRBV11-2 clonotypes of MIS-C patients were characterized by a diverse usage of associated TRBJ genes ( Figure S6C ) and a broad distribution of CDR3 length (Figure S6D) , arguing against oligoclonal expansions. We also observed that expansion of TRBV11-2 clonotypes in MIS-C positively correlates with levels of IL-15, IL-17, sCD25, IL-18, sTNFRI, and sTNFRII ( Figure S6E) . A similar correlation between frequency of TRBV11-2 clonotypes and biomarkers of inflammation has been previously reported [14] . Single cell gene expression analysis of these TRBV11-2 CD4+ T cell clones within MIS-C samples showed slightly higher average expression of effector genes like GZMK, PRF1, GZMA, and IL32 compared with other CD4+ MIS-C T cells ( Figure 5D ). In addition, MIS-C This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this this version posted September 27, 2021. ; https://doi.org/10.1101/2021.09.24.21263853 doi: medRxiv preprint expression of T cell co-stimulatory molecules CD28 and CD150 (SLAM) ( Figure 5E ). Furthermore, the transcriptional signature of TRBV11-2 CD4+ T cells was characterized by increased expression of genes involved in apoptosis ( Figure 5F ). Together, these results suggest that TRBV11-2 expressing T cells represent a population of cells poised to respond to activating signals and undergo apoptosis. This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. To investigate the presence of autoantibodies, we used the human proteomic (HuProt TM ) assay comparing 10 MIS-C samples (4 with and 6 without prior IVIG treatment) to 5 pCOVID-19. We detected several autoantibodies in MIS-C, including previously reported TROVE2/Ro60 and ATP4A [14] (Figure 5H ). However, positivity was mostly evident in MIS-C samples drawn after IVIG administration, suggesting that IVIG could be a confounding factor in the evaluation of the for use under a CC0 license. This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. presence of autoantibody in MIS-C. Pre-existing neutralizing autoantibodies targeting IFN-α and/or IFN-ω are frequently detected in critical aCOVID-19 [8] . To investigate whether such autoantibodies are also present in children, we screened serum from pHC (n=53), pCOVID-19 (n=70), and MIS-C (n=40) for autoantibodies against IFN-α, IFN-β, IFN-ω, and IFN-γ. Low levels of immunoreactivity were detected in 1 pHC, 5 MIS-C and 5 p-COVID-19, but none of these samples blocked IFN-α, IFN-β, IFN-ω, or IFN-γ-induced STAT1 phosphorylation in healthy donor peripheral blood mononuclear cells. Defining the pathophysiology underlying distinct phenotypes of SARS-CoV-2 related diseases in children represents an important medical need. Here, by applying a multi-omics approach to a relatively large cohort of patients from multiple centers and of various ethnicities, we have demonstrated that pCOVID-19 and MIS-C have distinctive biomarker signatures, with MIS-C presenting higher and qualitatively different inflammatory activation compared to pCOVID-19. Innate immune responses, and in particular type I IFN-dependent signaling, play a critical role in controlling replication of respiratory tract viruses early after infection [40] . Consistent with this, it has been shown that pre-activated antiviral innate immunity in the upper airways controls early SARS-CoV-2 infection in children [26] . Defective IFN-α responses have been demonstrated in severe aCOVID-19 [41, 42] , and it has been hypothesized that higher IFN-α levels may be responsible for the milder course of pCOVID-19 [43] . Our observations of intact frequencies of pDCs in pCOVID-19, associated with robustly elevated IFN-α2a levels and increased expression of type I IFN-dependent genes in peripheral blood samples collected within 7 days from onset of symptoms, support this hypothesis. None of the pCOVID-19 and MIS-C patients had neutralizing autoantibodies to type I IFN, which are associated with life-threatening aCOVID-19 [8] . However, for use under a CC0 license. This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this this version posted September 27, 2021. ; this issue should be revisited using a more sensitive assay that has been recently described to detect these autoantibodies [8] . The identification of decreased IL-33 and CCL22 levels as distinctive features of pCOVID-19 and MIS-C, respectively, represents a novel finding of our study which needs validation in other cohorts. IL-33 is a member of the IL-1 cytokine family and is released mainly by epithelial cells In our study of MIS-C, consistent with previous observations [11] [12] [13] 15] , we have demonstrated elevated levels of multiple soluble biomarkers associated with recruitment and activation of for use under a CC0 license. This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this this version posted September 27, 2021. ; 20 monocytes and neutrophils, vascular endothelium injury, matrisome activation, gastrointestinal and cardiac involvement, and septic shock. Activation of matrisome which encompasses proteins associated with the extracellular matrix including the endothelium [18] , and increased levels of This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. Consistent with this, a similar breadth of SARS-CoV-2 specific CDR3 clonotypes was also observed in children with MIS-C, several weeks after infection. While no skewing in the usage of TRBV genes was observed in pCOVID-19, we have confirmed and extended previous observations demonstrating increased TRBV11-2 usage in MIS-C, and the polyclonal nature of TRBV11-2 expressing T cells [30, 32] . Structural modeling has suggested that amino acid residues within the CDR2 region encoded by the TRBV11-2 gene may bind to a superantigen-like motif at the C-terminal region of the Spike S1 subunit in a CDR3-independent manner [30, 32] . In this regard, it is interesting that high levels of soluble Spike protein were found This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this this version posted September 27, 2021. ; https://doi.org/10.1101/2021.09.24.21263853 doi: medRxiv preprint 22 these biomarkers decreased within 1-2 weeks after use of glucocorticoids. It is known that glucocorticoids preferentially mediate apoptosis of activated T cells, predominantly through the mitochondrial pathway [58-60]. We hypothesize that this is the mechanism leading to rapid decrease of TRBV11-2 clonotypes following use of glucocorticoids in MIS-C. We have confirmed recent findings that MIS-C is associated with presence of the HLA-A2, C04, B35 alleles [30] . However, in our study, this association was independent of the severity of the disease and of the frequency of TRBV11-2 clonotypes. These observations argue for a genetic basis of susceptibility to MIS-C, but the precise mechanisms involved remain to be defined. and a broad spectrum of autoantibodies directed against gastrointestinal, endothelial, cardiac and immune cell antigens [11] [12] [13] [14] . In our study, we did not observe an increased frequency of IGHV4-34 clonotypes. However, we demonstrated an increased SHM rate in plasmablasts in MIS-C, correlating with increased expression of several activation markers on the cell surface of both memory B cells and plasmablasts. We also confirmed the presence of several autoantibodies. However, autoantibody positivity was largely restricted to samples obtained after IVIG administration, suggesting that use of IVIG is an important confounding factor. Similar This study has some limitations. Only a few children with severe pCOVID-19 were included in the study, and therefore the results obtained may not apply to rare cases of life-threatening disease in children. We did not have longitudinal samples available for all MIS-C patients. However, the for use under a CC0 license. This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this this version posted September 27, 2021. This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this this version posted September 27, 2021. ; 33 MIS-C, multisystem inflammatory syndrome in children; pCOVID-19, pediatric COVID-19; pHC, pediatric healthy controls. for use under a CC0 license. This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this this version posted September 27, 2021. ; 1 The study included a total of 181 pediatric patients (≤ 18 years) with clinically and laboratory confirmed MIS-C (n=76), pediatric COVID-19 (pCOVID-19, n=105) and pediatric healthy controls (pHC, n=60). All subjects were recruited following protocols approved by local Institutional Review Boards (IRBs) or from NIH NIAID IRB, as indicated below: This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this this version posted September 27, 2021. ; https://doi.org/10.1101/2021.09.24.21263853 doi: medRxiv preprint 3 disease was defined as follows: 1) asymptomatic, 2) mild, 3) moderate, 4) severe and 5) critical as per the NIH COVID-19 Treatment Guidelines [3] . Advisory case definition [4] , but only patients with evidence of prior SARS-CoV-2 infection (as determined by positive PCR ± anti-S/anti-N serology) were included. MIS-C patients were divided into moderate (MIS-C-M) and severe (MIS-C-S) groups based on clinical criteria. In particular, as previously reported [5] , MIS-C-S patients were defined if they manifested cardiac and/or pulmonary failure (requiring vasoactive medication and/or significant respiratory support with positive pressure or mechanical ventilation). Demographic, clinical and laboratory features are reported in Table 1 . Cytokine, chemokine and other biomarker analysis was performed on EDTA plasma or sera obtained from patients with pCOVID-19 (n=105), MIS-C (n=73) and pHC (n=53). Among the entire cohort, samples from children with pCOVID-19 within 7 days since onset of symptoms (n=55), and from children with MIS-C within 7 days from hospitalization (n=48) were considered to investigate the early phase on these conditions. Because of limited available volume, patient samples were analyzed as single determinations. Duplicate determinations of control samples and samples from pHC yielded coefficients of variation that were normally <20%. Blood samples were centrifuged, and serum or plasma samples frozen immediately in a -85 o C freezer prior to analysis. Levels of soluble biomarkers whose data were concordant for use under a CC0 license. This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this this version posted September 27, 2021. ; https://doi.org/10.1101/2021.09.24.21263853 doi: medRxiv preprint 4 between both plasma and sera, were measured as previously described [1] . The 50 soluble biomarkers included: interferon (IFN)-γ, IL-10, IL-2, IL-6, IL-8, TNF-α, GM-CSF, IL-12p40, IL-15, IL-16, IL-17, IL-5, LT-α ( TNF-β) This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. SARS-CoV-2 anti-spike and anti-nucleocapsid antibody testing was performed via luciferase immunoprecipitation systems assay, as previously described [6] . For the biomarker values that were below the lower limit of detection (LLOQ) the actual measured concentrations were used or, if unavailable and reported as zeros (for 26 of the 50 biomarkers), values were extrapolated as LLOQ divided by two. The exception was made for the comparison of pCOVID-19 and aCOVID-19, due to the absence of LLOQ for the biomarker measurements in adults. Therefore, only values over zeros were used for that analysis. The univariate analysis comparing biomarkers between or among the subject groups included only values measured in the first seven days from symptom onset for the acute COVID-19 cases or from hospitalization for subjects with MIS-C. The comparison was performed using Mann-Whitney U test (when two groups were compared) or Kruskal-Wallis test (corrected for multiple comparisons) when multiple groups were compared. Biomarkers differing significantly between or among groups were then included in the multivariate model together with age, sex and ethnicity. For the comparison of pCOVID-19 with pHC, allergic conditions (allergic rhinitis, asthma, atopic dermatitis) were also included as a variable in multivariate regression analysis, and. In addition, the analysis was repeated after the exclusion of allergic pHC, with unchanged for use under a CC0 license. This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. Trained with Python sklearn library's RandomForestClassifier object, using parameters: n_estimator=2000, random_state=42 for data set. Results represent the relative importance of the each 53 attribute provided by the model attribute RandomForestClassifier.feature_importances_. Attribute's direction of influence was based on the increase/decrease of its mean values between compared groups. For the comparison of pHC with pCOVID-19, the classification was then repeated after the exclusion of allergic pHC, with similar results. Total RNA was extracted from whole blood samples collected in PAXgene tubes (Qiagen, Germantown, MD). Gene expression of selected genes was determined by NanoString (NanoString Technologies, Seattle, WA) and a 28-gene type I IFN score and an 11-gene NF-κB score was calculated as previously described [8] . An IFN-γ score was calculated based on 15 IFN-γ-regulated genes [9] . Briefly, the 28-gene type I IFN score is the sum of the z-scores of 28 type I IFN response genes, the 11-gene NF-κB score is the sum of the z-scores of 11 NF-κB target genes and the 15-gene IFN-γ score is the sum of the z-scores of 15 response genes. Individual gene z-scores were calculated using the mean and standard deviation of the NanoString counts from pHC [2] . Nonparametric two-tailed Kruskal-Wallis test (corrected for for use under a CC0 license. This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this this version posted September 27, 2021. ; 7 multiple comparison) was used for group comparisons and p values below 0.05 were considered statistically significant. Statistical analyses were performed using GraphPad Prism version 8.00 for Mac OS. SOMAscan®, an aptamer-based proteomics assay, was used to measure 1305 human protein analytes in plasma (SomaLogic, Inc; Boulder, CO, USA). The platform technology is described in Candia et al. [10] . Sample data was normalized to remove hybridization variation within a run. Overall scaling was performed on a per-plate basis to remove overall intensity differences between runs. This was followed by median normalization across the different sample types to remove other assay biases within the run. The statistical analysis of SOMAscan® results was performed using R Studio (R Core Team, 2020) [11] , also using a specifically developed webtool for basic data plotting and analysis [12] . For each group comparison, top up-and downregulated 8 The potential interactions between all variables in the biomarker and timeline data (MIS-C samples only) were characterized by first scaling the values of each variable (with the scale function in R); then, Pearson correlation coefficients and random forest regression based interaction strengths between the variables were computed. The latter approach allowed us to integrate the biomarker levels with the timeline variables in a multivariate setting, while taking into account the potential linear and nonlinear interactions between all variables. Another advantage of this modeling approach is that it mitigates potential biases that could originate from model misspecification since random forest is designed to identify model structures (decision tree ensembles) optimal for any given data set. Pearson correlation coefficients were computed using the corr.test function (psych package in R). Biomarkers and the time interval variables are ordered by hierarchical clustering (with complete linkage) based on their overall correlation patterns that were visualized with the corrplot function (corrplot package in R). Random forest regression models were developed using the GENIE3 algorithm described in Huynh-Thu et al. [14] . Each model was composed of 1000 decision trees that collectively predict a given variable's value using all remaining variables as predictors. GENIE3 algorithm also identifies a predictive importance value of a given predictor in each predictor-target pair, which is also referred to as the interaction strength [14] . To assess the difference in the overall predictive importance (derived by GENIE3) of each variable with and without IVIG treatment, we summed the interaction strengths associated with each predictor-target pair for a given predictor variable in either treatment condition (steroid treatment applied in both conditions). The resulting values were visualized using the Complexheatmap package in R. The variables for use under a CC0 license. This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. sequencing on an Illumina NovaSeq 6000 using a S4 Reagent Kit (300 cycles) using 151+8+8+151 cycle run parameters. Primary sequencing data was demuxed using the Illumina HAS2.2 pipeline and sample-level quality control for base quality, coverage, duplicates, and contamination (FREEMIX < 0.05 by VerifyBamID) was conducted. All sequencing data were then processed with Burrows-Wheeler Aligner (BWA) and the Genome Analysis Toolkit (GATK) best-practice pipeline for alignment and variant call. Samples underwent whole genome sequencing at >=30X median depth. Raw fastq files were trimmed using Trimmomatic v0.39 [15] and mapped to the hg38 human reference genome using BWA-MEM v07.17. PCR Duplicates were marked using Samblaster v0.1.2.5 [16] and GATK4 v4.1.9.0 was used to perform BAM recalibration, and HLA*LA [17] was used to call HLA genotypes. for use under a CC0 license. This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. Hypervariable complementarity determining regions (CDR)3 of T-cell receptor (TCR) β-chains and B-cell receptor (BCR) immunoglobulin heavy (IGH) locus present in the human peripheral blood mononuclear cell (PBMC) samples were sequenced in a high-throughput manner using the immunoSEQ assay [18] after amplification of the extracted DNA in a bias-controlled multiplex PCR. The resulting CDR3 sequences were collapsed and filtered to quantify the absolute abundance and frequency of each unique CDR3 region with Adaptive Biotechnologies' pipeline [19] . Somatic hypermutation (SHM) rate was computed by first matching the germline sequences to IMGT gene identification, flagging the IGH assay mutations (mismatches) to V-gene segments as SHM in the same pipeline. Then, the number of detected SHMs was divided with the number of nucleotides in the region where each SHM set is observed (V gene region) to compute the fraction of clonotypes with >1% SHM rate per nucleotide. The alternative SHM rate metric, which is based on the mutation counts divided by the total number of nucleotides in the V gene region, produced similar results when we compared the SHM rate distributions between the pHC, pCOVID-19, and MIS-C samples. We computed the TCR repertoire statistics, including gene usage, using Immunarch [20] . SARS-CoV-2-specific breadth and depth of each sample was computed using the approach described in Snyder et al. [19] by utilizing the SARS-CoV-2-specific CDR3 sequences previously reported in the ImmuneCODE database [21] . The R package ggpubr was used for visualization of the results with violin, bubble, box, and density plots, whereas the non-parametric Wilcoxon rank sum and Kruskal-Wallis testing and Pearson correlation calculations (along with regression lines showing the 95% confidence for use under a CC0 license. This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this this version posted September 27, 2021. 11 intervals) were also performed with ggpubr. Random forest (RF) regression models were built to compute the interactions between biomarker levels, gene usage, and timeline variables with GENIE3 (GEne Network Inference with Ensemble of trees) [14] using scaled inputs. Frozen PBMC samples were thawed, recovered and washed using RPMI media with 10% FBS and 10mg/mL DNase I (STEMCELL) and then processed as previously described [22] for CITE- For each sample, 100,000-500,000 cells were processed in Trizol using the miRNAeasy micro kit (Qiagen, Germantown, MD) and standard RNA sequencing libraries were generated using for use under a CC0 license. This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this this version posted September 27, 2021. ; 13 Illumina Truseq library preparation kits. The results of bulk RNA sequencing were used for demultiplex of CITE-seq samples by generating SNP calls for each donor. Sequencing results were demultiplexed and converted to FASTQ format using Illumina bcl2fastq software. The sequencing reads were adapter and quality trimmed and then aligned to the human genome using the splice-aware STAR aligner and SNP calls were generated using the previously published protocol [23] . The software package demuxlet [24] was used to then match single cell gene expression data to each donor and identify empty droplets and doublets. Because multiple samples from different timepoints for each donor were collected and could not be demultiplexed by this method alone, 'hashtag' antibodies (Biolegend) were used to uniquely label the different time points. The single cell data processing, CITE-seq protein data denoise and clustering were performed as described before [22] . Specifically, CellRanger (10x Genomics) version 3.1.0 was used to map cDNA libraries to the hg19 genome reference and to count antibody tag features. Data were further processed using Seurat (v.3.1.0) [25] running in R v3.6.1. After filtering to single cell based on demuxlet output, we further demultiplexed the timepoints using the hashtag antibody staining. We removed cells with less than 250 or greater than 4,000 detected genes, greater than 20% mitochondrial reads, cell surface protein tag greater than 200,000, and hashtag antibody counts greater than 50,000. The protein data was normalized and denoised using the DSB This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this this version posted September 27, 2021. ; 14 10, denoise_counts = TRUE, use.isotype.control = TRUE. The DSB-normalized protein data, excluding the isotype control antibodies, were used to generate the Euclidean distance matrix computed for all single cells. Then the shared nearest neighbor (SNN) graph followed by knearest neighbors clustering were built using the FindNeighbors and FindClusters functions in Seurat (v3.1.0), respectively. Major cell clusters were then manually annotated using the surface protein together with gene expression. To compare the cell population frequencies directly with aCOVID-19 patients and avoid potential annotation batch effect, the previously published aCOVID-19 dataset [22] was projected onto CITE-seq data-query from this experiment in Seurat (v3.1.0) [25] using FindTransferAnchors function. Log normalization and first 30 PCs were used for the integration. Cell annotations were then predicted using TransferData function and the predicated labels were added to the metadata as predicated.id column. Pseudobulk gene differential expression analysis and gene set enrichment analysis were performed as described before [22] . Briefly, all unsorted cells in a given sample were computationally "pooled" according to their cluster assignment by summing all reads for a given gene. Pseudobulk libraries made up by few cells and therefore likely not modeled properly by bulk differential expression methods were removed from analysis for each cell-type separately to remove samples that contained fewer than 5 cells and less than 40000 unique molecular identifier counts detected after pooling. Lowly expressed genes were removed for each cell type for use under a CC0 license. This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this this version posted September 27, 2021. ; https://doi.org/10.1101/2021.09.24.21263853 doi: medRxiv preprint individually using the filterByExpr function from edgeR [26] . Differentially expressed genes were identified using the limma voom [27] workflow which models the log of the cpm (counts per million) of each gene. Scaling factors for library size normalization were calculated with the calcNormFactors function with method = "RLE". Genes were ranked using the moderated T statistics for the relevant coefficient from the limma voom model. Enriched gene sets were identified using the pre-ranked GSEA algorithm implemented in the FGSEA R package [https://www.biorxiv.org/content/10.1101/060012v3]. Gene set list used for enrichment assessment was the same list used in Liu et al., see [22] for details of gene sets. Selected pathways shown in figures were manually curated to select gene sets relevant to immunology and often enriched in several cell-types across the various differential expression comparisons. Using the pseudobulk limma voom workflow as described in "Pseudobulk differential expression and gene set enrichment analysis", differentially expressed genes between patient samples (with admission days < 41) and pHC were identified with a model with the following formula in R: ~ 0 + mis-c_vs_pediatric_healthy + age and ~ 0 + pediatric_covid_vs_healthy + age, where patient_vs_healthy is a factor variable with two levels. The contrasts.fit function was then used to compare the estimated means between patient and pediatric healthy controls. Similarly, differentially expressed genes between MIS-C samples (with admission days < 41) and pediatric COVID-19 were identified with a model with the following formula in R: ~ 0 + for use under a CC0 license. This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this this version posted September 27, 2021. ; 16 mis-c_vs_pediatric_covid + days_since_admission + age, where mis-c_vs_pediatric_covid is a factor variable with two levels, time effect was considered using the days_since_admission term. The contrasts.fit function was then used to compare the estimated means between MIS-C and pediatric COVID-19. Differentially expressed genes of MIS-C samples and pediatric COVID-19 samples associated with time respectively, were identified with a model with the following formula in R: ~ days_since_admission + age. The contrasts.fit function was then used to estimated changes associated with disease time course of MIS-C and pediatric COVID-19 respectively. Selected module scores (gene set signature score) representing enriched pathway activities were calculated for each sample as reported before [22] . Specifically, leading edge genes identified by GSEA from the MIS-C versus pediatric COVID-19 model above were used to enhance signal-tonoise ratio and highlight mainly the differences between MIS-C versus pediatric COVID-19. The This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. CellRanger (10x Genomics) version 3.1.0 was used to assemble V(D)J contigs (https://support.10xgenomics.com/single-cellvdj/software/pipelines/latest/algorithms/annotation). For TCR data, the V(D)J assignment and clonotype were from 10x CellRanger output of the filtered_contig_annotations.csv file. For BCR data, V(D)J sequencing contigs from 10x CellRanger output was processed using Immcantation v2.7.0 toolbox (https://immcantation.readthedocs.io/en/latest/index.html). IgBLAST and IMGT germline sequence databases and Change-O package [28] were used for sequence alignment and V(D)J annotations. BCR sequence genotype inference and mutation load quantification were performed with reference to the pipeline from Mathew et al. [29] using the TIgGER package [30] and SHazaM package [28] . The TCR and BCR sequence data, contig assignments and estimated BCR mutation frequencies were integrated with the single-cell RNA-seq Seurat object in the metadata. Autoantibody analysis was performed using HuProt™ v4.0 human protein microarrays and processed by CDI Laboratories (Baltimore, MD). IgG profiling was performed for 15 serum samples from 5 children with acute COVID19 and 10 children with MIS-C, of whom 4 had received IVIG. Briefly, the arrays were blocked and probed with the samples at a 1:1,000 dilution and incubated at room temperature for 1 hour. Then the arrays were washed and probed with Alexa-647-anti-human-IgG (Fc) for signal detection as previously described. Utilizing CDI software, quantile normalization of the raw signal intensities (F635 median for IgG; F532 median for IgA) was performed on all arrays. The data of several proteins that directly bind with secondary antibodies detected through buffer incubation without any serum were excluded (such for use under a CC0 license. This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. 18 as IGHG1, IGHG3 and so on) alongside the controls (such as Rhodamine+IgG64, Anti-human IgG, GST 10ng/ul etc.). The quantile normalized IgG binding intensities of the remaining 23040 protein targets were then visualized using the Morpheus, https://software.broadinstitute.org/morpheus. The T-test was used to compare the different groups and candidates were identified using the following criteria: the variance for the data points was greater than 1000000, the fold-change of average signal intensity was greater than 4 between the 2 groups, and the false discovery rate was <0.5. Plasma samples were screened for autoantibodies against IFN-α, IFN-β, IFN-ω and IFN-γ in a multiplex particle-based assay [31] , in which differentially fluorescent magnetic beads were covalently coupled to recombinant human proteins (2.5 ug/reaction). Beads were combined and incubated for 30 minutes with diluted plasma samples (1 to 100 dilution). Beads were then washed and incubated with PE-labeled goat anti-human IgG (1 ug/mL) for an additional 30 minutes. Beads were washed again, resuspended in assay buffer, and analyzed on a BioPlex X200 instrument. Plasmas with a fluorescence intensity > 1500 were tested for blocking activity. The blocking activity of autoantibodies was determined by assessing STAT1 phosphorylation in healthy control cells following stimulation with the appropriate cytokines in the presence of 10% healthy control or patient plasma. Surface-stained healthy control PBMCs were cultured in serum-free RPMI medium with 10% healthy control or patient plasma and were either left unstimulated or stimulated with 10 ng/mL of IFN-α, IFN-β, or IFN-ω or 400 units/mL of IFN-γ for 15 minutes at 37 o C. Cells were fixed, permeabilized, and stained for intranuclear pSTAT1 for use under a CC0 license. This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this this version posted September 27, 2021. ; 19 (Y701). Cells were acquired on a BD LSRFortessa cytometer, gated on CD14+ monocytes, and analyzed with FlowJo software. for use under a CC0 license. This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. for use under a CC0 license. This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetu pCOVID-19, pediatric COVID-19; MIS-C, multi-system inflammatory syndrome in children; pHC, pediatric healthy controls * Chi-square test except for age comparison ** Kruskall-Wallis test, pairwise comparison was significant only between acute COVID and healthy controls *** pHC subjects had allergy as the only co-morbidity **** pCOVID-19 patients negative for PCR had either positive immunoglobulin M or G for SARS-CoV-2 ALC, absolute lymphocyte count; ALT, alanine aminotransferase; ANC, absolute neutrophil count; CRP, C reactive protein; F, female; IQR, interquartile range; M, male; NA, not applicable; PCR, polymerase chain reaction; PLT, absolute platelet count; SARS-CoV-2, severe acute respiratory syndrome coronavirus 2. for use under a CC0 license. This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. Figure S1 -Higher interferon-⍺2a (IFN-⍺2a) levels in children with mild pCOVID-19, and differences in soluble biomarker levels among pediatric (pCOVID-19), adult COVID-19 (aCOVID-19), and pediatric and adult healthy controls (pHC, aHC) . A) Children with mild acute COVID-19 (n=38) in the first 7 days since symptom onset have significantly higher IFN-⍺2a levels compared to healthy pediatric controls [pHC] (n=15), healthy adult controls [aHC] (n=40), children with MIS-C (both in the first 7 days since hospitalization (n=36) and later in the course of the disease (n=32)), and adults with moderate acute COVID-19 (n=26). B-D) Comparison of soluble biomarkers measured within 7 days of symptom onset in children (n=8) and adults (n=26) with moderate acute COVID-19, as well as pHC (n=53) and aHC (n=45), both unadjusted (left graphs; Kruskal-Wallis test) and adjusted for the baseline differences in healthy subjects of the same age group (right graphs, Mann-Whitney test). B) Biomarkers whose serum levels were significantly different in healthy children and adults, but not in diseased subjects. C) Biomarkers that differed significantly in pediatric vs. adult acute COVID-19, but not between healthy children and adults. D) Biomarkers for which both age and SARS-CoV-2 infection independently contributed to differences in levels in children and adults. Right panels This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. Variables of high predictive importance in both data sets Variables of higher predictive importance in the "steroid only" data set Variables of very low predictive importance in either data set Variables of higher predictive importance in the "steroid+IVIG" data set Before glucocorticoids After glucocorticoids for use under a CC0 license. This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this this version posted September 27, 2021. This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) Infection in Children and Adolescents: A Systematic Review Characteristics and Outcomes of Children With Coronavirus Disease 2019 (COVID-19) Infection Admitted to US and Canadian Pediatric Intensive Care Units Hyperinflammatory shock in children during COVID-19 pandemic Multisystem Inflammatory Syndrome in Children in the United States. Reply Factors associated with COVID-19-related death using OpenSAFELY Immune dysregulation and autoreactivity correlate with disease severity in SARS-CoV-2-associated multisystem inflammatory syndrome in children The autoimmune signature of hyperinflammatory multisystem inflammatory syndrome in children Similarities and differences between the immunopathogenesis of COVID-19-related pediatric multisystem inflammatory syndrome and Kawasaki disease Treatment of Multisystem Inflammatory Syndrome in Children Development of a Validated Interferon Score Using NanoString References 1 An immune-based biomarker signature is associated with mortality in COVID-19 patients Distinct interferon signatures and cytokine patterns define additional systemic autoinflammatory diseases COVID-19) Treatment Guidelines. National Institutes of Health Centers for Disease Control and Prevention Immune dysregulation and autoreactivity correlate with disease severity in SARS-CoV-2-associated multisystem inflammatory syndrome in children Sensitivity in Detection of Antibodies to Nucleocapsid and Spike Proteins of Severe Acute Respiratory Syndrome Coronavirus 2 in Patients With Coronavirus Disease Released 2020. IBM SPSS Statistics for Windows Development of a Validated Interferon Score Using NanoString Technology Systematic identification of type I and type II interferon-induced antiviral factors Assessment of Variability in the SOMAscan Assay R: A language and environment for statistical computing. R Foundation for Statistical Computing Webtool For Basic Data Plotting And Analysis. An interactive open source software tool for a rapid first pass analysis Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles Inferring regulatory networks from expression data using treebased methods Trimmomatic: a flexible trimmer for Illumina sequence data SAMBLASTER: fast duplicate marking and structural variant read extraction HLA*LA-HLA typing from linearly projected graph alignments Comprehensive assessment of T-cell receptor beta-chain diversity in alphabeta T cells Magnitude and Dynamics of the T-Cell Response to SARS-CoV-2 Infection at Both Individual and Population Levels. medRxiv Immunarch, R package version 0.5.5.; An R Package for Painless Bioinformatics Analysis of T-cell and B-cell Immune Repertoire Data A large-scale database of T-cell receptor beta (TCRbeta) sequences and binding associations from natural and synthetic exposure to SARS-CoV-2. Res Sq Time-resolved systems immunology reveals a late juncture linked to fatal COVID-19 Assessment of kinship detection using RNA-seq data Multiplexed droplet single-cell RNA-sequencing using natural genetic variation Differential expression analysis of multifactor RNA-Seq experiments with respect to biological variation voom: Precision weights unlock linear model analysis tools for RNAseq read counts Change-O: a toolkit for analyzing large-scale B cell immunoglobulin repertoire sequencing data Single-cell BCR and transcriptome analysis after influenza infection reveals spatiotemporal dynamics of antigen-specific B cells Automated analysis of high-throughput B-cell sequencing data reveals a high frequency of novel immunoglobulin V gene segment alleles Determination of human anticytokine autoantibody profiles using a particle-based approach We thank the patients, their families, and healthy donors for placing their trust in us.We thank Andrea Catzola, MD, and Luca Pierri, MD, from the Pediatric COVID19 HUB at the This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available preprint (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.