Main

After infection with SARS-CoV-2, most children develop mild and self-limiting symptoms of COVID-19 (ref. 1), although severe cases and fatal outcomes have been also reported2. However, 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 MIS-C3,4,5.

The mechanisms underlying the different picture of pCOVID-19 and MIS-C remain ill-defined. Older age, male sex, obesity, co-existing comorbidities, genetic defects of Toll-like receptor 3 (TLR3)-dependent and TLR7-dependent type I IFN pathways and neutralizing auto-antibodies against type I IFNs are associated with more severe clinical outcomes in adults with COVID-19 (aCOVID-19)6,7,8,9. More limited information is available on the immune response to acute SARS-CoV-2 infection in children10. Elevated serum levels of several inflammatory biomarkers, an expansion of T cell clonotypes expressing the T cell receptor (TCR) TRBV11-2 gene (possibly in response to a SARS-CoV-2 superantigen) and presence of auto-antibodies directed against several self-antigens have been reported in MIS-C11,12,13,14,15.

The magnitude of the inflammatory response in MIS-C correlates with disease severity13,16, and use of glucocorticoids and intravenous immunoglobulins (IVIGs) improves clinical outcome17, whereas limited data are available on the efficacy of biologics, such as IL-1 receptor (IL-1R) and tumor necrosis factor-α (TNF-α) antagonists and tocilizumab18,19. Nevertheless, the temporal trajectory of inflammatory markers in response to treatment during the course of the disease has not been elucidated. In this study, we used a multi-omics approach (with analysis of soluble biomarkers, proteomics, single-cell gene expression profile, T and B cell receptor repertoires and auto-antibodies) to comparatively assess longitudinal changes of innate and adaptive immune responses of pCOVID-19 and MIS-C, and we identified distinct signatures associated with pCOVID-19 and MIS-C that might help define the pathophysiology of these disorders and guide treatment.

Results

Characteristics of the study cohorts

We included a total of 186 pediatric patients (110 with pCOVID-19 and 76 with MIS-C) and 76 pHCs. The demographic, clinical and laboratory characteristics of patients and pHCs are reported in Table 1, and the number of patients analyzed with various assays is outlined in Fig. 1.

Table 1 Demographic, clinical and laboratory features
Full size table
Fig. 1: Study cohort and outline of the multi-omics approach.
figure 1

a, Schematic representation of subject cohorts and workflow, with the number of individuals included in each analysis shown in the table. The figure was created with BioRender.com. b, c, The number of patients with pCOVID-19 (b) and MIS-C (c) analyzed by various combination of assays is shown by vertical bars on the top of the diagrams. The total number of patients analyzed with each assay is indicated by horizontal bars on the right of each panel.

Source data

Full size image

Soluble biomarkers in the early phase of pCOVID-19 and MIS-C

To explore early immune and inflammatory responses, we measured levels of 50 soluble biomarkers in serum or plasma obtained from 57 children with pCOVID-19 within 7 days since the onset of symptoms (median, 2 days, interquartile range (IQR), 1–3 days) and in 48 children with MIS-C within 7 days from hospitalization (‘MIS-C Early’; median, 2 days, IQR, 1–4 days). Soluble biomarkers were also measured in 60 patients with MIS-C more than 7 days after admission (‘MIS-C Late’; median, 14 days, IQR, 10.25–31 days) and in 53 pHCs. Distinctive signatures characterized pCOVID-19 and MIS-C. Higher levels of IFN-α2a were detected in pCOVID-19 (Fig. 2a), especially in children with mild disease (Extended Data Fig. 1a). High levels of IFN-α2a in pCOVID-19 were associated with a higher type I IFN score, as determined by a NanoString assay capturing expression of 28 type I IFN-stimulated genes in both myeloid and lymphoid cells.20 (Fig. 2b). In addition, pCOVID-19 was also characterized by low levels of IL-33, an epithelial and endothelial cell alarmin, and by increased levels of some inflammatory biomarkers, whose levels rapidly declined over time (Fig. 2a and Extended Data Fig. 2a,b). However, NanoString analysis of the expression of 15 type II IFN-dependent and 11 nuclear factor (NF)-κB-responsive genes did not reveal differences between pCOVID-19 and pHC (Fig. 2c).

Fig. 2: Blood biomarker analysis in pCOVID-19 and MIS-C.
figure 2

a, Comparison of serum biomarker levels in children with MIS-C Early (n = 48) (within 7 days since admission) and MIS-C Late (more than 7 days after admission, n = 60), pCOVID-19 (n = 57) within 7 days from symptom onset and pHCs (n = 53). b, Comparison of type I IFN score in paired MIS-C Early and MIS-C Late (n = 11), pHC (n = 12) and pCOVID-19 (n = 15) with elevated (pCOVID-19hi, n = 6) and lower (pCOVID-19low, n = 9) IFN-α2a levels. c, Comparison of NF-κB score and type II IFN score in paired MIS-C Early and MIS-C Late (n = 11), pCOVID-19 (n = 15) and pHC (n = 12). d, Random forest classification comparing pCOVID-19 within 7 days from symptom onset (n = 57) to pHC (n = 53). e, Random forest classification comparing MIS-C Early (n = 48) to pHC (n = 53). f, Random forest classification comparing MIS-C Early (n = 48) to pCOVID-19 within 7 days from symptom onset (n = 57). g, Serum spike protein levels in MIS-C (n = 21), pCOVID-19 (n = 9) and pHC (n = 16). Maxima of box plots in a, b, c and g represent median values, and bars represent IQR. Statistical analysis in a–c and g was performed by Kruskal–Wallis test with adjustment for multiple comparisons. P values are marked as follows: *P < 0.05, **P < 0.01, ***P < 0.001 and ****P < 0.0001.

Source data

Full size image

To investigate whether age plays an important role in modulating inflammatory responses (including attenuated inflammation in pCOVID-19 compared to aCOVID-19), we compared levels of soluble biomarkers measured in moderate forms of pCOVID-19 (n = 9) and aCOVID-19 (n = 26)21 as well as in pHCs (n = 53) and adult healthy controls (aHCs, n = 45). For most biomarkers (38/50), blood levels differed between pHC and aHC (Extended Data Table 1), indicating that age plays an important role in setting baseline immune status. Adjustment for these baseline differences is necessary when interpreting the influence of COVID-19 (Extended Data Fig. 1b–d).

Analysis of MIS-C samples obtained within 7 days of hospitalization in 48 patients demonstrated a significant increase in biomarkers related to type II IFN signaling (IFN-γ, CXCL9 and CXCL10), macrophage activation (IL-6, sTNFRI, IL-10, sCD25, IL-17, TNF-α, sCD163, CCL2, CCL3, CCL4, ferritin and IL-15), endothelial injury and activation (VEGF, sVCAM-1/sCD106 and sE-Selectin/sCD62E), neutrophil activation (MPO and lactoferrin), matrisome-related inflammation (MMP-9, sST2/sIL-33R and CX3CL1) and septic shock (LBP) and low levels of CCL22 (Fig. 2a and Extended Data Fig. 2b). The SARS-CoV-2 polymerase chain reaction (PCR) status around the time of admission had no significant effect on the clinical presentation and on the levels of soluble biomarkers (Supplementary Table 1 and Supplementary Fig. 1). For most biomarkers, levels tended to decrease at later time points (MIS-C Late) during hospitalization (Fig. 2a and Extended Data Fig. 2b), concurrent with clinical improvement. Consistent with this broad inflammatory signature, NanoString analysis of 15-gene type II IFN-dependent and 11 NF-κB-responsive genes revealed significantly higher scores in paired samples obtained from patients with MIS-C at earlier versus later time points during hospitalization (Fig. 2c), and a similar pattern was observed also for type I IFN score (Fig. 2b).

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-α, ferritin, CCL2, MPO, IL-15, IFN-α2a, soluble VCAM-1 (sVCAM-1) and IL-10 as the most important parameters distinguishing pCOVID-19 from pHC (Fig. 2d). Using the same approach, elevated levels of several inflammatory biomarkers, and low levels of CCL22, emerged as the most important parameters distinguishing MIS-C Early from pHC (Fig. 2e). Furthermore, random forest classification identified molecules involved in matrisome (sST2/sIL-33R), intestinal inflammation and myocardial damage (Reg3A) and T cell homeostasis (CCL22) as the most important factors distinguishing MIS-C from pCOVID-19 (Fig. 2f). Multivariate regression analysis identified IL-33 as the only biomarker whose levels were significantly different in pCOVID-19 versus pHC, whereas CCL3 and IL-15 distinguished MIS-C from pHC and pCOVID-19, respectively, with a role also for CCL22 in both cases (Extended Data Table 2). The prominent inflammatory signature of MIS-C was associated with significantly elevated levels of soluble spike protein (Fig. 2g). Of note, among 15 patients in whom spike protein levels higher than 40 pg ml−1 were detected within 7 days after admission, only two tested positive for SARS-CoV-2 mRNA in nasopharyngeal swabs. Finally, anti-spike (anti-S) and anti-nucleocapsid (anti-N) antibody levels were significantly higher in MIS-C than in pCOVID-19 (Extended Data Fig. 2c), consistent with the limited time interval between onset of symptoms and sample collection in the pCOVID-19 group.

Proteomic analysis of immunopathological signatures

To gain additional insights into the inflammatory signature of MIS-C and pCOVID-19, we performed proteomic profiling of a subgroup of patients using SOMAscan22. In ten patients with pCOVID-19, we observed a limited number of upregulated and downregulated proteins (26 and 25, respectively) relative to four pHCs, including increased levels of myeloid activation-associated proteins (MPO, IL18R1, TNFAIP6 and ACP5) and SIGLEC7, an inhibitor of natural killer (NK) cell pyroptosis and inflammasome activation23 (Fig. 3a,b). Gene set enrichment analysis (GSEA) revealed molecular signatures of immune activation, compatible with active SARS-CoV-2 viral infection (Fig. 3a).

Fig. 3: Proteomic analysis in MIS-C compared to pCOVID-19.
figure 3

a, b, Upregulated (a) and downregulated (b) plasma proteins obtained from the comparison between pCOVID-19 (n = 10) and pHC (n = 4). c, d, Top 25 upregulated and downregulated plasma proteins obtained from the comparison between MIS-C (within the first 7 days of hospitalization, n = 16) and pHC (n = 4). e, f, Top 25 upregulated and downregulated plasma proteins obtained from the comparison between MIS-C (within the first 7 days of hospitalization, n = 16) and pCOVID-19 (n = 10). g, Median predictive importance values derived from random forest regression of soluble biomarker values in a group of 101 samples obtained at various time points after hospitalization from 38 patients with MIS-C who received both systemic glucocorticoids and IVIG and in another group of 57 samples from 25 patients with MIS-C who received systemic glucocorticoids only. In each random forest regression model (composed of 1,000 decision trees with one model per target), the predictive importance value for each predictor–target pair is computed using the algorithm described in ref. 63. In a–f, top upregulated and downregulated proteins were identified by selecting all proteins with FDR < 0.05 and P < 0.05 (two-tailed t-test) and then ordering them according to increased or decreased fold changes expressed in a log2 scale. Heat maps show the most significantly enriched pathways for the group comparison, and the statistical significance is expressed as −log(P value). FC, fold change; RF, random forest.

Source data

Full size image

A marked inflammatory profile was observed in patients with MIS-C, with a high number of significantly increased (n = 242) and decreased (n = 158) proteins compared to pHC (Fig. 3c,d). Patients with MIS-C had increased levels of several inflammatory biomarkers (serum amyloid A (SSA1), CRP, ferritin, CXCL10, sST2/sIL-33R and CXCL9) and of B natriuretic peptide (NPPB.1), the latter consistent with cardiac involvement in MIS-C. GSEA showed hyperactivation of the matrisome-associated response. Overall, the inflammatory activation observed in MIS-C appeared to be higher and qualitatively different from pCOVID-19 (Fig. 3e,f).

Longitudinal evolution of blood biomarkers in MIS-C

We hypothesized that the differences in soluble biomarker levels detected at early and later time points during the course of MIS-C (Fig. 2a and Extended Data Fig. 2b) could be due to early intervention with systemic glucocorticoids and IVIG17. However, how these interventions modulate the inflammatory response has not been elucidated. The timeline of initiation of therapeutic intervention with various classes of drugs and blood sampling compared to day of admission in patients with MIS-C is reported in Extended Data Fig. 3. We identified 12 patients for whom biomarker levels were measured both before (median, 0 day; IQR, −1 to 0 days) and after (median, +5 days; IQR, +4 to +7.5 days) glucocorticoid administration. Two of these patients had previously received IVIG, and eight additional patients received IVIG in the interval. Biomarkers associated with type II IFN response (IFN-γ and CXCL9), T cell activation (sCD25), cell adhesion (sE-Selectin/sCD62E) and monocyte/macrophage activation (sTNFRII, M-CSF, ferritin and IL-6) decreased after treatment (Extended Data Fig. 4a). To investigate how rapidly treatment with glucocorticoids and/or IVIG might impact on the inflammatory phenotype, we re-analyzed the MIS-C Early cohort, segregating patients into two groups: those whose blood samples were drawn before (n = 12) or after (n = 36) therapeutic intervention. A significant difference of blood levels between untreated and pre-treated patients with MIS-C Early was observed for four biomarkers (lymphotoxin-α (LT-α), lactoferrin, IL-12p70 and IL-5), and a similar trend was present for several other proteins (Extended Data Fig. 4b). Furthermore, treatment before blood sampling was not among the top ten most important variables when introduced in the random forest regression analysis comparing MIS-C Early versus pHC (Extended Data Fig. 4c). Altogether, these data indicate that treatment did not entirely obscure the hyperinflammatory phenotype that characterizes MIS-C early in the course of the disease. However, longitudinal analysis during the entire course of hospitalization revealed 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 (Supplementary Fig. 2a,b). Random forest regression analysis identified several biomarkers that were of higher median predictive importance in patients who received glucocorticoids (Fig. 3g); concurrent use of IVIG had a more specific effect on IL-1R antagonist (IL-1Ra), MPO, sIL-2Rα, sTNFRI, LBP, sICAM-1, CCL3 and sCD163.

Multimodal single-cell profiling of MIS-C and pCOVID-19

To better understand and compare the cell-type-specific gene expression profile of MIS-C and pCOVID-19, we performed single-cell cellular indexing of transcriptomes and epitopes by sequencing (CITE-seq)24,25 in peripheral blood mononuclear cells (PBMCs) from seven patients with MIS-C, eight patients with pCOVID-19 and seven age- and sex-matched pHCs. Two longitudinal samples were available for three patients with MIS-C. We also performed CITE-seq profiling on sorted non-naive T and B cells to enhance TCR and B cell receptor (BCR) clonality analysis. Unsupervised clustering identified 24 annotated coarser-level cell populations (Fig. 4a). Integrating the CITE-seq data with previously published aCOVID-19 datasets25 yielded largely concordant cell clusters (Extended Data Fig. 5a). The frequency of non-classical monocytes was reduced in patients with MIS-C; a similar pattern was observed in aCOVID-19 and especially in patients with more severe disease25 (Disease Severity Matrix (DSM)_high in Extended Data Fig. 5b). Reduced frequencies of plasmacytoid dendritic cells (pDCs) were detected in MIS-C but not in pCOVID-19 compared to pHC. Another characteristic of pCOVID-19 was the increased frequency of CD8+ memory T cells, which was also noted in adults with less severe disease (DSM_low in Extended Data Fig. 5b).

Fig. 4: Multimodal single-cell profiling of MIS-C and pCOVID-19.
figure 4

a, UMAP visualization of single-cell clusters based on protein expression profiles (see Methods for cell type acronyms). b, GSEA of MIS-C versus pHC (left) and MIS-C versus pCOVID-19 (right) at time points within 40 days of admission. Selected gene sets are grouped into functional/pathway categories. Dot color denotes normalized gene set enrichment score, and size indicates −log10(adjusted P value). P values were from GSEA test of the whole gene sets (Methods) and adjusted using the Benjamini–Hochberg method. The sample size for each group was as follows: MIS-C, n = 8 (two patients with two time points); pCOVID-19, n = 7; and pHC, n = 7. Further details of statistical analysis are described in the Methods. c, GSEA result of pCOVID-19 (top) and MIS-C (bottom) based on the association with time (days since admission), showing only the type I IFN-related response signatures. The sample size for each group was as follows: MIS-C, n = 10 (three patients with two time points); pCOVID-19, n = 8; and pHC, n = 7. d, Heat map of HALLMARK_TNFa_Signaling_via_NFκB gene set in CD4+ memory T cells and classical monocytes. Heat map showing the scaled average mRNA expression (row Z-score) of LE genes from the GSEA analysis of MIS-C versus pCOVID-19. Shared LE genes and selected top LE genes from both cell types are labeled by gene symbol. The shared LE genes are annotated on the right column. Each column represents a sample. Patients are grouped by pHC, pCOVID-19 and MIS-C classes, and columns are ordered by days since admission; also shown are the days since admission of each sample (top of the heat maps). e, Per-sample gene set signature scores of the HALLMARK_TNFα_Signaling_via_NFκB gene set in selected cell populations. Gene set scores were calculated using the Gene Set Variation Analysis (GSVA) of LE genes from the MIS-C versus pCOVID-19 model (Methods). P values shown are adjusted P values from GSEA result in b. Box plot showing the median, first and third quantiles (lower and upper hinges) and smallest (lower hinge − 1.5× IQR) and largest (upper hinge + 1.5× IQR) values (lower and upper whiskers). Sample size was as follows: MIS-C, n = 8 (two patients with two time points), and pCOVID-19, n = 7. See Methods for details of some low representative populations.

Full size image

We next systemically assessed cell-type-specific transcriptional changes among pHC, pCOVID-19 and MIS-C using the cell clusters derived from surface proteins (Fig. 4b and Extended Data Fig. 5c). Strong T and B cell activation signatures and increased antigen presentation in both innate and adaptive cell populations were observed in both pCOVID-19 and MIS-C groups compared to pHC (Fig. 4b and Extended Data Fig. 5c). Consistent with a recent report13, we observed enrichment of the gene set ‘KEGG_Natural_Killer_cell_mediated_cytotoxicity’ in CD16hi NK cells from patients with MIS-C but not from patients with pCOVID-19 (Fig. 4b and Extended Data Fig. 5c).

Type I IFN signatures (including gene signatures induced by live viral challenge or vaccination26,27) were strongly elevated in almost all immune cell subsets in pCOVID-19 but only in a few MIS-C adaptive cell populations and pDCs (Fig. 4b and Extended Data Fig. 5c); MIS-C exhibited broadly lower type I IFN signatures across cell types compared to pCOVID-19 (Fig. 4b). Consistent with our prior CITE-seq analysis in adults25, time effect analysis hinted that the type I IFN signature in pCOVID-19 decreased over time in most cell types (Fig. 4c, top), although we caution that the number of longitudinal samples was small.

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 (Extended Data Fig. 5d, left). 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 (Extended Data Fig. 5d, middle and right). However, classical monocytes from patients with MIS-C showed repressed inflammatory signatures (HALLMARK_TNFα_via_NFκB signaling and HALLMARK_inflammatory response pathways/gene sets) compared to both pCOVID-19 and pHC (Fig. 4b,d,e). Intriguingly, the lymphocytes (CD4+ and CD8+ T cells and NK cells) and dendritic cell (DC) populations tended to have lower inflammatory signatures instead in pCOVID-19 than in both MIS-C and pHC (Fig. 4b,d,e and Extended Data Fig. 5c). This repressed inflammatory gene signatures in non-monocyte populations in pCOVID-19 could point to differences in the systemic immune responses in children compared to adults, as also recently reported by others28.

To validate these observations, we interrogated an independent published cohort with single-cell data13 and observed similarly strong signatures of T and B cells, NK and CD8+ T cell cytotoxicity and enhanced type I IFN response (mainly seen in T and B cell populations) in patients with MIS-C (Supplementary Fig. 3a). The repressed inflammatory signatures of monocytes were also seen in this validation cohort, with overlapping leading edge (LE) genes driving these repressed signatures (Supplementary Fig. 3a,b). We next visually assessed these LE genes from the MIS-C versus pCOVID-19 comparison in our cohort by plotting the cell-type-specific expression heat maps of these genes using data from the validation cohort13. This revealed that these genes indeed tend to have lower expression in classical monocytes in MIS-C compared to pHC, although this trend appeared less significant in memory CD4+ T cells (Supplementary Fig. 3b).

TRBV11-2 usage over time in MIS-C CD4+ T cells

Bulk high-throughput sequencing of TCRβ (TRB) repertoire was performed to analyze the breadth of the SARS-CoV-2-specific TCR repertoire, representing the fraction of TRB clonotypes that are SARS-CoV-2 specific in each repertoire. A modest increase in the breadth of SARS-CoV-2-specific clonotypes was observed in pCOVID-19 and MIS-C compared to pHC (Extended Data Fig. 6a).

Analysis of TRBV gene usage revealed markedly increased frequency of TRBV11-2 clonotypes in MIS-C (Fig. 5a), confirming previous reports13,14,15,29,30. 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 proportion of TRBV11-2 clonotypes was observed thereafter (Fig. 5b), as also reported by others15. 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 (Fig. 5c) of patients with MIS-C. Computational analysis revealed enrichment of unique SARS-CoV-2-specific CDR3 clonotypes among TRBV11-2-positive clonotypes in all groups (MIS-C, pCOVID-19 and pHC); however, such enrichment was significantly lower in patients with MIS-C compared to pCOVID-19 and pHC (Extended Data Fig. 6b). Moreover, TRBV11-2 clonotypes of patients with MIS-C were characterized by a diverse usage of associated TRBJ genes (Supplementary Fig. 4a) and a broad distribution of CDR3 length (Supplementary Fig. 4b), arguing against oligoclonal expansions.

Fig. 5: High-throughput sequencing and CITE-seq analysis of T and B cell repertoire.
figure 5

a, TRBV gene usage in MIS-C (n = 96 samples from 58 patients), pCOVID-19 (n = 21 samples from 21 patients) and pHC (n = 13 samples from 13 individuals). Clonotypes with ambiguous gene assignments are excluded from the figure. For each gene, non-parametric Kruskal–Wallis test with unadjusted P values was used to compare the three groups. NS: P > 0.05 (not significant), *P ≤ 0.05, **P ≤ 0.01, ***P ≤ 0.001, ****P ≤ 0.0001. b, TRBV11-2 gene usage observed in patients with MIS-C within the first 7 days (in blue, n = 36 samples from 35 patients) and at later time points (in yellow, n = 59 samples from 44 patients) during hospitalization. Pearson correlation coefficient (number of days from admission versus TRBV11-2 gene usage) and its P value are shown for both time intervals. The inset plot in the figure provides a comparison between the TRBV11-2 gene usage distributions in these two time intervals and a P value derived from two-tailed Wilcoxon rank-sum test. Box plots show the median, first and third quantiles (lower and upper hinges) and smallest (lower hinge − 1.5× IQR) and largest (upper hinge + 1.5× IQR) values (lower and upper whiskers). c, Upper panel, TRBV11-2 usage (TRBV11-2 ratio among each sample) in CD4+ T cells among three groups (pHC, n = 7; pCOVID-19, n = 7; and MIS-C, n = 8 (two patients with two time points)) within 40 days of admission. P values shown are from two-sided Wilcoxon test between indicated two groups. Lower panel, TRBV11-2 usage frequency in MIS-C CD4+ T cells over time (days since admission, n = 10). Pearson correlation (R) and associated P values are shown. The shaded area represents standard error. Each dot indicates a sample. Box plot elements are the same as in Fig. 4e. d, Mutation quantification of plasmablasts in the three groups (pHC, n = 7; pCOVID-19, n = 8; and MIS-C, n = 7). P values shown were obtained using two-sided Wilcoxon test between indicated two groups. Each dot indicates a cell. Box plot elements are the same as in Fig. 4b. e, Heat map showing auto-antibodies with the highest variance ordered by fold change, using a cutoff of 4 fold change (Methods). Comparisons were made among pCOVID-19 (n = 5), MIS-C that did not receive IVIG (n = 6) and MIS-C after IVIG administration (MIS-C_IVIG, n = 4).

Source data

Full size image

The frequency of TRBV11-2 clonotypes in MIS-C positively correlated with levels of several inflammatory biomarkers (Extended Data Fig. 6c), consistent with previous observations14. Single-cell CITE-seq gene expression analysis showed slightly higher average expression of genes associated with T cell activation (HCST and DUSP2) and effector function (GZMK, PRF1, GZMA and IL32), immune cell synapse and adhesion formation (PSMB1, HAVCR2, SIRPG, CTLA4, RAC2, MSN, ITGB2 and SELL) and IL-2 and IL-15 signaling response pathways (SIRPG, IL2RB and IL2RG) in TRBV11-2 CD4+ T cell clones compared to other CD4+ MIS-C T cells (Extended Data Fig. 6d). Differential expression analysis on the cell surface markers (CITE-seq antibody data) revealed higher expression of T cell co-stimulatory molecules CD28 and CD150 (SLAM) (Extended Data Fig. 6e). Furthermore, the transcriptional signature of TRBV11-2 CD4+ T cells was characterized by increased expression of genes involved in apoptosis and lymphocyte activation (Extended Data Fig. 6f and Supplementary Table 2).

Interestingly, time elapsed from the first use of glucocorticoids negatively correlated with the frequency of TRBV11-2 clonotypes (Extended Data Fig. 6g) and was predictive of reduced TRBV11-2 gene usage over time (Extended Data Fig. 6h). This suggests that the use of glucocorticoids might have contributed to the apoptosis transcriptional signatures that we saw in the single-cell analysis above (Extended Data Fig. 6f), even though this could also reflect the contraction of CD4+ T cell subsets during the resolution of disease31,32,33.

It was previously shown that patients with MIS-C with a severe clinical phenotype and expansion of TRBV11-2 share the combination of human leukocyte antigen (HLA) class I alleles A*02, B*35 and C*04, indicating a possible contribution of HLA-mediated restriction in the process of TRBV11-2 expansion29. To determine whether a similar association was found in our patient cohort and to avoid confounding factors due to different frequencies of HLA alleles in different ethnic groups, we analyzed a subcohort of Italian patients only (MIS-C, n = 9; pCOVID-19, n = 64; pHC, n = 44), which we further restricted by selecting individuals of predicted European ancestry (MIS-C, n = 7; pCOVID-19, n = 45; pHC, n = 35). The A*02, B*35 and C*04 allele combination was present in five of the seven patients with MIS-C, in two of the 35 pHCs and in none of the 45 patients with pCOVID-19, reaching statistical significance (Extended Data Table 3). Of note, the combination of HLA A*02, B*35 and C*04 alleles was not associated with severity of the MIS-C phenotype, as it was found in four patients with moderate disease and one child with severe disease.

B cell activation and repertoire in MIS-C

Previous studies have documented B cell abnormalities in severe aCOVID-19 and in MIS-C, with increased number of IgD−CD27−CD11c+ cells in the former34 and increased number of plasmablasts in both conditions34,35, along with increased frequency of IGHV4-34 and IGHV4-39 clonotypes14,34 and presence of auto-antibodies against a variety of self-antigens12,13,14. High-throughput sequencing of the IGH repertoire in 13 pHCs, 15 patients with pCOVID-19 and 19 patients with MIS-C did not reveal major differences in the usage of IGHV genes (Extended Data Fig. 7a,b) but demonstrated an increased rate of somatic hypermutation (SHM) among MIS-C IGHV clonotypes (Extended Data Fig. 7c). CITE-seq analysis revealed a significantly increased frequency of SHM in plasmablasts in MIS-C compared to pCOVID-19 (Fig. 5d), and a similar trend was observed in memory B cells (Extended Data Fig. 7d). Several surface markers associated with B cell activation correlated with mutation frequencies within memory B cells (lower IgD, CD305 and IgM and higher CD27, CD95, CD71 and CD99; Extended Data Fig. 7e) and plasmablasts (CD95, CD99 and HLA-DR; Extended Data Fig. 7f).

To investigate the presence of auto-antibodies, we used the human proteomic (HuProt) assay comparing ten MIS-C samples (four with and six without prior IVIG treatment) to five pCOVID-19 samples. We detected several auto-antibodies in MIS-C, including previously reported TROVE2/Ro60 and ATP4A14 (Fig. 5e). However, positivity was mostly evident in MIS-C samples drawn after IVIG administration, suggesting that IVIG might represent an important confounding factor in the evaluation of the presence of auto-antibody in MIS-C. Pre-existing neutralizing auto-antibodies targeting IFN-α and/or IFN-ω are frequently detected in critical aCOVID-19 (ref. 9). To investigate whether such auto-antibodies are also present in children, we screened serum from pHC (n = 53), pCOVID-19 (n = 70) and MIS-C (n = 40). Borderline levels of positive immunoreactivity against IFN-α and/or IFN-ω were detected in a few patients with MIS-C and pCOVID-19, and no neutralizing activity was detected (Supplementary Fig. 5).

Discussion

Defining the pathophysiology underlying distinct SARS-CoV-2-related diseases in children represents an important medical need. Type I IFN-dependent responses play a critical role in controlling replication of respiratory tract viruses early after infection36. Defective type I IFN responses have been shown in severe aCOVID-19 (refs. 37,38). 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, contrast with findings in aCOVID-19 and are consistent with the demonstration that pre-activated antiviral innate immunity in the upper airways controls early SARS-CoV-2 infection in children28,39.

We identified reduced induction of systemic inflammatory responses as another important feature distinguishing pCOVID-19 versus aCOVID-19 (refs. 21,25), as shown by lower levels of inflammatory biomarkers and decreased transcriptional inflammatory signatures of lymphocyte and DC populations in the former.

The identification of decreased IL-33 levels in pCOVID-19 represents a finding that needs validation in other cohorts. IL-33 is a member of the IL-1 cytokine family and is released mainly by epithelial cells upon infection, cell damage or exposure to allergens40,41. High IL-33 levels are increased in children with severe viral and bacterial infections42,43,44,45. The low IL-33 levels detected in pCOVID-19 might be indicative of modest respiratory epithelium cell damage, whereas high levels of IL-33 were previously demonstrated by our group in critical, but not moderate, aCOVID-19 (ref. 21).

Analysis of soluble biomarker levels in MIS-C revealed low levels of CCL22, a homeostatic chemokine that promotes regulatory T cell migration and function46. By dampening regulatory T cell responses, low CCL22 levels in MIS-C might favor uncontrolled inflammation. Notably, both IL-33 and CCL22 are involved in Th2 responses47,48, and both are negatively regulated by IFN-γ49,50. Along with increased levels of IFN-γ in MIS-C (and, to a lesser extent, in pCOVID-19), these observations indicate that pCOVID-19 and MIS-C are characterized by prominent Th1 and suppressed Th2 responses.

Consistent with previous observations11,12,13,16, we showed that patients with MIS-C had elevated levels of soluble biomarkers associated with recruitment and activation of 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 endothelium51, and increased levels of biomarkers indicative of endothelial cell damage in MIS-C mirror what is observed in various vasculitides, including KD52.

In addition, CITE-seq analysis revealed an MIS-C monocyte signature characterized by increased expression of several members of the S100A family of alarmins and of the scavenger receptor CD163. However, in comparison to pCOVID-19, MIS-C monocytes had lower type I IFN and NF-κB/inflammatory signatures and repressed antigen presentation genes, which were phenotypically similar to the MS1 monocyte cell state reported in severe aCOVID-19 and in bacterial sepsis53,54. These reduced inflammatory signatures of monocytes in MIS-C might have been contributed by the routine administration of glucocorticoids and IVIG early in the course of the disease.

In our study, elevated levels of soluble spike protein were detected in 15 of 21 patients with MIS-C. A previous study correlated elevated spike protein levels in MIS-C to persistence of SARS-CoV-2 in the gastrointestinal tract55. Although we did not investigate the presence of SARS-CoV-2 mRNA in stool samples, only two of these 15 patients with MIS-C had a positive PCR on nasopharyngeal swab within 7 days after admission, indicating that elevated spike protein levels were not due to persistent respiratory tract infection.

Analysis of the T cell and B cell repertoires revealed other important features of pCOVID-19 and MIS-C. The modest increase in the breadth of SARS-CoV-2-specific CDR3 clonotypes in children with pCOVID-19 and MIS-C compared to pHCs is consistent with previous studies showing that younger individuals have pre-existing CD4+ T cells to human endemic β-coronaviruses that are cross-reactive to SARS-CoV-2 spike protein56,57 and that might help contain virus replication, limiting the development of a larger pool of newly generated SARS-CoV-2-specific T cells in infected children.

We confirmed previous observations showing an expansion of TRBV11-2-positive polyclonal T cells in MIS-C, possibly driven by a superantigen-like motif within the C-terminal region of the spike S1 subunit15,29,30,58. TRBV11-2 CD4+ T cells expressed high levels of CD150 and CD28 on their surface, and their transcriptional profile was characterized by expression of genes involved in cell adhesion and activation and of the mitochondrial pathway of apoptosis. Together, these results suggest that TRBV11-2-expressing T cells represent a cell population poised to respond to activating signals and undergo apoptosis. The proportion of TRBV11-2 clonotypes positively correlated with levels of various inflammatory biomarkers, and both the frequency of TRBV11-2 clonotypes and levels of most of these biomarkers decreased within 1–2 weeks after use of glucocorticoids. We postulate that the rapid decrease of TRBV11-2 clonotypes was contributed by the use of glucocorticoids, which are known to mediate apoptosis of activated T cells, predominantly through the mitochondrial pathway59,60,61.

Notably, by selecting patients of homogeneous predicted ancestry, we validated the recent demonstration of the association of MIS-C with the combination of the HLA-A*02, B*35 and C*04 alleles29, arguing for a genetic basis of susceptibility to MIS-C.

Analysis of the B cell compartment of patients with MIS-C showed an increased SHM rate in plasmablasts, correlating with increased expression of several activation markers on the cell surface of both memory B cells and plasmablasts. On the other hand, although auto-antibodies have been reportedly detected in patients with MIS-C also before IVIG administration11,12,13,14, we detected them at higher frequency in samples collected after IVIG administration, indicating that use of IVIG is an important confounding factor. Similar observations have been recently obtained in KD62.

This study has some limitations. Only a few children with severe pCOVID-19 were investigated, and no cases of acutely ill children with conditions other than COVID-19 were included. The transcriptional signature of PBMCs was analyzed in a limited number of patients. Nonetheless, we were able to detect early and late signatures of the disease, and the characteristic gene expression profile identified in our cohort correlated with what has been observed by others13. The vast majority of patients with MIS-C received treatment with glucocorticoids (alone or in association with IVIG) shortly upon hospitalization, so it was not possible to define the relative role of therapeutic interventions and natural history of the disease on the dynamic changes of biomarkers analyzed. However, we postulate that timely therapeutic intervention played a critical role in facilitating resolution of inflammatory complications and favorable clinical outcome in all patients included in the study. Too few patients received IVIG alone (n = 4) or various biologics (n = 12) to allow definition of the specific effects of these treatments. Finally, all blood samples were collected at the time when only the ancestral Wuhan strain, the B1.177 (European lineage) variant and the 1.1.7 (Alpha) variant were circulating at the centers where the patients were enrolled. Therefore, the effect of the Delta and Omicron variants on innate and adaptive immune responses in children with pCOVID-19 and MIS-C remains to be studied.

Relatively few studies have explored immune responses to SARS-CoV-2 in children, most often in a limited number of patients. By applying a multi-omics approach to a large cohort of patients, we have shown important differences in the response to acute SARS-CoV-2 infection in children and adults and established that pCOVID-19 and MIS-C have distinctive immunopathological signatures, which might help better characterize the pathophysiology of these disorders and guide optimal treatment.

Methods

Statistics and reproducibility

This was a natural history study of consecutive cases of patients with pCOVID-19 and MIS-C enrolled at the referring institutions. Informed consent was provided by the parents or guardians and assent by the minor, when appropriate. No statistical method was used to predetermine sample size. Investigators analyzing biomarker levels were blinded to the characteristics of the patients from whom the blood samples had been obtained.

Study population

The study included 186 pediatric patients (≤18 years old) with clinically and laboratory confirmed MIS-C (n = 76) and pCOVID-19 (n = 110) and pHCs (n = 76), whose blood samples were collected between 30 March 2020 and 8 February 2021, upon informed consent and according to protocols approved by local institutional review boards (IRBs): Comité Ético Científico Facultad de Medicina Clínica Alemana Universidad del Desarrollo, Santiago, Chile (protocol 2020-41); Ethics Committee of the Fondazione IRCCS Policlinico San Matteo, Pavia, Italy (protocol 20200037677); Comitato Etico Interaziendale A.O.U. Città della Salute e della Scienza di Torino, Turin, Italy (protocol 00282/2020); Ethics Committee of the University of Naples Federico II, Naples, Italy (protocol 158/20); Comitato Etico Provinciale, Brescia, Italy (protocol NP-4000); University of Milano Bicocca-San Gerardo Hospital, Monza, and Ethics Committee of the National Institute of Infectious Diseases ‘Lazzaro Spallanzani’, Italy (protocol 84/2020); Hadassah Medical Organization IRB, Jerusalem, Israel (protocol HMO-235-20); and National Institute of Allergy and Infectious Diseases (NIAID), National Institutes of Health (NIH), Bethesda, Maryland, USA (protocols NCT04582903, NCT03394053 and NCT03610802).

Clinical datasets from international sites were translated, checked for consistency and transformed to the same scale and units as needed using Python libraries (NumPy, pandas and dateutil), and outliers were manually reviewed. The data harmonized across all sites were collected in LabKey (LabKey Server, Enterprise Edition, version 21.11.4) where final curation was performed by the clinical research team.

The severity of pCOVID-19 was defined as follows: (1) asymptomatic, (2) mild, (3) moderate, (4) severe and (5) critical, as per the NIH COVID-19 Treatment Guidelines64. The clinical severity was not affected by age, sex or ethnicity, and there were no fatal outcomes.

MIS-C diagnosis was based on the Centers for Disease Control and Prevention Health Advisory case definition5, but only patients with evidence of prior SARS-CoV-2 infection (as determined by positive PCR with or without anti-S/anti-N serology) were included. Patients with MIS-C were divided into moderate (MIS-C-M; n = 52, 68%) and severe (MIS-C-S; n = 24, 32%) groups, as previously described13. All patients with MIS-C improved markedly during the hospitalization and were eventually discharged.

For the comparison of pCOVID-19 and aCOVID-19, we used previously published data from our group on biomarkers in aCOVID-19 (ref. 21) as well as a cohort of healthy adults. For NanoString and spike protein levels, pHC samples were obtained from a cohort of healthy children studied by NIAID Translational Autoinflammatory Disease Studies.

Measurement of soluble biomarkers

Analysis of soluble biomarker levels was performed on plasma or serum obtained from patients with pCOVID-19 (n = 110), patients with MIS-C (n = 73) and pHCs (n = 53), including 57 patients with pCOVID-19 and 48 patients with MIS-C whose samples were obtained within 7 days since the onset of symptoms or hospitalization, respectively. Because of limited available volume, patient samples were analyzed as single determinations. Duplicate determinations of samples from pHCs yielded coefficients of variation that were normally less than 20%. Blood samples were centrifuged, and serum or plasma samples were immediately frozen at −85 °C before analysis. Levels of soluble biomarkers whose data were concordant between both plasma and sera were measured as previously described21. Depending on the nature of the analyte, measurements were obtained using the V-PLEX Human Cytokine 30-Plex Kit (Meso Scale Discovery) and analyzed on a MESO QuickPlex SQ 120 reader (Meso Scale Discovery) or using a customized, magnetic bead-based, multiplex assay (R&D Systems), according to the manufacturers’ specifications for standards and dilutions, and the magnetic beads were analyzed on Bio-Plex 3D instrumentation (Bio-Rad). Standard curves were analyzed using non-linear curve fitting, and unknowns were calculated based on the derived equation. Samples that exceeded the highest standards were reanalyzed at higher dilution until the values fell within the range of the known standards. Two control plasma samples and a control sample spiked with a known quantity of each analyte were analyzed on each plate to assess the inter-plate variation and to determine the effect of the biological matrix on the measurement of each analyte. For most analytes, the control samples had less than 25% variation from plate to plate, and the recoveries were generally more than 70%.

For the biomarker values that were below the lower limit of quantification (LLOQ), the actual measured concentrations were used, or, if unavailable and reported as 0 (for 26 of the 50 biomarkers), values were extrapolated as LLOQ divided by 2. The exception was made for the comparison of pCOVID-19 and aCOVID-19, owing to the absence of LLOQ for the biomarker measurements in adults. Therefore, only values over 0 were used for that analysis.

The univariate analysis of biomarker levels 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 and atopic dermatitis) were also included as a variable in multivariate regression analysis. These analyses were completed with IBM SPSS Statistics version 27 and GraphPad Prism version 9 software.

For the random forest classification, we used Python version 3.8.10 and the following libraries: pandas 1.1.2, numpy 1.18.5, scikit-learn 0.23.2 and matplotlib 3.3.2. Three models were trained with 53 attributes: (training set size/validation set size/accuracy) pHC versus MIS-C (n = 78/n = 20/95%), MIS-C versus pCOVID-19 (n = 82/n = 21/100%) and pHC versus pCOVID-19 (n = 87/n = 22/100%), trained with Python sklearn library’s RandomForestClassifier object, using the following parameters: n_estimator = 2,000 and random_state = 42 for dataset. Results represent the relative importance of each of the 53 attributes 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.

Spike protein measurement

Patient serum was collected and analyzed for the concentrations of spike protein using the COVID-19 S-Protein (S1RBD) ELISA kit (ab284402, Abcam). Recombinant SARS-CoV-2 S1 + S2 ECD (S-ECD) protein (RP01283LQ, ABclonal) was spiked at increasing concentrations into pre-COVID serum from healthy controls and was used as standard for the calculation of the spike protein concentration. Pre-COVID-19 pediatric (n = 7, age 7–18 years) and adult (n = 9, age 19–63 years) serum samples were used as controls.

SARS-CoV-2 antibody testing

SARS-CoV-2 anti-S and anti-N antibody testing was performed via luciferase immunoprecipitation systems assay, as previously described65.

NanoString assay

Total RNA was extracted from whole blood samples collected in PAXgene tubes (Qiagen). Gene expression of selected genes was determined by NanoString (NanoString Technologies), and 28-gene type I IFN and 11-gene NF-κB scores were calculated as previously described20. An IFN-γ score was calculated based on 15 IFN-γ-regulated genes66. In brief, 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. Non-parametric two-tailed Kruskal–Wallis test (corrected for multiple comparison) was used for group comparisons, and P values less than 0.05 were considered statistically significant. Statistical analyses were performed using GraphPad Prism version 8.00 for Mac OS.

SOMAscan proteomic discovery platform analysis

SOMAscan (SomaLogic), an aptamer-based proteomics assay, was used to measure 1,305 human protein analytes in plasma. The platform technology is described in Candia et al.22. Sample data were 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), also using a specifically developed webtool for basic data plotting and analysis63. For each group comparison, top upregulated and downregulated proteins were identified by selecting all the proteins with false discovery rate (FDR) < 0.05 and P < 0.05 and then ordering them according to increased or decreased fold change, expressed in a log2 scale. Pathway enrichment analysis was performed on differentially expressed biomarkers among the groups (pCOVID-19, MIS-C and pHC), using the Molecular Signatures Database version 7.4, part of the GSEA software, a joint project of the University of California, San Diego and the Broad Institute.

Biomarker interaction analysis

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 non-linear interactions between all variables.

Pearson correlation coefficient values were computed using the corr.test function (psych package in R). Biomarkers and the time interval variables were 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 built to compute the interactions among biomarker levels, gene usage and timeline variables with GENIE3 (Gene Network Inference with Ensemble of trees)67 using scaled inputs. Each model was composed of 1,000 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 strength67. The median predictive importance value (derived by GENIE3) was extracted from the importance distribution associated with each predictor versus all its targets in either treatment condition (glucocorticoids alone or glucocorticoids + IVIG). The resulting values were visualized using the pheatmap and Complexheatmap packages in R. The variables were clustered based on the median interaction strength (or predictive importance) per variable, by implementing agglomerative hierarchical clustering with Euclidean distance and average linkage.

HLA typing

Genomic DNAs were extracted from patients’ whole blood using the QIAsymphony DNA Midi Kit and quantified using a fluorescence dye-based assay (PicoGreen dsDNA reagent) by a microplate reader (Molecular Devices SpectraMax Gemini XS). Whole-genome sequencing libraries are generated from fragmented DNA using the Illumina TruSeq DNA PCR-Free HT Library Preparation Kit with minor modifications for automation (Hamilton STAR Liquid Handling System) and IDT for Illumina TruSeq DNA UD Indexes (96 Indexes, 96 Samples) adapters. Sequencing libraries were quantified using the KAPA qPCR Quantification Kit (Roche Light Cycler 480 Instrument II) and combined as 24-plex pools after normalization and sequencing on an Illumina NovaSeq 6000 using a S4 Reagent Kit (300 cycles) using 151+8+8+151 cycle run parameters. Primary sequencing data were 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 and the Genome Analysis Toolkit (GATK) best practice pipeline for alignment and variant call. Samples underwent whole-genome sequencing at ≥30× median depth. Raw FASTQ files were trimmed using Trimmomatic version 0.39 (ref. 68) and mapped to the hg38 human reference genome using BWA-MEM version 07.17. PCR duplicates were marked using SAMBLASTER version 0.1.2.5 (ref. 69); GATK4 version 4.1.9.0 was used to perform BAM recalibration; and HLA*LA70 was used to call HLA genotypes. Ethnicity was computed from whole-genome sequencing data by Peddy using 2,504 genome samples from The 1000 Genomes Project as background.

Bulk TCR and BCR repertoire

The CDR3 regions of TRB and IGH rearrangements present in PBMC samples were sequenced in a high-throughput manner using the immunoSEQ assay 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 the Adaptive Biotechnologies pipeline71.

SHM rate was computed by first matching the germline sequences to IMGT gene identification and flagging the IGH assay mutations (mismatches) to V-gene segments as SHM in the same pipeline. Then, the number of detected SHMs was divided by the number of nucleotides in the region where each SHM set is observed (V gene region) to compute the fraction of clonotypes with greater than 1% SHM rate per nucleotide.

We computed the bulk TCR and BCR repertoire statistics, including gene usage, using immunarch72. Gene usage was defined as the fraction of unique clonotypes per sample in which a given gene is present. SARS-CoV-2-specific breadth and depth of each sample was computed using the approach described in Snyder et al.71 by using the SARS-CoV-2-specific CDR3 sequences previously reported in the ImmuneCODE database73.

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 intervals) were also performed with ggpubr. The reported P values and significance levels are based on two-tailed testing.

CITE-seq experimental methods

Single-cell CITE-seq processing

Frozen PBMC samples were thawed, recovered and washed using RPMI media with 10% FBS and 10 mg ml−1 of Dnase I (STEMCELL Technologies) and then processed as previously described25 for CITE-seq staining. In brief, samples from different donors were pooled, and different time points from the same donor were pooled separately so that each pool contains only one time point from one donor. PBMC pools were Fc blocked (Human TruStain FcX, 1:10 dilution, BioLegend) and stained with biotinylated SARS-CoV-2 S1 protein (0.4 µg, Acro Biosystems), TotalSeq-C human ‘hashtag’ antibodies (1:100 dilution, BioLegend) and TotalSeq-C PE Streptavidin (1:500 dilution, BioLegend) and then washed with staining buffer (2% BSA in PBS). A fraction of the combined cells was used for sorting non-naive T and B cells (see below). For the unsorted cell fraction, hashtagged PBMC pools were combined, and cells were stained with a cocktail of TotalSeq-C human lyophilized panel (BioLegend) of 188 surface proteins (plus four isotype controls; see repository file 10, 50-µl reconstitution for 1 million cells staining). Then, cells were washed, resuspended in PBS and counted before proceeding immediately to the single-cell partition step.

Sorting of non-naive B and T cell populations

Pooled PBMC samples from different donors were washed with PBS and incubated with Zombie Red Fixable viability dye (1:1,000 in PBS, BioLegend) for 20 min at 4 °C protected from light. Then, cells were washed with flow staining buffer (10% FBS in PBS) and Fc blocked (Human TruStain FcX, BioLegend) for 15 min on ice. The fluorescence-labeled antibody cocktail against human CD45 (APC/Cyanine7, CD3 (AF488), CD19 (APC), CCR7 (BV786), CD95 (BV650), IgD (PerCP-Cy5.5) and CD27 (PE/Cyanine7); all antibodies were obtained from BioLegend, and all were used at 1:20 dilution) was added at the end of blocking and incubated for 20 min at 4 °C in the dark. Cells were washed and sorted on a BD Aria sorter (BD Biosciences) in a Biosafety Level 3 laboratory. Non-naive B cell populations were gated by CD45+CD19+IgD− or CD27+, and non-naive T cell populations were gated by CD45+CD3+CCR7low or CD95+.

Single-cell RNA sequencing

PBMC samples were partitioned into single-cell gel bead in emulsion (GEM) mixed together with the reverse transcription (RT) mix using 5′ Chromium Single Cell Immune Profiling Next GEM (version 1.1 Chemistry) (10x Genomics), as previously described25. The RT step was conducted in the Veriti Thermo Cycler (Thermo Fisher Scientific). Single-cell gene expression, cell surface protein and TCR and BCR libraries were prepared as instructed by 10x Genomics user guides. All libraries were quality controlled using Bioanalyzer (Agilent) and quantified using Qubit Fluorometric (Thermo Fisher Scientific). 10x Genomics 5′ single-cell gene expression, cell surface protein tag and TCR and BCR libraries were pooled and sequenced on an Illumina NovaSeq platform (Illumina) using the sequencing parameters recommended by the 10x Genomics 5′ version 1.1 user guide.

Bulk RNA sequencing and single-cell sample demultiplexing

For each sample, 100,000–500,000 cells were processed in TRIzol using the miRNAeasy Micro Kit (Qiagen), and standard RNA sequencing libraries were generated using Illumina TruSeq library preparation kits. The results of bulk RNA sequencing were used for demultiplex of CITE-seq samples by generating single-nucleotide polymorphism (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 protocol74. The software package demuxlet was used to then match single-cell gene expression data to each donor and identify empty droplets and doublets. Because multiple samples from different time points 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.

CITE-seq quantification and statistical analysis

Single-cell data processing and clustering

Single-cell data processing, CITE-seq protein data denoise and clustering were performed as described previously25. Specifically, Cell Ranger (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 (version 3.1.0) running in R version 3.6.1. After filtering to single cell based on demuxlet output, we further demultiplexed the time points using the hashtag antibody staining. We removed cells with less than 250 or more than 4,000 detected genes, more than 20% mitochondrial reads, cell surface protein tag greater than 200,000 and hashtag antibody counts greater than 50,000. The protein data were normalized and denoised using the DSB method75. The following parameters were used in the DSB normalization function: define.pseudocount = TRUE, pseudocount.use = 10, denoise_counts = TRUE and 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 graph followed by k-nearest neighbors clustering were built using the FindNeighbors and FindClusters functions in Seurat (version 3.1.0), respectively. Major cell clusters were then manually annotated using the surface protein together with gene expression. Major cell clusters identified based on protein expression profile and shown in Fig. 4a included: B_Mem: Memory B cells; CD4_Mem: Memory CD4 T cells; CD8_Mem: Memory CD8 T cells; CD4_isobinding: isotype antibodies binding CD4 T cells; cDC: Conventional dendritic cells; cKit+CD3- activated: cKit high cells with enrichment of activated T cell signatures but lacking surface CD3, CD4 and CD8 expression; dim: low quality, cell subsets with high mitochondria/ribosome genes and most surface markers lowly expressed; DNT: Double negative T cells; HSC: Hematopoietic stem or progenitor cells; MAIT: Mucosal-associated invariant T cells; Mono_Classical/Intermediate/NonClassical: Classical/Intermediate/NonClassical Monocytes; and NK_CD16hi/NK_CD56hiCD16lo: CD16 highly expressed/CD56 highly and CD16 lowly expressed NK cells.

Label transfer for cell annotations

To compare the cell population frequencies directly with patients with aCOVID-19 and to avoid potential annotation batch effect, the previously published aCOVID-19 dataset25 was projected onto CITE-seq data—query from this experiment in Seurat (version 3.1.0) using the FindTransferAnchors function. log-normalization and the first 30 principal components were used for the integration. Cell annotations were then predicted using the TransferData function, and the predicated labels were added to the metadata as predicated.id column.

Pseudo-bulk differential expression and GSEA

Pseudo-bulk gene differential expression analysis and GSEA were performed as described previously25. In brief, all unsorted cells in a given sample were computationally ‘pooled’ according to their cluster assignment by summing all reads for a given gene. Pseudo-bulk 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 five cells and fewer than 40,000 unique molecular identifier counts detected after pooling. Lowly expressed genes were removed for each cell type individually using the filterByExpr function from edgeR76. Differentially expressed genes were identified using the limma voom77 workflow, which models the log of the counts per million (CPM) 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. Gene set lists used for enrichment assessment (including GO BP, KEGG, Reactome, the Molecular Signatures Database’s Hallmark collection, Blood Transcriptomic Modules and a few published datasets) were the same as described in Liu et al.25. P values were adjusted using the Benjamini–Hochberg method for the whole gene set list. 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.

Models used for differential expression: patients with MIS-C and patients with pCOVID-19 versus pHCs

Using the pseudo-bulk limma voom workflow as described in ‘Pseudo-bulk differential expression and GSEA’, differentially expressed genes between patient samples (with admission days <41) and pHCs 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 patients and pHCs.

Models used for differential expression: patients with MIS-C versus patients with pCOVID-19

Similarly, differentially expressed genes between MIS-C samples (with admission days <41) and pCOVID-19 samples were identified with a model with the following formula in R: ~0 + mis-c_vs_pediatric_covid + days_since_admission + age, where mis-c_vs_pediatric_covid is a factor variable with two levels, and 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 pCOVID-19.

Models used for differential expression: time effect on gene expression in patients with MIS-C and patients with pCOVID-19

Differentially expressed genes of MIS-C samples and pCOVID-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 pCOVID-19, respectively.

Gene set module scores calculation

Selected module scores (gene set signature score) representing enriched pathway activities were calculated for each sample as reported previously25. Specifically, LE genes identified by GSEA from the MIS-C versus pCOVID-19 model above were used to enhance signal-to-noise ratio and highlight mainly the differences between MIS-C versus pCOVID-19. The pseudo-bulk gene counts were normalized with the varianceStabilizingTransformation function from DESeq2 (ref. 78) for the score calculation. The scores were generated using the gene set variation analysis (GSVA) method from the GSVA R package.

TCR and BCR data processing

Cell Ranger (10x Genomics) version 3.1.0 was used to assemble V(D)J contigs (https://support.10xgenomics.com/single-cell-vdj/software/pipelines/latest/algorithms/annotation). For TCR data, the V(D)J assignment and clonotype were from 10x Cell Ranger output of the filtered_contig_annotations.csv file. For BCR data, V(D)J sequencing contigs from 10x Cell Ranger output were processed using Immcantation version 3.0.0 toolbox (https://immcantation.readthedocs.io/en/latest/index.html). IgBLAST and IMGT germline sequence databases and Change-O package79 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.80 using the TIgGER R package81 and the ShazaM R package79. The TCR and BCR sequence data, contig assignments and estimated BCR mutation frequencies were combined, respectively, using scRepertoire R package and integrated with the single-cell RNA sequencing Seurat object in the metadata.

CITE-seq data visualization

For heat maps showing pseudo-bulk gene expression profiles, the log of counts per million (CPM) for each sample and gene for a given cell type was calculated by pooling cells as described in ‘Pseudo-bulk differential expression and GSEA’. Library size normalization was performed without additional scaling factors, and heat maps were scaled to Z-score among samples for each gene. ComplexHeatmap82 and pheatmap were used for plotting heat maps using R. The ggplot2 and ggpubr R packages were used for box, bubble and scatter plot visualization.

Validation of gene set enrichments in external single-cell RNA sequening data from ref. 13

Single-cell data from the cohort of Ramaswamy et al.13 were downloaded from fastgenomics (the Ramaswamy2021_MIS-C_10x_PBMC dataset). Using the pre-annotated cell clusters from the original publication, single-cell gene expression data were pooled into pseudo-bulk libraries, and differential expression and GSEAs of patients with MIS-C versus pHCs were done as described in ‘Pseudo-bulk differential expression and GSEA’ and ‘Models used for differential expression: patients with MIS-C and patients with pCOVID-19 patients versus pHCs’. Age was included in the model as a covariate.

HuProt auto-antibody analysis

Auto-antibody analysis was performed using HuProt version 4.0 human protein microarrays and processed by CDI Laboratories. IgG profiling was performed for 15 serum samples from five children with pCOVID-19 and ten children with MIS-C, of whom four had received IVIG. In brief, the arrays were blocked and probed with the samples at 1:1,000 dilution and incubated at room temperature for 1 h. Then, the arrays were washed and probed with Alexa Fluor 647 anti-human IgG (Fc) for signal detection as previously described. Using CDI software, quantile normalization of the raw signal intensities (F635 median for IgG and 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 as IGHG1, IGHG3 and so on) alongside the controls (such as rhodamine+IgG64, anti-human IgG and GST 10 ng µl−1). The quantile-normalized IgG binding intensities of the remaining 23,040 protein targets were then visualized using 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 10,000,000; the fold change of average signal intensity was greater than 4 between the two groups; and the FDR was less than 0.5.

Multiplex particle-based anti-cytokine auto-antibody screening assay and functional evaluation

Plasma samples were screened for auto-antibodies against IFN-α, IFN-β, IFN-ω and IFN-γ in a multiplex particle-based assay83, in which differentially fluorescent magnetic beads were covalently coupled to recombinant human proteins (2.5 μg per reaction). Beads were combined and incubated for 30 min with diluted plasma samples (1:100 dilution). Beads were then washed and incubated with PE-labeled goat anti-human IgG (1 µg ml−1) for an additional 30 min. Beads were washed again, resuspended in assay buffer and analyzed on a BioPlex X200 instrument. Plasma samples with a fluorescence intensity greater than 1,500 were tested for blocking activity. The blocking activity of auto-antibodies was determined by assessing STAT1 phosphorylation in healthy control cells after 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−1 of IFN-α, IFN-β, or IFN-ω or 400 units per milliliter of IFN-γ for 15 min at 37 °C. Cells were fixed, permeabilized and stained for intranuclear pSTAT1 (Y701). Cells were acquired on a BD LSRFortessa cytometer, gated on CD14+ monocytes, and analyzed with FlowJo software.

Reporting Summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.