SARS-CoV-2 induced immune perturbations in infants vary with disease severity and differ from adults’ responses

During the first year of the COVID-19 pandemic, SARS-CoV-2-related morbidity and mortality disproportionately affected adults. However, as the pandemic evolved and vaccines became available first for adults, an increasing number of symptomatic infections were observed in children and infants, in part related to the emergence of new variants. This shift challenged the initial belief that children were mostly resistant to SARS-CoV-2 infections. Indeed, recent data from the CDC indicate that infants <6 months of age with SARS-CoV-2 infection had high rates of hospitalization during the later phase of the pandemic¹.

Studies in adults with severe COVID-19 have shown increased plasma concentrations of pro-inflammatory cytokines, including IL-1β and the interferon (IFN)-induced chemokine CXCL10 (IP-10)², together with significant alterations in transcriptional profiles and leukocyte composition such as lymphopenia and increased neutrophil counts^3,4. In adults with mild disease, a marked early interferon-stimulated gene (ISG) signature was identified, though this signature was absent in more severe cases⁵. Conversely, other studies reported that interferons (i.e., Type I and Type III) and ISG signatures were induced only at late stages of severe infection, suggesting that IFNs might play a detrimental rather than protective role during COVID-19^6,7.

Comparing systemic immune responses to SARS-CoV-2 infection between children and adults revealed increased numbers of naive T and B lymphocytes and reduced numbers of NK cells in children, whereas adults exhibited significantly increased cytotoxic and ISG^high T cell subsets⁸. In addition, school-age children with COVID-19 showed enhanced mucosal expression of innate immune pathways, including the NLRP3 inflammasome, and increased production of mucosal cytokines such as IFN-α2, IFN-γ, IP-10, IL-8, and IL-1β⁹, compared to adults. Another study highlighted enhanced local IFN responses in the respiratory tract of children, in contrast to a more significant systemic IFN response in adults¹⁰. Higher baseline expression of pattern recognition receptors MDA5 and RIG-I was observed in upper airway epithelial cells, macrophages, and dendritic cells of older children (median age ~9 years) compared with adults, suggesting that children might be primed to produce IFN in response to viral triggers¹⁰. In summary, while most studies suggest that school-age children exhibit a more effective local innate immune response to SARS-CoV-2 than adults, it remains unclear whether they also mount stronger systemic IFN responses.

Despite significant advances in understanding the impact of COVID-19 on adults and school age children, we still know surprisingly little about how SARS-CoV-2 affects the immune system of infants. A recent study¹¹ comprehensively analyzed the immune responses induced by SARS-CoV-2 infection (including both Omicron and non-Omicron variants) in infants and young children with mild disease (median age ~9 months), who were not hospitalized and did not require high level of care. That study showed the emergence of spike antibodies 4–5 days post infection, which persisted for up to 300 days. Infants and young children showed a robust Th17 response in the nasal mucosa, but no significant detection of pro-inflammatory molecules/cytokines in plasma. However, that cohort did not include infants hospitalized with severe disease, leaving a significant gap in our understanding of how developing immune system handles severe SARS-CoV-2 infection, particularly during the critical first weeks of life. To address this question, we recruited 26 infants hospitalized with SARS-CoV-2 infection during the first weeks of life (median age ~1.6 months) and 14 healthy matched controls. We analyzed their blood immune cell transcriptome at the single-cell level and evaluated serum cytokine profiles, viral loads, and antibody responses. We subsequently compared our findings with published data from both infected adults and infants with mild disease not requiring hospitalization. Unique features of the immune responses to SARS-CoV-2 infection through a wide range of clinical severity during the first year of life are reported.

Clinical categorization of infected infants and changes in PBMC compositions

We collected blood samples and processed serum and peripheral blood mononuclear cells (PBMCs) from 26 infants (median, IQR age: 1.63 [0.93–7.59] months) hospitalized with SARS-CoV-2 infection (pediatric COVID; pCoV) and 14 matched pediatric healthy controls (pHC; median IQR age: 2.01 [1.86–4.28] months) (Fig. 1a). Patient demographics, clinical and laboratory data and treatment are summarized in Supplementary Data 1a, b. We categorized pCoV patients into three groups (G1, G2 and G3; Fig. 1b) based on (i) the time (days) since exposure (DSE) to a known COVID-19 case, (ii) disease severity, and (iii) semi-quantitative viral loads (VL) measured within a median of 17 [11–36] hours of hospitalization. Infants in G1 (n = 10) were diagnosed with COVID-19 as part of the diagnostic work-up for infants prompting hospitalization mainly for fever and/or respiratory symptoms. Infants in G1 had a median of 21 days since they were exposed to an index case and had lower SARS-CoV-2 viral loads (subacute group). Infants in G2 (n = 11) had moderate disease and high SARS-CoV-2 loads, but did not require significant medical interventions such as supplemental oxygen. One infant in G2 received intravenous fluids. G3 (n = 5) infants also had high SARS-CoV-2 viral loads but required a higher level of care, such as oxygen supplementation, intensive care unit (ICU) admission or mechanical ventilation (Supplementary Data 1a–d).

a Study design. b Viral load quantification across the three pCoV groups (Patient groups; G1, G2 and G3). Days since exposure (DSE), cycle threshold (Ct) values and disease severity are shown for each patient group (see Supplementary Data 1 for more details). c UMAP plot representing 203,402 cells colored by annotated immune subsets. Top 100 differentially expressed genes in each subset can be found in Supplementary Data 2. d Immune subset annotation. Dot plot representing expression values of selected genes (x-axis) across each cluster (y-axis). Dot size represents the percentage of cells expressing the marker of interest. Color intensity indicates the mean expression within expressing cells. e UMAP plot representing Interferon-stimulated genes (ISG) score (gene list available in Supplementary Data 3). f UMAP plots representing cell compositional changes (cell density) across healthy controls (pHC) and patient groups (G1, G2 and G3). g Bar plot showing the cell frequency of each cluster across the 40 individuals.

Full size image

To assess the immune response to SARS-CoV-2 in these infants and to determine how these responses correlated with their varying clinical phenotypes, we performed: (i) single cell RNA-seq (scRNA-seq) of PBMCs, (ii) serum cytokine concentration assessment using the 92-Olink platform with a focus on inflammatory pathways, (iii) quantification of antibodies against coronaviruses (anti-Spike and anti-RBD) and anti-IFN activity, and (iv) semi-quantitative viral load measurement in nasopharyngeal (NP) swabs (Fig. 1a, b and Supplementary Fig. 1a).

For scRNA-seq data analysis, raw data from pHC PBMCs (38% of the pool) and pCoV PBMCs (62% of the pool) were integrated (Fig. 1a). Hybrid transcriptomes (multiplets) were identified using scrublet¹² and excluded from the rest of the analysis. After filtering steps, pHC and pCoV samples yielded, respectively, a mean (±S.D.) of 6559 (±2016) cells and 5981 (±1807) cells per individual, and a mean of 1411 and 1603 genes per cell (Supplementary Fig. 1b, c, Supplementary Data 2). scRNA-seq profiles that passed quality control (n = 203,402 cells) were then corrected for batch effect using BBKNN¹³ (Supplementary Fig. 1d, e). Unsupervised clustering of the corrected data, followed by a two-dimensional uniform manifold approximation and projection (UMAP), yielded 25 clusters independently of 10x run (Supplementary Fig. 1e) and subject (Supplementary Fig. 1f) batch effects. Principle Variance Component Analysis which quantifies the contribution of different factors (e.g., biological groups, sex, age, technical batches, etc) into the total variance, showed that “patient groups” was the major contributor (Supplementary Fig. 1g). Using marker genes, clusters were first mapped to major immune cell types, including CD4⁺ and CD8⁺ T cells, B cells, NK cells, CD14⁺ and CD16⁺ monocytes, conventional dendritic cells (cDCs), plasmacytoid dendritic cells (pDCs), plasma cells (PCs), hematopoietic stem cells (HSCs) and Erythroblasts/Erythrocytes (Eryth) (Fig. 1c, d, Supplementary Fig. 1h, i). To analyze the cellular distribution of interferon stimulated genes (ISGs) within different immune subsets, we calculated an ISG score based on expression levels of genes in interferon modules¹⁴ (Supplementary Data 3). This calculation revealed increased ISG scores in most cell types, including monocytes, T cells and B cells (Fig. 1e). The comparison of PBMCs from acutely infected infants (combined G2 and G3) vs. healthy controls further confirmed a broad ISG upregulation across all cell subsets (Supplementary Data 5a). Cell compositional analyses in pHC, G1, G2 and G3 pCoV groups revealed comparable distribution of pHC and G1 groups, contrasting with striking differences observed in G2 and G3 groups (Fig. 1f, g). Overall, these high-level analyses revealed significant alterations in the PBMC composition of infants upon SARS-CoV-2 infection and a robust ISG signature in most immune cell types in response to the infection. To refine these alterations, we further subclustered and analyzed the different immune cell subsets.

Increased numbers of ISG^hi IL-1B⁺ CD14⁺ monocytes in infants with mild, moderate and severe disease

Previous studies have reported alterations in the CD14⁺ compartment of adults with severe COVID-19, including decreased HLA-DR surface expression and the accumulation of immature circulating myeloid cells, suggesting dysregulated myelopoiesis/innate immune response^15,16. To investigate this in our infant cohort, we further clustered CD14⁺ monocytes (n = 16,906) into four subclusters (SCs) (Fig. 2a). CD14⁺ mo-SC0 (n = 5495) expressed antigen presentation-related transcripts (e.g., HLA-DRA and CD74) together with small and large ribosomal subunits (e.g., RPS3A and RPL3), while CD14⁺ mo-SC2 (n = 4590) expressed high levels of alarmins S100A8 and S100A9. CD14⁺ mo-SC1 (n = 5017) displayed an ISG signature, including ISG15, IFI6, and MX1. CD14⁺ mo-SC3 (n = 1804) was characterized by the upregulation of IL1R2, CD163, FKBP5, as well as ISGs (Fig. 2b, Supplementary Data 6a). CD14⁺ mo-SC1 was significantly increased in G2, while CD14⁺ mo-SC0 and CD14⁺ mo-SC2 were predominantly present in pHC and G1 (Fig. 2c). CD14⁺ mo-SC3 was almost entirely contributed by 3/5 patients with severe disease (G3), two of whom had received systemic steroids (Supplementary Data 1a). In fact, the CD14⁺ mo-SC3 transcriptional profile (Fig. 2b, c) was reminiscent of what was previously described in adults with SARS-CoV-2 infection receiving steroid treatment^17,18,19. The comparison of CD14⁺ monocytes from acutely infected (G2 and G3) infants versus healthy controls revealed a clear upregulation of ISGs (Fig. 2d and Supplementary Data 5b). The ISG^hi subcluster (CD14⁺ mo-SC1) was infrequent in pHC and significantly increased in 2/10 of G1, 11/11 of G2, and 2/5 of G3 patients (Fig. 2e). As both inflammation and type I IFN pathways were previously found activated in monocytes from adults with SARS-CoV-2 infection⁵, we analyzed the co-expression of inflammatory transcripts and ISGs at the single cell level in these infants. Approximately, 49% of SC1 monocytes co-expressed IL1B and ISG15 (Supplementary Fig. 2a, b).

a UMAP plots representing CD14⁺ monocyte subclusters (SCs; n = 4, left panel), patient groups (G1, G2, G3) or pHC. b Density plot, generated by the Nebulosa⁵¹ package, showing the expression of selected marker genes. c Boxplot comparing the proportion of each CD14⁺ monocyteSC (n = 4) across the individuals (n = 40). pCoV patient groups (G1, G2 and G3) are shown in purple, and pHC in green. P values were calculated using a two-sided t-test comparing the mean (G1, G2 or G3) with mean (pHC). ns: non-significant. The upper and lower bounds represent the 75% and 25% percentiles, respectively. d Heat map representing scaled expression values of the most significant genes in CD14⁺ monocytes from acutely infected (G2 and G3) infants, compared with healthy controls. e Bar plot showing the cell frequency of each CD16⁺ monocyte SC across the 40 individuals, as categorized by patient groups (G1, G2 and G3; in purple) and pHC in green. f As in (a) in CD16⁺ monocyte SCs (n = 3). g As in (c) in CD16⁺ monocytes. h As in (d) in CD16⁺ monocytes. i As in (a) in cDC SCs (n = 7). j Violin plot showing median expression values for selected marker genes across the cDC SCs. k As in (c) in cDC SCs.

Full size image

Thus, infants with moderate and severe COVID-19 were separated in two groups: the first included 2/10 G1, all G2 and 2/5 G3 infants and was characterized by increased numbers of CD14⁺ monocytes co-expressing inflammatory molecules and ISGs with low HLA gene expression. The second group consisted of 3/5 G3 infants whose monocytes expressed ISGs and IL1R2 along with CD163, likely reflecting the influence of steroid treatment.

Infants with SARS-CoV-2 infection displayed increased interferon expression but normal frequency of CD16⁺ monocytes

Contraction of CD16⁺ monocytes has been previously described as a feature of severe COVID-19 in adults¹⁶. Subclustering of CD16⁺ monocytes in infants (n = 5987) generated three SCs (Fig. 2f). CD16⁺ mo-SC0 (n = 3009) expressed an ISG signature (e.g., ISG15, IFI6, and MX1). CD16⁺ mo-SC1 (n = 2217) expressed ribosomal-associated genes RPL5, RPL10A, or RPL3 and CD16⁺ mo-SC2 (n = 761) was characterized by the expression of inflammatory molecules (e.g., IL1B and CCL3) (Supplementary Fig. 2c, Supplementary Data 6b). CD16⁺ mo-SC0 (n = 3009) was significantly increased in G2 and G3; SC1 was enriched in pHC and G1, while CD16⁺ mo-SC2 showed a mild increase in G2/G3 (Fig. 2g). In contrast to the CD14⁺ population, only about 13% of CD16⁺ SC0 monocytes co-expressed IL1B and ISG15 (Supplementary Fig. 2d). The comparison of CD16⁺ monocytes from acutely infected (G2 and G3) infants with those from healthy controls revealed an enrichment in ISGs (Fig. 2h and Supplementary Data 5c). The ISG^hi subcluster SC0 was rare in pHC and significantly increased in 9/11 of G2 and 4/5 of G3 infants (Supplementary Fig. 2e). Unlike severe COVID-19 in adults, which was characterized by CD16⁺ monocyte depletion²⁰, infants with moderate to severe disease (G2/G3) exhibit normal frequencies of monocytes due to the increased numbers of ISG^hi CD16⁺ monocytes.

To explore how our findings compare with previous studies, we integrated our scRNA-seq data with a recent study by Wimmers et al.¹¹, which analyzed PBMCs from older infants with mild disease (average 9 months of age, not hospitalized and untreated) during both Omicron and non-Omicron variants of SARS-CoV-2 infection. This cohort included: (i) nine samples collected before (pre, n = 9), during (acute; n = 10) and after infection (convalescent or conv; n = 9) with non-Omicron variants, (ii) eight samples during acute Omicron infection (acute omicron) and (iii) seven matched controls. Our integrative analysis of CD14⁺ and CD16⁺ monocytes from 83 infants in total confirmed the presence of an ISG signature in infants with moderate to severe disease (10/10 G2 and 5/5 G3 patients in our cohort), which contrasted with only 5/10 acute non-Omicron and 1/10 acute Omicron patients in the Wimmers cohort (Supplementary Fig. 2f–i). These differences are likely attributable to variations in the time from exposure to the virus, the age of the infants, and, most importantly, the clinical disease severity between the two cohorts.

SARS-CoV-2 infected infants showed reduced frequency of conventional dendritic cells (DCs)

Similar to classical monocytes, DCs from SARS-CoV-2 infected adults exhibit lower levels of HLA-DR compared to healthy controls²¹. Subclustering conventional DCs (cDCs) in infants (n = 440) generated seven SCs (Fig. 2i and Supplementary Fig. 3a), which were further classified into (i) cDC1 (SC5; expressing CLEC9A and XCR1), (ii) cDC2 (SC0/SC3; expressing CD1C and CLEC10A), (iii) monocyte-derived cDCs (Mo-DCs; SC1/SC4; expressing CD14 and S100s) and (iv) AXL⁺ DCs (SC6; expressing AXL and DAB2) (Supplementary Data 6c). The AXL⁺ DCs (SC6) expressed the highest levels of ISGs (Fig. 2j and Supplementary Fig. 3b). G2/G3 groups showed a marked reduced numbers of cDC1, cDC2 and Mo-DCs, while two G2 patients showed an increase in AXL⁺ ISG^hi SC (SC6) (Fig. 2k and Supplementary Fig. 3c). The comparison of cDCs from acutely infected (G2 and G3) infants with those from healthy controls revealed a pronounced upregulation of ISGs (Supplementary Data 5e). Overall, our data showed that profiles of G2 and G3 infants were characterized by a decrease in cDC1, cDC2, and Mo-DCs.

The frequency of pDCs in peripheral blood was reduced in adults with COVID-19²², possibly due to pDC apoptosis¹⁷. Subclustering of pDCs in infants (n = 386) generated two SCs (Supplementary Fig. 3d). pDC-SC0 (n = 226) was significantly reduced in G2/G3 patients, while pDC-SC1 (n = 160) was significantly increased in these infants and characterized by the upregulation of ISGs (Supplementary Fig. 3e, f). The comparison of pDCs from acutely infected (G2 and G3) infants with those from healthy controls revealed a pronounced upregulation of HLA related molecules and ISGs (Supplementary Data 5f). The depletion of ISG^low pDCs is consistent with reports in adults with COVID-19; however, a subset of ISG^hi circulating pDCs was increased in G2/G3 infants.

Circulating megakaryocytes (MGKs) have previously been linked to COVID-19 outcomes in adults²³. Two of the infant’s SCs (in C22, Supplementary Fig. 1h, i) were annotated as megakaryocyte (MGK; n = 408; expression of PPBP, CLU and PF4) and hematopoietic stem cells (HSCs; n = 305; expression of CD34 and ITM2C) (Supplementary Fig. 3g, h, Supplementary Data 6e). MGKs (n = 408 cells) were significantly increased while HSC (n = 305 cells) were significantly decreased in G1 and G3 groups (Supplementary Fig. 3i).

ISG^hi naive CD4⁺ T cells are increased in SARS-CoV-2 infected infants

Activation and proliferation of both CD4⁺ and CD8⁺ T cells have been associated with COVID-19 disease severity^4,17, though no specific T cell signatures have been identified as primary correlates in this context. In this infant study, subclustering CD4⁺ T cells (n = 100,482) yielded four SCs (Fig. 3a). CD4⁺T-SC0 (n = 43,567), expressing ribosomal protein-associated genes such as RPL5, RPL10A and RPL3, was significantly decreased in G2 and G3 infants. CD4⁺T-SC1 (n = 22,847), which exhibited an ISG signature (e.g., IFI44L, ISG15 and XAF1), was increased in G2 and G3 infants. CD4⁺T-SC2 (n = 14,383), characterized by expression of JUNB/JUND, CD69 and CXCR4 did not differ between infants with COVID-19 and pHC. Finally, CD4⁺T-SC3 (n = 14,383), expressing the memory marker S100A4, was increased in all three pCoV groups (Fig. 3b, Supplementary Data 6f). The ISG^hi CD4⁺T-SC1 was present in 2/10 G1, 8/11 G2, and 5/5 G3 infants. Notably, more than 90% of CD4⁺ T cells from a fraction of G2 and G3 patients (e.g., pCoV11 and pCoV26) switched to an ISG^hi state (Fig. 3c). Overall comparison of CD4⁺ T cells from acutely infected (G2 and G3) infants with those from healthy controls revealed an enrichment in ISGs (Fig. 3d, Supplementary Data 5g). CD4⁺T-SC0 and CD4⁺T-SC1 exhibited naive (CCR7, LEF1 and SELL) markers, while CD4⁺T-SC3 displayed memory markers (Fig. 3b). Our analysis also revealed the expression of the transcription factor SOX4 within the naive compartment, which is consistent with an immature/naive phenotype (Supplementary Fig. 4a). While the majority of ISG^hi T cells originated from infected patients, we also detected these cells in pHC, indicating constitutive expression of ISGs in young infants (Supplementary Fig. 4b, c).

a UMAP plots representing CD4⁺ T cell subclusters (SCs; n = 4) and patient groups (G1, G2, G3 or pHC). b Heat map representing scaled expression values of selected marker genes. c Boxplot comparing the proportion of each CD4⁺ T cell SC (n = 4) across the individuals (n = 40). pCoV patient groups (G1, G2 and G3) are shown in purple, and pHC in green. P values were calculated using a two-sided t-test comparing the mean (G1, G2 or G3) with mean (pHC). ns: non-significant. The upper and lower bounds represent the 75% and 25% percentiles, respectively. d Heat map representing scaled expression values of the most significant genes in CD4⁺ T cells from acutely infected (G2 and G3) infants, compared with healthy controls. e UMAP plot representing memory CD4⁺ T cell subsets. f Density plot, generated by the Nebulosa⁵¹ package, showing the expression of selected marker genes. g As in (c) in memory CD4⁺ T cell SCs. h As in (f) in Tregs. i As in (g) in Tregs.

Full size image

Further analysis of CD4⁺T-SC3, resulted in the separation of memory (Tmem, six subsets) and regulatory T cells (Tregs; two subsets). Tmem included seven subsets: (a) central memory (CM), (b) SOX4⁺ CM, (c) ISG^hi, (d) GZMK⁺ Tmem (Gzk), (e) KLRB1⁺ Tmem, (f) TH2 (PTGDR2, GATA3-AS1) and (g) TH17 (CTSH) (Fig. 3e, f). Infected infants exhibited increased numbers of ISG^hi Tmem, along with reduced numbers of GzK and TH17 cells (Fig. 3g). Treg were separated into naive and memory (Fig. 3e, h). G3 infants exhibited lower naive Tregs and higher memory Treg frequencies (Fig. 3i). Within the memory compartment, our analysis revealed the presence of cells expressing both GMZK and Temra-associated transcripts such as PRF1 and CCL5, an unusual finding in the adult CD4⁺ T cell memory compartment (Fig. 3f).

Analysis of granzyme (GZM) family expression revealed a broad upregulation of GZMM, contrasting with an absence of GZMB and restriction of GZMK to CD4⁺T-SC3.

The integration analysis with the Wimmers et al. dataset further confirmed an augmented ISG signature in all G2 and G3 individuals in our cohort, contrasting with 4/10 acute non-Omicron and 1/10 acute Omicron individuals in the Wimmers cohort (Supplementary Fig. 5a, b). Overall, these data indicate that the blood CD4⁺ T cells of COVID-19 infected infants fall into two main compartments: (i) naive, as expected the predominant population, that had switched to an ISG^hi state, and (ii) memory, which encompassed a mixture/combination of central and effector memory CD4⁺ T cells, T helper subsets and Tregs that was altered in SARS-CoV-2 infected infants.

ISG^hi CD8⁺ T cells are increased in infants with COVID-19

COVID-19 associated lymphopenia³ is reported to preferentially impact CD8⁺ T cells⁴, possibly due to T cell exhaustion²⁴. However, the overall CD8⁺ T cell compartment (n = 34,366), was not decreased in infants with COVID-19 and was distributed into five SCs (Fig. 4a). CD8⁺ T-SC0 (naive, n = 23,768) expressed ribosomal-associated genes (e.g., RPL34 and RPL39), and was decreased in G2 and G3 infants; CD8⁺T-SC1 (ISG^hi; n = 3483) exhibited an ISG signature (e.g., IFI44L, ISG15 and XAF1) and was increased in G2/G3 infants. CD8⁺ T-SC2 (n = 2967) displayed similar frequencies in infected infants and pHC and expressed cytotoxic genes, including PRF1, GZMA and GNLY, inflammatory chemokines (CCL4 and CCL5), S100A4 and CX3CR1 consistent with a TEMRA phenotype. CD8⁺ T-SC3 (GzK; n = 2201) expressed GZMK and CCL5 and was mostly present in G1, but reduced in G2/G3 infants. Finally, CD8⁺T-SC4 (Prolif.; n = 1947), which expressed proliferation markers such as MKI67 and HMGBs, was increased in G2 and G3 infants (Fig. 4b, c, Supplementary Data 6g). The CD8⁺ T ISG^hi SC was rare in pHC but present in 2/10 G1, 8/11 G2, and 3/5 G3 infants. Notably, more than 50% of CD8⁺ T cells in some G2 (pCoV13, pCoV24, pCoV26) and G3 (pCoV2 and pCoV11) infants demonstrated an ISG^hi state (Fig. 4d). A direct comparison of CD8⁺ T cells from acutely infected (G2 and G3) infants with those from healthy controls revealed an upregulation of ISGs (Supplementary Data 5h). The proliferative SC4 included a subset of cells expressing TEMRA markers, such as CX3CR1, FCGR3A and GZMB (Fig. 4e, f). This SC represented 15–20% of CD8⁺ T cells in 6/16 of the G2/G3 infants (Fig. 4d). Further clustering of proliferative CD8⁺ T cells revealed three subsets, including a naive (CCR7, SOX4), a granzyme K (GZMK) and cytotoxic (PRF1 and GZMB) subsets (Fig. 4e), the latter being significantly increased in G2 (Fig. 4g). Comparison with published data¹¹ on infants with mild disease further confirmed the presence of a robust ISG signature in all G2 and G3 patients in our cohort, contrasting with 4/10 acute non-Omicron and 1/10 acute Omicron individuals in the Wimmers cohort (Supplementary Fig. 5c, d). In summary, infants with more severe disease exhibited increased numbers of ISG^hi CD8⁺ T cells and induction of CX3CR1⁺ proliferating CD8⁺ T cells with a cytotoxic program, a finding that contrasts with previous reports in adults²⁵.

a UMAP plots representing CD8⁺ T cell subclusters (SCs, n = 5) and patient groups (G1, G2, G3 or pHC). b Heat map representing scaled expression values of selected ISGs and marker genes across CD8⁺ T cell SCs. c Boxplot comparing the proportion of each CD8⁺ T cell SC (n = 5) across the individuals (n = 40). pCoV patient groups (G1, G2 and G3) are shown in purple, and pHC in green. P values were calculated using a two-sided t-test comparing the mean (G1, G2 or G3) with mean (pHC). ns: non-significant. The upper and lower bounds represent the 75% and 25% percentiles, respectively. d Bar plot showing the cell frequency of each CD8⁺ T cell SC across the 40 individuals, as categorized by patient groups (G1, G2 and G3; in purple) and pHC in green. e UMAP plot representing proliferative CD8⁺ T cell SCs. f Density plot, generated by the Nebulosa⁵¹ package, showing the expression of selected marker genes. g As in (c) in proliferative CD8⁺ T cell SCs. h As in (a) in NK cells. i As in (b) in NK cells. j As in (c) in NK cells. k As in (d) in NK cells.

Full size image

ISG^hi NK cells are increased in infants with severe SARS-CoV-2 infection

A recent study reported an increase of “adaptive” NK cells in the blood of adults with severe COVID-19²⁶. However, little is known about the role of these cells in infants. Herein, NK cell (n = 13,403) subclustering yielded four SCs (Fig. 4h). NK-SC0 (n = 5013), reduced in G2 and G3 infants, expressed ribosomal-associated genes. NK-SC1 (n = 3393), reduced in G2 infants, expressed CD160 and KLRB1. NK-SC2 (n = 2721), increased in G2/G3 infants, exhibited an ISG signature. NK-SC3 (n = 2276) displayed similar frequencies in pHC and infected infants and upregulated transcripts such as XCL1, XCL2 and GZMK (Fig. 4i, j). While very rare in pHC, the ISG^hi SC2 was present in 2/10 G1, 11/11 G2, and 5/5 G3 infants. Notably, more than 70% of NK cells from 9/11 G2 and all 5 G3 infants switched to an ISG^hi state (Fig. 4k, Supplementary Data 6h). The direct comparison of NK cells from acutely infected (G2 and G3) infants with those from healthy controls confirmed the enrichment in ISGs (Supplementary Data 5i).

The comparison with previously reported cohorts confirmed the presence of a strong ISG signature in all G2 and G3 patients in our cohort, contrasting with 5/10 acute non-Omicron and 1/10 acute Omicron patients in the Wimmers cohort (Supplementary Fig. 5e, f).

SARS-CoV-2 infection in infants is associated with increased ISG^hi transitional and naive B cells

In adults with COVID-19, activated B cells and plasmablasts are expanded^4,27, and an increased frequency of extrafollicular (Tbet⁺CD11c⁺CXCR5^neg) B cells correlates with disease severity²⁸. In our study, further clustering of infant’s B cells (n = 30,119) generated five SCs (Fig. 5a). Differential marker expression distinguished three compartments: (i) naive B cells (NBC; IGHD, CCR7 and SELL), (ii) Transitional B cells (TrBC; MME (encoding CD10), CD9, and CD24) and (iii) Memory B cells (MBCs; CD27, IGHA1, IGHG1 and S100A4). An ISG signature (e.g., IFI44L, ISG15 and IRF7) was detected in both the NBC and TrBC compartments (Fig. 5b). In G2/G3 infants, B-SC1 (n = 5228) and B-SC2 (n = 4985) were increased, while B-SC0 (n = 13,820), B-SC3 (n = 4857) and B-SC4 (n = 1229) were decreased (Fig. 5c). Thus, NBCs included B-SC0 (ISG^low; increased in pHC/G1) and B-SC2 (ISG^hi; increased G2/G3); TrBCs encompassed B-SC3 (ISG^low; increased in pHC/G1) and B-SC1 (ISG^hi; increased in G2/G3). B-SC4 included transcripts characteristic of double negative (DN)2 memory cells (e.g., ITGAX, TFEC and FGR) (Fig. 5b, Supplementary Fig. 4d, Supplementary Data 6i). The ISG^hi naive B-SC2 expressed TLR7, a marker of both extrafollicular activated naive and DN2 cells²⁹ (Fig. 5b). In addition to memory markers (conventional, switched and extrafollicular), B-SC4 also expressed B cell survival and plasma cell differentiation markers such as TNFRSF13B (encoding TACI) and TNFRSF17 (encoding BCMA), which could correspond to cells prone to differentiate into plasma cells. Although rare in pHC, the ISG^hi SC was present in 2/10 G1, 11/11 G2, and 5/5 G3 infants. Overall, >70% of naive (B-SC2) and transitional (B-SC1) B cell compartments from 13/16 of the G2/G3 infants switched to an ISG^hi state (Fig. 5d). The comparison of B cells from acutely infected (G2 and G3) infants with those from healthy controls revealed an enrichment in ISGs (Supplementary Data 5j). While the majority of the ISG^hi B cells originated from infected patients, we also identified an ISG signature in both naive and transitional B cells from healthy control infants (Supplementary Fig. 4e, f).

a UMAP plots representing B cell subclusters (SCs; n = 5) and patient groups (G1, G2, G3 or pHC). b Heat map representing scaled expression values of selected marker genes defining each of the B cell SC (Supplementary Data 6i). c Boxplot comparing the proportion of each B cell-SC (n = 5) across the individuals (n = 40). pCoV patient groups (G1, G2 and G3) are shown in purple, and pHC in green. P values were calculated using a two-sided t-test comparing the mean (G1, G2 or G3) with mean (pHC). ns: non-significant. The upper and lower bounds represent the 75% and 25% percentiles, respectively. d Bar plot showing the cell frequency of each B cell SC across the 40 individuals, as categorized by patient groups (G1, G2 and G3; in purple) and pHC in green. e As in (a) in plasma cells. f UMAP plots representing scaled expression values of selected marker genes. g As in (c) in plasma cells.

Full size image

The comparison between cohorts demonstrated the presence of a strong ISG signature in all G2 and G3 patients in our cohort, contrasting with 10% and 40% of acute Omicron and non-Omicron individuals, respectively, in the Wimmers cohort (Supplementary Fig. 5g, h).

Finally, plasma cells (PC, n = 544) subclustering generated three SCs (Fig. 5e), which were distinguishable based on IGH expression (Fig. 5f). Although the low number of cells is a limitation, cell compositional analysis did not reveal significant differences between COVID-19 infant groups and pHC (Fig. 5g), which contrasts with what was reported in adults with COVID-19³⁰.

Patient stratification using cell frequencies

To integrate cell subcluster analysis data (Supplementary Fig. 6a–d) with clinical disease severity, we performed an unsupervised clustering based on the abundance of detected subclusters (n = 40) across infants with COVID-19 (n = 26) and pHC (n = 14). The study subjects were accordingly clustered into two main sets (Fig. 6a). The first set, which included 10 infants with COVID-19 (3/5 G3, 6/11 G2, and 1/10 G1 infants), was characterized by the expansion of ISG^hi SCs. The second set, which included the remaining 16 COVID-19 infants together with the 14 pHC, displayed ISG^low SCs including naive CD4⁺, CD8⁺ T cells and B cells (Fig. 6b). Interestingly, the two G3 patients who clustered with pHC/G1 (pCoV17 and pCoV21) had received systemic steroids, consistent with the modulatory effect of corticosteroids on type I IFN expression³¹. These data therefore support that clinically severe SARS-CoV-2 infection in infants is associated with a robust ISG response.

a Heatmap representing the abundance of each of the subclusters (n = 40) across the individuals (n = 40). Clinical information are highlighted on top of the heatmap. b Subcluster annotation. Dot plot represents expression values of selected genes (x-axis) across each cluster (y-axis). Dot size represents the percentage of cells expressing the marker of interest. Color intensity indicates the mean expression within expressing cells. c Heatmap representing the most significant (two-sided ANOVA, padj < 0.05) inflammatory-related cytokine (CK) measurements (n = 33 CK) across the individuals (n = 54; G1 n = 7; G2 n = 20, G3 n = 7 and matched HC n = 20; Supplementary Data 7). d Dot plot representing selected CK across the individuals, as categorized per patient groups and healthy controls. *P < 0.05; **P < 0.01; ***P < 0.001: ****P < 0.0001.

Full size image

SARS-CoV-2 infected infants had increased serum concentrations of inflammatory cytokines

To complement our transcriptional profiling studies, we measured serum concentrations of cytokines in 34 children with COVID-19 (7 G1, 20 G2, and 7 G3) and 20 age-matched pHCs using the Olink inflammatory panel (n = 92 analytes). This cohort included 24 of the 26 pCoV infants profiled with scRNAseq (Supplementary Fig. 1a, Supplementary Data 1e). Overall, 72 cytokines were consistently detected in serum samples (Supplementary Fig. 6e); and 33 of them showed significantly different concentrations between disease severity groups (Fig. 6c). IFNγ concentrations were increased in infants with COVID-19 across all disease groups. Cytokines and inflammatory proteins such as IL6, IL8, IL17C, IL18R1 and CXCL10 (IP-10) were particularly increased in G3 infants (Fig. 6d), while soluble CD6 and TNFS12/TWEAK were markedly decreased in this group of severe disease compared with pHCs. Patients in the Wimmers cohort did not show significant increases in plasma inflammatory cytokines possibly due to differences in time since infection and/or clinical disease severity. In summary, SARS-CoV-2 infected infants displayed increased serum concentrations of inflammatory cytokines, especially those infants with severe disease (G3).

Infected infants develop anti-SARS-CoV-2 IgG antibodies but not anti-seasonal coronavirus or anti-IFN responses

To assess the infants’ humoral response to SARS-CoV-2 and seasonal coronaviruses, we measured antibody levels against: (i) SARS-CoV-2 full-length spike, RBD and S2 subunit, (ii) seasonal coronaviruses full-length spike, (iii) S1 subunit of beta-coronaviruses OC43 and HKU1, and (iv) full-length spike of alpha-coronavirus 229E. Samples from 16 infants were obtained at the time of hospitalization, and 13 of those had paired follow-up samples four weeks later. The 13 infants with paired samples, who were also included in the scRNAseq cohort (Supplementary Fig. 1a, Supplementary Data 1f), predominantly (12/13) showed a significant increase in IgG antibody titers against SARS-CoV-2 antigens at follow-up, including spike, RBD and S2 subunit (Supplementary Fig. 7a, b). While 10 infants increased their anti-S2 and anti RB levels by 20 and ~8 fold, respectively, no increase in antibody titers against any of the seasonal alpha- or beta-coronaviruses was observed, despite the presence of detectable IgG levels at baseline (T0); (Supplementary Fig. 7c). Given the young age of this cohort (median 3.2 months), we hypothesized that the detected antibodies against seasonal coronaviruses are likely of maternal origin.

Finally, although anti-IFN autoantibodies have been identified in adults with severe COVID-19³², none of the patients (pCoV, n = 16) or pHCs (n = 6) in this pediatric cohort demonstrated anti-IFN activity in either the acute or convalescent serum samples (Supplementary Fig. 7d).

Lymphocytes from SARS-CoV-2 infected infants express a broader ISG signature compared to those from infected adults

We next compared our infant data with a scRNA-seq PBMC dataset (GSE161918) from a reported cohort of 28 adults with COVID-19 and controls¹⁷. For this comparison, we included data from the adult patients’ first time-point (T0), as well as from 11 healthy matched controls (aHC)¹⁷. After concatenating adult (n = 39) and infant (n = 40) datasets, we applied our pipeline (see methods for details) on doublet-cleaned data, which included pre-processing, batch correction (using Harmony³³), and unsupervised clustering. Based on the criteria we used for the infant cohort (excluding viral loads, which were not available for the adult dataset), we categorized adult patients into three clinical groups, from lower to higher severity (aG1, aG2 and aG3; Fig. 7a). The first round of clustering generated 28 clusters and seven cell-types (Fig. 7b, c) that were then analyzed separately.

a Pediatric and adult cohorts’ information. b UMAP plot representing clusters (n = 28) from 408,281 PBMCs and 79 samples (40 infants and 39 adults) c Cluster annotation. Dot plot represents expression values of selected genes (x-axis) across each cluster (y-axis). Dot size represents the percentage of cells expressing the marker of interest. Color intensity indicates the mean expression within expressing cells. d UMAP plots representing CD14⁺ monocyte subclusters (SCs; n = 4) and their respective annotations. e Bar plot showing the cell frequency of each CD14⁺ monocyte SC across the 79 individuals, as categorized by patient groups: pG1, pG2 and pG3; in purple, pHC in dark green, aG1, aG2 and aG3; in blue and aHC in light green. f As in (d) in CD4⁺ T cell SCs (n = 6). g As in (e) in CD4⁺ T cell SCs. h As in (d) in B cell SCs (n = 7). i As in (e) in B cell SCs. j Boxplot comparing the proportion of ISG^hi CD14⁺ monocyte (SC1), CD4⁺ T cell (SC1) and B cell (SC3) subclusters across the individual (n = 79) in both pediatric (n = 40) and adult (n = 39) cohorts. P values were calculated using a two-sided t-test comparing the mean (G1, G2 or G3) with mean (pHC). ns: non-significant. The upper and lower bounds represent the 75% and 25% percentiles, respectively.

Full size image

CD14⁺ monocytes included four subclusters (SCs), which were categorized into: (i) ISG^low (SC0), (ii) ISG^hi (SC1), (iii) ISG^hi inflam⁺ (SC2) and (iv) CD163⁺ IL1R2⁺ (SC3) (Fig. 7d and Supplementary Fig. 8a, b). Cell composition analysis of CD14⁺ monocytes showed a decrease of SC0 (ISG^lo) and a switch to SC1 (ISG^hi) and/or SC2 (IL1B⁺ ISG^hi) in pG2/pG3 and aG2/G3 groups compared to their respective age-matched controls. Interestingly, we replicated the increase in SC3 (IL1R2⁺ CD163⁺) in patients receiving steroid therapy (Fig. 7e). The analysis of CD16⁺ monocytes showed a switch to an ISG^hi state in both G2/G3 infants and G3 adults (Supplementary Fig. 8c–e).

Subclustering CD4⁺ T cells generated six SCs (Fig. 7f), annotated as: (i) naive (SC0 and SC1; SELL and CCR7), with SC0 being ISG^low and SC1 ISG^hi, (ii) memory (SC2 and SC4, S100A4), (iii) SOX4⁺CD4⁺ T cells (SC3; SOX4) and (iv) regulatory T cells (SC5, FOXP3 and IL2RA). SC4 (CTLs) exhibited the cytotoxic cell-associated chemokine CCL5 (Supplementary Fig. 8f). Cell composition analysis showed: (i) an increase in numbers of naive cells in pediatric healthy controls (pHC) compared to adult healthy controls (aHC), (ii) a reduction of naive cells in pG2/pG3 compared to pHC/pG1, (iii) an increase of memory cells (SC2) in aHC compared to pHC, (iv) increase of CTLs (SC4) in aHC, and (v) an increase of SOX4⁺ CD4⁺ T cells in pHC. Total Tregs (SC5) did not significantly change upon infection in either infants or adults. The ISG^hi CD4⁺ T cell SC (SC1) was present in both infected infants (2/10 pG1, 8/11 pG2, and 5/5 pG3) and adults (6/21 aG3 patients), although their proportion was higher in infants (Fig. 7g).

B cell subclustering generated seven SCs (Fig. 7h), annotated as: (i) naive (SC0, CCR7), (ii) activated (SC1, CD69) (iii) memory (SC2, S100A4) and (iv) transitional (SC4; CD9, MME (CD10)). SC3 encompassed ISG^hi B cells and SC5 upregulated CD83 (Supplementary Fig. 8g). Cell composition analysis showed: (i) increase of naive B cells in pHC compared to aHC, (ii) an increase of memory B cells in aHC compared to pHC, and (iii) a subtle ISG^hi increase in infected adults (aG3) contrasting with a near-complete ISG^hi status in the infected infants with severe disease (pG2/pG3) (Fig. 7i).

CD8⁺ T cell subclustering generated five SCs, classified as: (i) naive (SC0; RPLs), (ii) TEMRA (SC1; GZMB), (iii) Gzk (SC2, GZMK), (iv) MAIT (SC3, KLRB1 and ZBTB16) and (v) proliferative CD8 T cells (SC4, MKI67) (Supplementary Fig. 9a). As expected, cell composition analysis showed increased naive CD8⁺ T cells in pHC and TEMRA/MAIT cells in aHC (Supplementary Fig. 9b). NK cell subclustering generated two SCs (Supplementary Fig. 9c). SC0 upregulated FCGR3A (CD16⁺ NK), and SC1 was XCL1⁺; no disease severity-associated differences were detected (Supplementary Fig. 9d).

Overall, these comparative analyses (Supplementary Fig. 10a–d) showed that CD4⁺ T and B cells from infants with severe disease almost entirely switch to an ISG^hi state—while this is only observed in a minority of T and B cells in adults with severe disease (Fig. 7j).

In this study, we characterized the immune response to SARS-CoV-2 infection during the first weeks of life. By assembling a cohort of young infants with varying disease severity—ranging from subacute (G1) to moderate (G2), and severe (G3) disease—we were able to perform an in-depth analysis of disease severity-associated changes in cell composition and transcriptional profiles at the single cell level, serum cytokine concentrations, and antibody responses. Unlike previous studies¹¹, our cohort focused on very young infants, covering an age range of 3 days to 19 months, and median age of ~2 months. We found disease severity-associated cell composition and transcriptional changes in almost all immune cell types investigated, along with increased concentrations of inflammatory cytokines, IFNγ, and IFN-inducible cytokines. Importantly, these infants mounted consistent antibody responses to SARS-CoV-2 antigens in the absence of anti-IFN autoantibodies.

During the initial stages of the pandemic, it was believed that infants and children were particularly resistant to COVID-19. However, subsequent studies showed that they not only get infected but can also develop moderate to severe disease³⁴. Despite this, these patient populations have been largely understudied, therefore there is limited information on the immune responses of infants who develop severe versus mild disease. A recent study¹¹ analyzed the mucosal and systemic immune response to SARS-CoV-2 infection in infants and young children with mild disease not requiring hospitalization. Children in that cohort showed durable antibody responses, and a robust mucosal immune response characterized by the production of inflammatory cytokines, type I IFN, as well as detection of T helper (Th) 17 and neutrophil-derived proteins (e.g., IL-17, IL-8, and CXCL1). In the systemic compartment, these children had increased chemokine and IFNα concentrations that correlated with viral loads and with marked ISG expression in myeloid cells¹¹. However, no inflammatory cytokines were detected in the blood.

Our study revealed that infants with moderate/severe disease have a distinct immune response, including a prominent ISG transcriptional response across all PBMCs. These transcriptional differences were partially anticipated, as bulk RNA-seq analyses have shown a similar correlation between disease severity and transcriptional IFN responses in infants with other respiratory viral infections^35,36.

The prominent ISG signature in our cohort was expressed in both the myeloid and lymphoid compartments. Among CD14⁺ monocytes, three distinct cell states were identified according to the co-expression of ISGs and pro-inflammatory IL-1B transcripts: (i) ISG^low + IL1-B, which was increased in pHC and G1, (ii) ISG^high + IL1-B, which was increased in G2 patients, and (iii) ISG^high only, which was increased in G3 patients. Several studies suggest a cross-regulation between IL-1B and IFN³⁷, but monocytes co-expressing ISGs and IL-1 were reported in SARS-CoV-2 infected adults³⁸, as well as in pediatric patients with systemic lupus erythematosus (SLE)^39,40. This suggests that SARS-CoV-2 may, directly or indirectly, alter this cross-regulation, a finding that needs further mechanistic exploration.

Lymphocytes from COVID-19 infants showed an ISG signature far more prominent than that observed in adults. This suggests that infants might either produce more type I IFN than adults or respond more robustly to it. Such an enhanced IFN response might represent an innate mechanism of protection for infants against viral infections, particularly while their adaptive immune system is still developing. Although the ISG^hi CD4⁺ T cell subset (CD4_SC1) was mostly contributed by moderate to severe (G2/G3) patients, it was also present at low frequency in every pHC, perhaps indicating a constitutive activation of the IFN pathway in the steady state of young infants (Supplementary Fig. 4a, b).

Circulating CD14⁺ monocytes with macrophage markers (CD163, IL1R2, and FKBP5) were another notable feature observed in G3 infants. IL-1 Receptor antagonist (IL1R2) is a secreted decoy receptor for IL1A, IL1B and IL1RN that modulates IL1 activities⁴¹ and reduces IL-1-mediated inflammation. A similar CD163⁺ IL1R2⁺ population was increased in 7 out of 16 adult patients receiving corticosteroid treatment¹⁷. Indeed, two infants on our cohort who were treated with systemic steroids (dexamethasone) exhibited this CD163⁺ IL1R2⁺ CD14⁺ monocyte subset. It is tempting to speculate that these cells could contribute to a “return to homeostasis”. COVID-19 infected infants also displayed a reduced frequency of pDCs, as found in adult patients. However, unlike adult pDCs, which expressed a signature of apoptosis¹⁷, the pDCs in infants exhibited an ISG signature.

Analysis of the infant’s lymphoid cell populations showed that, as expected given their young age, the majority displayed a naive phenotype. However, a memory compartment was detected within CD4⁺ T cells, including Tregs and CD4⁺ T cells expressing cytotoxicity-associated markers. The distribution of GZM transcripts was particularly intriguing, with broad GZMM expression, co-expression of GZMK and PRF in a small memory subset, and lack of detection of GZMB. Contrary to previous reports⁴², we did not detect an increase in the frequency of CXCR5⁺ CD4⁺ Tfh cells in infants with moderate to severe disease. While this discrepancy might be due to methodological variations, our data suggest an enhanced role of the extrafollicular B cell differentiation pathways, similar to what was reported in adults^28,43.

Regarding CD8⁺ T cells, the most remarkable finding was the expansion of cycling MKI67 CX3CR1⁺ CD8⁺ T cells up-regulating a cytotoxic program including GZMB and PRF1 in infected infants, a cell population that has, until now, only been identified in children with multisystem inflammatory syndrome in children (MIS-C)²⁵. While in children with MIS-C the presence of these cells was reported to correlate with the use of vasopressor drugs²⁵, only one infant in our cohort received such medication, suggesting alternative induction pathways. Further, studies assessing the SARS-CoV-2-specificity of these unique T cells could provide valuable insights into the origin and role of this population.

The B cell compartment was as expected predominantly composed of naive B cells, with limited numbers of memory cells. Yet, the magnitude of the B cell ISG signature was remarkable. ISG^hi subclusters were identified in 65% (17/26) of the infected infants, though also detected at low frequency in every pHC, perhaps indicating a tonic activation of the IFN pathway in the steady state (Supplementary Fig. 4c, d). In addition, a small subset of memory B cells including conventional switched, extrafollicular memory as well as plasma cell markers (B-SC4) was observed in these infants. This might be indicative of an infant-specific memory population that is prone to differentiate into plasma cells. Overall, infants with COVID-19 showed a robust antibody response to the different SARS-CoV-2 antigens, regardless of disease severity. We also detected antibodies against beta coronaviruses, likely of maternal origin, however, contrary to previous hypotheses, these pre-existing antibodies from endemic coronaviruses did not appear to prevent the infection in these young infants. In addition, none of the infants in our cohort displayed anti-interferon antibodies.

In conclusion, our analysis of infants with severe disease revealed a robust ISG signature across most immune cell types, as well as a strong systemic inflammatory response evidenced by increased transcription and high serum concentrations of inflammatory cytokines. Co-expression of ISGs and inflammasome-related genes, especially in monocytes, was especially a feature of infants with severe disease¹¹.

At this point, however, it is yet unclear whether the strong ISG signature observed in young infants has a protective or pathologic role, or possibly both. The presence of a tonic ISG^hi signature in young healthy infants is appealing to hypothesize that it might have a protective role, and that the stronger signature observed in infants with severe disease suggest a stronger effort to control the infection. Earlier studies using bulk transcriptome analyses showed that influenza infection induced a stronger ISG signature than respiratory syncytial virus (RSV) and rhinovirus infections, and that influenza infection was associated with higher fever but generally less severe pulmonary disease⁴⁴. Further longitudinal studies in infants infected with different respiratory viruses that analyze both the systemic and mucosal responses are needed to provide improved mechanistic understanding of the functional significance of the interferon responses in early life immunity.

Our study has the limitations associated with the difficulty of studying infants longitudinally, especially those severely ill, which explains the modest cohort size of the group with severe disease. Other studies have described biphasic early and late immune responses to SARS-CoV-2, pointing out to the importance of timing in immune studies¹⁷. For this reason, we accounted for duration of viral exposure and/or symptoms, which allowed us to classify infants in three groups and investigate the differences in immune responses by multidimensional analyses. Infants in G1 were enrolled later after SARS-CoV-2 exposure compared with infants in G2 and G3, which could partly explain the more prominent ISG response in these two latter groups. Nevertheless, accounting for differences in the duration of infection allowed us to document that ISG responses are influenced by time, and not only by disease severity. Except for follow-up serologic studies, samples for scRNA-seq and serum cytokine analyses were obtained at a single time point during the hospitalization. Finally, we did not measure the impact of SARS-CoV-2 infection at the mucosal level, the primary site of infection.

Despite these limitations, our study detected and underscored distinct immune alterations triggered by SARS-CoV-2 infection in young infants and identified important differences compared with infants with mild disease and with infected adults. The magnitude and significance of the ISG signature across most immune cell populations, coupled with the capacity to mount a significant antibody response to the virus despite a predominantly naive adaptive immune compartment, are among the most remarkable and unexpected findings. Whether these differences are generalizable to other viruses such as influenza or RSV warrants future studies. These insights may have important implications for developing more effective preventive and therapeutic strategies for this vulnerable population.

In summary, our multidimensional approach leveraging small volume blood samples from infants with COVID-19 and variable disease severity provides valuable data to better understand the complex immune alterations that take place during SARS-CoV-2 infection in this age group. Future studies including correlations with mucosal responses are warranted.

Study design

A convenience sample of children <2 years of age hospitalized with COVID-19 were prospectively enrolled at Nationwide Children’s Hospital (NCH) in Columbus, Ohio, USA. Blood (2.5–5 mL) and nasopharyngeal samples were obtained within a median (25–75% interquartile range [IQR]) of 31 [12–57] hours of hospitalization for multiomic analyses and SARS-CoV-2 quantitation by real-time polymerase chain reaction (PCR)⁴⁵. Children with multisystem inflammatory syndrome (MIS-C), or those in whom confirmatory nasopharyngeal SARS-CoV-2 PCR was negative were excluded from the study. At enrollment, we collected demographic and clinical information using a standardized questionnaire designed for the study, and information transferred to a secure database (REDCap). The information collected included duration of illness or time since exposure to a COVID-19 case, and standard parameters of disease severity including oxygen administration, pediatric intensive care unit (PICU) admission, duration of hospitalization, and administration of COVID-19 directed therapies including remdesivir or systemic steroids. During the study period all children hospitalized at NCH, irrespective of the presence of symptoms, underwent nasopharyngeal SARS-CoV-2 testing using a PCR assay per standard of care as described⁴⁵. As a refence for all immune assays, we also included in the study a cohort of age-matched healthy control infants with no respiratory symptoms or treated with antibiotics within two weeks of enrollment. All healthy controls were enrolled pre-pandemic. Healthy controls were typically enrolled in the operating room while undergoing minor scheduled surgical procedures not involving the respiratory tract, or at the primary care offices during well-child visits. The study was approved by the Nationwide Children’s Hospital (NCH) IRB (18-00591). Written informed consent was obtained from all children’s guardians before study participation. Further clinical details of all participants are summarized in Supplementary Data 1.

Patient categorization

Infants in G1 were mostly hospitalized with fever and/or respiratory symptoms, all had low viral loads at admission, and a median of 21 days since exposure to an index SARS-CoV-2 positive case. Infants in G2 were also hospitalized with fever and respiratory symptoms but had high viral loads, ~2 days of symptoms at the time of hospitalization, and a median of 7 days since exposure to an index SARS-CoV-2 positive case. None of these children required supplemental oxygen and one received intravenous fluid. Infants in G3, were similar to G2 in terms of duration of illness at admission and exposure to an index case, they also had high viral loads but demonstrated a more severe clinical phenotype requiring PICU care and/or respiratory support. Two of these children also received steroids and remdesivir for the management of severe COVID-19.

PCR assays for SARS-CoV-2 viral loads

Nasopharyngeal (NP) swabs were collected at enrollment, placed in viral transport media, transported immediately to the laboratory, aliquoted and stored at −80 °C. Similarly, blood samples were collected in EDTA tubes (BD Vacutainer, Franklin Lakes, NJ, USA), centrifuged, and aliquots of plasma were stored at −80 °C until processed in batches. Viral RNA was extracted from 200 microliters of NP or plasma samples using the QIAcube HT instrument (Qiagen Inc, Germantown, MD, USA) and eluted into 100 microliter volume. In brief, SARS-CoV-2 viral loads were measured using a two-step reverse-transcription (RT) quantitative PCR assay targeting the conserved region of the N1 gene using described primers and probe⁴⁶. A standard curve was generated from a dilution series of a known concentration(10⁷ copies/mL) of the N1 gene of SARS-CoV-2. Standards and negative controls were included and tested with each PCR assay. The lower limit of detection of the assay was 1–2 copies/reaction (~200 copies/mL).

Cytokine assays

Samples were analyzed at Nationwide Children’s Hospital, following the instructions of the manufacturer as described⁴⁷. We applied the Olink (Boston, MA) inflammation panel that includes 92 proteins (Supplementary Data 7). Samples were analyzed in the Olink Boston Laboratory.

Enzyme-linked immunosorbent assay (ELISA)

Ninety-six-well microtiter plates (Thermo Fisher) were coated at a concentration of 2 μg/mL of recombinant protein overnight (4 °C): SARS-CoV-2 full-length spike, RBD and S2 subunit. For seasonal coronaviruses, full-length and S1 subunit was used for beta- OC43 and HKU1; and full-length spike for alpha- 229E. The next day, plates were washed three times with PBS (phosphate-buffered saline; Gibco) containing 0.1% Tween-20 (T-PBS, Fisher Scientific) using an automatic plate washer (BioTek). After washing, the plates were blocked for 1 h at room temperature with 200 ml blocking solution (PBS-T with 3% (w/v) milk powder (American Bio)) per well. The blocking solution was removed, and serum samples diluted to a starting concentration of 1:80, serially diluted 1:3 in PBS-T supplemented with 1% (w/v) milk powder and incubated at room temperature for 2 h. The plates were washed three times with PBS-T and 50 ml anti-human IgG (Fc-specific) horseradish peroxidase antibody (HRP, Sigma, #A0170) diluted 1:20,000 in PBS-T containing 1% milk powder was added to all wells and incubated for 1 h at room temperature. After washing, 100 μL of peroxidase substrate (3,3′,5,5′-Tetramethylbenzidine, TMB, Rockland) was added and incubated at room temperature for 30 min. The reaction was stopped with 1 N sulfuric acid solution (Fisher Science) and the plates read at a wavelength of 450 nm with a plate reader (BioTek). The results were recorded in Microsoft Excel and AUC values were computed by plotting normalized optical density (OD) values against the reciprocal serum sample dilutions for ELISAs in GraphPad Prism.

Functional evaluation of anti-IFN-I auto-antibodies

A549 cells (ATCC; CRM-CCL-185) were cultured in Dulbecco’s modified Eagle medium (DMEM; Gibco) supplemented with 10% fetal bovine serum (FBS; PEAK) and penicillin-streptomycin (Gibco), at 37 °C, under an atmosphere containing 5% CO₂. Cells were checked periodically and tested negative for mycoplasma contamination. A549 cells were seeded in 96-well plates at a density of 3 × 10³ cells/well. The next day, plasma samples or a commercial anti-IFNα2 antibody (catalog number 21100-1; R&D systems) were serially diluted (10-fold) and incubated with 50 pg/ml recombinant IFN-α2 (catalog number 11101-2; R&D systems) for 1 h at 37 °C (starting concentration: plasma samples = 1/100 and anti-IFN-α2 antibody = 1/1000). Following this incubation period, the cell culture medium in the 96-well plates was removed by aspiration and replaced with the plasma/antibody-IFN-α2 mixture. Each sample was tested once, in triplicate. The plates were incubated overnight and the plasma/antibody-IFN-α2 mixture was removed by aspiration. The cells were then washed three times with PBS to remove potential anti-influenza neutralizing antibodies and infected with a recombinant Cal/09 virus expressing NS1-mCherry (CalNSmCherry) at MOI = 0.5. One day after infection, cells were fixed with 4% formaldehyde, washed twice with PBS, and stained with DAPI. The percentage of infected cells was quantified using the Celigo (Nexcelcom) imaging cytometer. All antibody assays were performed at Mt. Sinai School of Medicine.

Blood preparation for single cell RNA sequencing (scRNA-seq)

PBMCs were thawed quickly at 37 °C and into DMEM supplemented with 10% FBS. Cells were quickly spun down at 400 × g, for 10 min. Cells were washed once with 1 x PBS supplemented with 0.04% BSA and finally re-suspended in 1 x PBS with 0.04% BSA. Viability was determined using trypan blue staining and measured on a Countess FLII. Briefly, 12,000 cells were loaded for capture onto the Chromium System using the v2 single cell reagent kit (10X Genomics). Following capture and lysis, cDNA was synthesized and amplified (12 cycles) as per manufacturer’s protocol (10X Genomics). The amplified cDNA was used to construct an Illumina sequencing library and sequenced on a single lane of a HiSeq 4000.

Single-cell Raw data processing and data combining

Illumina basecall files (*.bcl) were converted to fastqs using cellranger v3.0.2, which uses bcl2fastq v2.17.1.14. FASTQ files were then aligned to hg19 genome and transcriptome using the cellranger v3.0.2 pipeline, which generates a gene—cell expression matrix. The samples were combined using cellranger aggr from cellranger, which aggregates outputs from multiple runs, normalizing them to the same sequencing depth (normalize=mapped) and then re-computing the gene-barcode matrices and analysis on the combined data. The sequencing information for individuals included in the study (e.g., number of reads per cell or UMI counts per cell) are shown in Supplementary Data 2.

Scrublet for multiplet prediction and removal

Generally, we expected about 2 to 8% of the cells to be hybrid transcriptomes or multiplets, occurring when two or more cells are captured within the same microfluidic droplet and are tagged with the same barcode. Such artifactual multiplets can confound downstream analyses. We applied Scrublet¹² python package to remove the putative multiplets. Scrublet assigns each measured transcriptome a ‘multiplet score’, which indicates the probability of being a hybrid transcriptome. Multiplet scores were determined for each individual (using the raw data). Before/after multiplet removal, the number of cells was 313,076/300,708.

Single-cell pre-processing, dimension reduction, graph-based clustering, and cluster annotation

The cleaned (after multiplet removal using scrublet¹²) aggregated matrices were fed into the Python-based Scanpy⁴⁸ workflow (https://scanpy.readthedocs.io/en/stable/), which includes pre-processing, visualization, clustering and differential expression testing. The pipeline we applied was inspired by the Seurat R package workflow and can be found here: https://github.com/dnehar/Infants_Cov19.

Quality control and cell filtering

We applied the following filtering parameters: (i) all genes that were not detected in ≥3 cells were discarded, using pp.filter_genes function, (ii) cells with less than 400 total unique transcripts were removed prior to downstream analysis using pp.filter_cells function, (iii) cells in which >25% of the transcripts mapped to the mitochondrial genes were filtered out, as this can be a marker of poor-quality cells and (iv) cells displaying a unique gene counts >4500 genes were considered outliers and discarded.

Data normalization, log transformation and scaling

We normalized the data using the pp.normalize_per_cell function. Thus, library-size normalization was performed based on gene expression for each barcode by scaling the total number of reads per cell to 10,000. We log-transformed the data using the pp.log1p function and scaled to unit variance using pp.scale function (with the following parameters: max_value=10). The 1362 highly variable genes (HVG) were identified using filter_genes_dispersion function (with the following parameters: min_mean=0.0125, max_mean=3, min_disp=0.5).

Linear dimensional reduction using PCA

To reduce the dimensionality of the data, we ran principal component analysis (PCA) using tl.pca function, which reveals the main axes of variation and denoises the data. The contribution of each PCs to the total variance was assessed using pl.pca_variance_ratio function.

Neighborhood graph computing, embedding and clustering

The neighborhood graph of cells was computed based on the PCA representation of the data matrix, using pp.neighbors function (with the following parameters: n_neighbors = 10, n_pcs = 40). The neighborhood graph was then embedded using UMAP⁴⁹ (tl.umap function) and visualized using pl.umap function. We finally used the Leiden graph-based clustering using tl.leiden function with resolution = 1.2 (stored in ‘Res1_2_BC’ column).

Batch effect correction

To account for technical source of variation, such as ‘10X runs’ or ‘Years’, we applied a batch effect correction using BBKNN.

bbknn.bbknn function was ran using the following parameters: metric=“angular”, approx=True, neighbors_within_batch=5, n_pcs=20, trim=50 and copy=True.

All the analyses presented in this study were done on the BBKNN corrected data. UMAP visualization and Leiden clustering were then computed on the BBKNN corrected using the following parameters:

• sc.tl.umap(BBKNN_corrected_object, min_dist=0.3, n_components=3)

• sc.tl.leiden(BBKNN_corrected_object, Resolution=1.2, key_added=‘Res1_2_AC’)

To account for batch effects in the pediatric-adult (paCov) dataset, we used harmony³³.

IFN score calculation

ISGs (from previously described interferon-related modules⁵⁰, Supplementary Data 3) were used to score IFN expression in each cluster/subcluster generated in the pCov40 dataset. To do so, we calculated the mean expression for each cell, within each cluster/subcluster using the h5ad object (adata), as follow: adata.obs[‘ISG_score’] = adata.X[:,IFN_markers].mean(1).

Finding marker genes/evaluation of cluster identity

To annotate the clusters generated from the BBKNN corrected object, we used both differential expression analysis between clusters and classification based on putative marker gene expression. We applied tl.rank_genes_groups function to compute a ranking for the differential genes in each cluster/Subcluster, comparing each cluster to the rest of the cell using Wilcoxon test. We only considerate clusters/Subclusters that showed distinct transcriptomic programs. The top 100 marker genes for cluster and subclusters, in pCoV are included in Supplementary Data 4 and 6, respectively. The top 100 marker genes defining infected infants compared with healthy across PBMCs and immune subsets are included in Supplementary Data 5. The top 10 marker genes were visualized using the sc.pl.rank_genes_groups_matrixplot() function.

Subclustering parameters and data cleaning

A script showing the subclustering process can be found here: https://github.com/dnehar/Infants_Cov19. We only considered SC defined by distinct gene sets, by merging similar ones. Based on the number of cells and to avoid over-clustering, we used different clustering resolutions. The following parameters of resolution were used: 0.4 (CD14⁺ and CD16⁺ monocytes, NK cells, CD4⁺ T cells, MGK), 0.6 (B cells, CD8⁺ T cells and cDCs), 0.8 allowed the separation of PC and pDC (C24_PCs_pDC in Fig. 1d). We discarded two small SC of CD14⁺ monocytes (one PF4⁺ and the other was IGKC⁺). cDC (n = 504 cells) which initially mingled with CD14⁺ monocytes were separated after the subclustering step. Suspected doublets in NK cells (PF4⁺ SC), B cells (NKG7⁺, GNLY⁺ SC), MGK (T and B markers), Eryth (T, B, NK markers) were discarded from the analysis. Overall, the number of cells after multiplet removal (n = 300,708), the filtration (n = 256,766), and visual inspection steps was 203,402 cells. A similar strategy was applied to the pediatric and adult combined (paCoV) cohort. Number of cells before filtration was 499,996, after filtration 440,494 cells and after visual inspection steps was 408,281 cells.

Subclustering of memory CD4⁺ T cells and proliferative CD8⁺ T cells

Memory CD4⁺ T cells (CD4-SC3) and proliferative CD8⁺ T cells (CD8-SC4) were extracted and underwent a clustering pipeline as described previously, including the redefinition of variable genes, PCA analysis, Neighborhood graph computing, embedding and clustering and Batch effect correction.

Adult COVID-19 data Integration (paCoV dataset)

We first downloaded scRNA-seq raw counts (GSE161918¹⁷) from an adult cohort, which included 28 patients (aCov; we only selected the time 0), as well as 11 healthy controls (aHC). After converting the Seurat object into an anndata object (.h5ad), we concatenated adult (28 aCov collected at time 0 and 11 aHC) and pediatric (14 pHC and 26 pCov) raw counts, resulting in an object including 79 samples. We then cleaned the data from doublets using Scrublet¹² and applied a computational pipeline, including pre-processing, batch correction (using Harmony³³) and non-supervised clustering.

Comparison with Wimmers et al. cohort¹¹

We downloaded the single cell Multiome dataset (cell ranger outputs) from GEO (GSE239799¹¹). We then concatenated the raw counts (n = 43) with our pediatric cohort (n = 40), resulting in anndata object (.h5ad) including 83 samples. We then cleaned the data from doublets using Scrublet¹² and applied a computational pipeline, including pre-processing, batch correction (using Harmony³³) and non-supervised clustering. The Wimmers et al. cohort included 9 longitudinal samples (pre-infection or “pre” – acute non-Omicron or “acute” and convalescent or “conv”), 1 acute non-omicron, 8 acute Omicron and 7 healthy controls¹¹. Pseudo-bulk analysis was performed after aggregating raw counts per donor across different immune subsets (Supplementary Data 5). ISG counts were then plotted in each donor, as categorized by clinical groups.

Statistical analysis

Statistical analysis was performed using R/4.0.2. Tests were used to determine data distribution and depending on the normality of the data, comparisons were performed using two-sided Student t tests (for two groups, parametric) or the non-parametric the Wilcoxon signed rank test (for two groups, paired) with two-tailed P values unless otherwise stated. Differences were considered significant when P < 0.05 (*), P < 0.01 (**), P < 0.001 (***) and P < 0.0001 (****).

Python module versions

scanpy==1.7.1 anndata==0.7.5 umap==0.4.6 numpy==1.19.2 scipy==1.5.2 pandas==1.1.3 scikit-learn==0.24.1 statsmodels==0.12.2 python-igraph==0.8.3 leidenalg==0.8.3.

Reporting summary

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