Systems biological assessment of the temporal dynamics of immunity to a viral infection in the first weeks and months of life

Cohort and study design

Specimens for this analysis were collected as part of the ongoing IMPRINT Influenza Cohort, an NIH-funded longitudinal birth cohort of healthy mothers and children in Cincinnati, Ohio. Pregnant women were screened for enrollment at three local delivery hospitals, and mothers were enrolled in the third trimester of pregnancy. Starting at two weeks of life, research coordinators trained mothers to complete twice weekly SMS surveys for their child to report on the presence or absence of symptoms common to respiratory infections. Mothers were also trained to collect midturbinate nasal swabs from themselves and from their children, which were analyzed at an IMPRINT laboratory. Study participants were asked to collect at least one nasal swab each week. An additional swab was submitted if either the mother or the child was symptomatic of respiratory illness. Among children enrolled in this cohort, we included 19 pre-omicron SARS-CoV2 cases and 22 omicron SARS-CoV2 cases. Cases were selected based on the timing of pre-infection, acute and convalescent blood collections, and availability of PBMCs at these time points. A set of 30 children, age- and sex-matched to pre-Omicron child cases, were selected as controls; age-matching was based on the age of the case at the time of convalescence. COVID-19 infections were considered symptomatic if the child experienced any symptom at least once, within −7/21+ days of their first COVID-19 positive swab, as reported on weekly surveys. Illness severity was evaluated by the presence of cough, fever or difficulty breathing, duration of cough and/or fever, and receipt of medical attention. The presence or absence of symptoms on the day of blood collection was determined by responses on the weekly survey completed closest to the blood collection date (−3/+3 days of blood collection). If a weekly survey response was unavailable, responses from the visit questionnaire were used (n=1). Demographic data is shown in DataS1.

SARS-CoV-2 RT-PCR

Nasal swab samples underwent nucleic acid extraction using a custom protocol and the QIAamp 96 Virus Kit on the QIAGEN QIAcube HT instrument. Extracted RNA was run on Applied Biosystems 7500 or 7500 Fast instruments using CDC-developed real-time RT-PCR assays for Influenza ^{48, 49}, SARS-CoV-2 ^{50, 51}, or Influenza and SARS CoV-2 in the Flu SC2 Multiplex assay ^{52, 53}. Nasal swabs were run on the influenza assay from November 2019-March 17, 2020. SARS-CoV-2 PCR was run for nasal swabs collected starting February 1, 2020, and was used in the lab through December 31, 2021. Starting on January 3, 2022, the Flu SC2 assay has been used for SARS-CoV-2 testing. In addition, all infant nasal swabs underwent multiplex detection of 21 respiratory pathogens using the Luminex NxTAG-RPP assay ⁵⁴, including Influenza A and B. All samples positive for Influenza A and/or B underwent additional RT-PCR detection for Influenza A subtype and Influenza B lineage ^{48, 49}.

SARS-CoV-2 sequencing

Nasal swab samples underwent nucleic acid extraction using the QIAGEN Viral RNA Mini Kit (Qiagen, Inc) on the QIAGEN QIAcube Connect instrument following the manufacturer’s recommendations. RNA samples were subjected to direct sequencing using the modified ARTIC3 protocol ((https://artic.network/ncov-2019) with the addition of primer booster sets as implemented in the Qiagen QiaseqDirect protocol. RNA was subjected to random primed cDNA synthesis followed by amplification in two pools of multiplexed primer sets resulting in overlapping amplicons spanning the entire genome. Subsequently, 24 cycles of polymerase chain reaction were utilized to add dual index primers and amplify SARS-Cov-2 amplicons. DNA concentrations were normalized, samples were pooled and then subjected to sequencing to a depth of at least 100,000 reads per sample utilizing paired 150 nucleotide reads on an Illumina NextSeq 500 sequencing machine (Illumina, Inc).

Raw sequence data were demultiplexed and then aligned against the ancestral Wuhan-1 genome (Accession MN908947) ⁵⁵ using bwa-mem ⁵⁶. Samtools commands “sort”, “index”, “view”, and “mpileup” ⁵⁷ were applied sequentially, and the ivar “consensus” command ⁵⁸ was used to output a consensus sequence. Variant identification and lineage calling were performed with the software Phylogenetic Assignment of Named Global Outbreak Lineages (Pangolin) version 4.0.4 ⁵⁹.

Sample processing

Peripheral blood was collected from study children at two and six weeks of life, every summer and acutely following “events.” An event included receipt of influenza or COVID-19 vaccines or having a nasal swab test positive for influenza or SARS-CoV2; participants completed an additional symptom survey at event visits. Up to 16mL of collected blood was deposited into sodium citrate Mononuclear Cell Preparation tubes (CPT) and promptly delivered to the laboratory for processing. After the initial spin of the CPT tubes, plasma was collected, and up to ten aliquots of plasma at 0.5mL were stored at −65°C or colder. Cells were collected, washed, and counted using the automated Vi-CELL XR Viability Analyzer. Aliquots were made at concentrations of either 2.0×10⁶ or 5.0×10⁶ cells in 1mL of Fetal Bovine Serum (FBS) +10% Dimethyl Sulfoxide (DMSO), depending on the number of cells obtained. Once aliquoted, the cryovials were placed into a Mr. Frosty cooler filled with isopropyl alcohol and placed into a −65°C or colder freezer for 24-72 hours. Samples were then placed into liquid nitrogen storage units for long-term storage.

Adult Cohort

The samples of the adult participants included in the study were collected in 2020 at the start of the pandemic, prior to the advent of vaccines. All samples were collected at the Hope Clinic of the Emory Vaccine Center or at Stanford University. Healthy controls were asymptomatic adults from whom samples were collected prior to the widespread circulation of SARS-CoV-2 in the community. Patients with COVID-19 were defined according to the original WHO guidance and positive SARS-CoV-2 RT-PCR testing by nasopharyngeal swabs as described in our original study ²¹. Patients with COVID-19 were classified as acute (less than 4 weeks from symptom onset or symptomatic at the time of sample collection) or convalescent (more than 4 weeks since symptom onset and asymptomatic or negative SARS-CoV-2 RT-PCR testing and asymptomatic). The severity of inpatient COVID-19 cases was classified based on the adaptation of the Sixth Revised Trial Version of the Novel Coronavirus Pneumonia Diagnosis and Treatment Guidance. Mild/moderate cases were defined as respiratory symptoms with radiological findings of pneumonia. Severe cases were defined as requiring supplemental oxygen. Critical cases were organ failure necessitating intensive care unit (ICU) care. The study received approval from the appropriate Institutional Review Board at Emory (#00022371) and Stanford (#55689) University. All samples were collected after informed consent. Demographic data is shown in DataS1.

Anti-Spike electrochemiluminescence (ECL) binding ELISA

Anti-Spike IgG titers were measured using V-plex COVID-19 panel 23 from Mesoscale Discovery (Cat #K15567U). The assay was performed as per the manufacturer’s instructions. Briefly, the multi-spot 96 well plates were blocked in 0.15 ml of blocking solution with shaking at 700 rpm at room temperature. After 30 min of incubation, 50 μl of plasma samples serially diluted in antibody diluent solution and serially diluted calibrator solution was added to each plate in the designated wells and incubated at room temperature for 2 h with shaking. Plasma samples were assayed at a 1:100 starting dilution and 6 additional 5-fold serial dilutions. After 2 h of incubation, the plates were washed, 50 μl of Sulfo-tag conjugated anti-IgG was added, and the plates were incubated at room temperature for 1 h. After incubation, the plates were washed, and 0.15 ml of MSD-Gold read buffer was added. The plates were immediately read using the MSD instrument. The Meso scale arbitrary light unit signal was used for calculating the area under curve (AUC) in Prism v.9.4.1.

Pseudovirus production and neutralization assay

VSV-based GFP/nanoluciferase-expressing SARS-CoV-2 pseudoviruses were produced as described previously ⁶⁰. VSV-ΔG-GFP/nanoluciferase and plasmids encoding spike genes of SARS-CoV-2 Wuhan (SΔ19), Delta (B.1.617), and Omicron (B.1.529) were provided by Dr. Gene S. Tan (J. Craig Venter Institute, La Jolla, CA). To perform the neutralization assay, Vero E6-TMPRSS2-T2A-ACE2 cells (BEI Resources, NIAID; NR-54970) were seeded at a density of 1×10⁴ per well in half area 96-well black opaque plates (Greiner Bio-One) and were grown overnight at 37°C in a 5% CO2 atmosphere. Serum samples were 5-fold serially diluted using the infection medium (DMEM supplemented with 2% FBS and 100 U/mL Penicillin-Streptomycin) in duplicates. Diluted serum samples were then mixed with an equal volume of Wuhan, Delta, or Omicron pseudoviruses, diluted in infection medium at an amount of 200-400 focus-forming units/mL per well, followed by incubation at 37°C for 1 hour. Subsequently, immune complexes were added onto the monolayers of PBS-washed Vero E6-TMPRSS2-T2A-ACE2 cells and incubated at 37°C. At 18 hours post-incubation, supernatants were removed, cells were washed once with PBS, and nanoluciferase enzymatic activities were measured using the Nano-Glo Luciferase Assay System (Promega) and a SpectraMax iD3 multi-mode microplate reader. Percent inhibition values were calculated by subtracting the percent infection from 100. Non-linear curves and IC₅₀ values were determined using GraphPad Prism.

Neutralization assay - Omicron subtyping

Human samples were evaluated in a qualified SARS-CoV-2 pseudovirion neutralization assay (PsVNA) using SARS-CoV-2 WA1/2020 strain and variants. SARS-CoV-2 neutralizing activity measured by PsVNA correlates with PRNT (plaque reduction neutralization test with authentic SARS-CoV-2 virus) in previous studies ^61–63. Pseudovirions were produced as previously described ⁶¹. Briefly, human codon-optimized cDNA encoding SARS-CoV-2 spike glycoprotein of the WA1/2020 and variants was synthesized by GenScript and cloned into eukaryotic cell expression vector pcDNA 3.1 between the BamHI and XhoI sites. Pseudovirions were produced by co-transfection Lenti-X 293T cells with psPAX2(gag/pol), pTrip-luc lentiviral vector, and pcDNA 3.1 SARS-CoV-2-spike-deltaC19, using Lipofectamine 3000. The supernatants were harvested at 48h post-transfection and filtered through 0.45µm membranes, and titrated using 293T-ACE2-TMPRSS2 cells (HEK 293T cells that express ACE2 and TMPRSS2 proteins). Neutralization assays were performed as previously described ^{36, 62–65}. For the neutralization assay, 50 µL of SARS-CoV-2 S pseudovirions (counting ∼200,000 relative light units) were pre-incubated with an equal volume of medium containing serial dilutions of samples at room temperature for 1h. Then 50 µL of virus-antibody mixtures were added to 293T-ACE2-TMPRSS2 cells (10⁴ cells/50 μL) in a 96-well plate. The input virus with all SARS-CoV-2 strains used in the current study was the same (2 x 10⁵ relative light units/50 µL/well). After a 3 h incubation, fresh medium was added to the wells. Cells were lysed 24 h later, and luciferase activity was measured using One-Glo luciferase assay system (Promega, Cat# E6130). The assay of each sample was performed in duplicate, and the 50% neutralization titer was calculated using Prism 9 (GraphPad Software). Controls included cells only, virus without any antibody and positive sera.

Auto-antibody analysis with Bead-based array

We created two different custom, bead-based antigen arrays (DataS2) modeled on similar arrays that we previously used to study autoimmune diseases, immunodeficiency disorders, and COVID-19 disease in adults ^{11, 66–73}. Antigens were selected based on our published datasets, literature searches that have implicated specific antigens in COVID-19, potential for mechanistic contribution to COVID-19 pathogenesis, and compatibility with bead-based platforms. The autoantigen array contained 60 commercial protein antigens associated with connective tissue diseases and tissues associated with tissue inflammation. The cytokine array comprised 61 proteins including interferons, interleukins, and other cytokines. The arrays were constructed as previously described ⁷³. Antigens were coupled to carboxylated magnetic beads (MagPlex-C, Luminex Corp.) such that each antigen was linked to beads with unique barcodes ^{68, 70}. Briefly, 8 μg of each antigen or control antibody was diluted in phosphate-buffered saline (PBS) and transferred to 96-well plates. Diluted antigens and control antibodies were conjugated to carboxylated magnetic beads per ID. Beads were distributed into 96-well plates (Greiner BioOne), washed, and re-suspended in phosphate buffer (0.1 M NaH2PO4, pH 6.2) using a 96-well plate washer (Biotek). The bead surface was activated by adding 100 μl of phosphate buffer containing 0.5 mg 1-ethyl-3(3-dimethylamino-propyl)carbodiimide (Pierce) and 0.5 mg N-hydroxysuccinimide (Pierce). After 20 min incubation on a shaker, beads were washed and resuspended in activation buffer (0.05 M 2-N-Morpholino EthaneSulfonic acid, MES, pH 5.0).

Diluted antigens and control antibodies were incubated with beads for 2 h at room temperature. Beads were washed three times in 100 μl PBS-Tween, re-suspended in 60 μl storage buffer (Blocking reagent for Enzyme-Linked Immunosorbent Assay, ELISA, Roche) and stored in plates at 4 °C. Immobilization of some antigens and control antibodies on the correct bead IDs was confirmed by analysis using commercially available monoclonal antibodies or antibodies specific for epitope tags. In addition, prototype human plasma samples derived from participants with autoimmune diseases with known reactivity patterns (from ImmunoVision, Stanford Autoimmune Diseases Biobank, and OMRF) were used for validation. Plasma samples were tested at 1:100 dilution in 0.05% PBS-Tween supplemented with 1% (w/v) bovine serum albumin (BSA) and transferred into 96-well plates in a randomized layout. The bead array was distributed into a 384-well plate (Greiner BioOne) by transfer of 5 μl bead array per well. 45 μl of the 1:100 diluted sera were transferred into the 384-well plate containing the bead array.

Samples were incubated for 60 min on a shaker at room temperature and then stored at 4 °C overnight. Beads were washed with 3 × 60 μl PBS-Tween on a plate washer (EL406, Biotek) and 50 μl of 1:500 diluted R-phycoerythrin (R-PE) conjugated Fcγ-specific goat anti-human IgG F(ab’)2 fragment (Jackson ImmunoResearch, Cat # 109-116-098) was added to the 384-well plate for detection of bound human IgG. After incubation with the secondary antibody for 30 min, the plate was washed with 3 × 60 μl PBS-Tween and re-suspended in 60 μl PBS-Tween prior to analysis using a FlexMap3DTM instrument (Luminex Corp.) and Luminex xPONENT® version 4.2 software. A minimum of 50 events per bead ID were counted. Binding events were displayed as Mean Fluorescence Intensity (MFI). To ensure reproducibility and rigor, all samples were run in duplicate in each experiment. Prototype autoimmune sera as described above, were used as positive controls. All data analyses and figure generation were performed using Python, R and R studio. For normalization, average MFI values for “bare bead” IDs were subtracted from average MFI values for antigen-conjugated bead IDs. The average MFI for each antigen was calculated using samples from healthy infants who have not had COVID and have not tested positive for another upper respiratory infection in the last 7 days. Antibodies were considered “positive” if MFI was >5 SD above the average MFI and >3000 units, which is a threshold that we have used previously for publication and is more stringent than commonly published in related literature ^{11, 73}.

Spike-specific memory B cell staining

Cryopreserved PBMCs were thawed and washed twice with 10 mL of FACS buffer (1 x PBS containing 2% FBS and 1 mM EDTA) and resuspended in 100 uL of 1x PBS containing Zombie UV live/dead dye at 1:200 dilution (BioLegend, 423108) and incubate at room temperature for 15 minutes. Following washing, cells were incubated with an antibody cocktail for 1 hour protected from light on ice. The following antibodies were used: IgD PE (BD Biosciences, 555779), IgM PerCP-Cy5.5 (BioLegend, 314512), CD20 APC-H7 (BD Biosciences, 560734), CD27 PE-Cy7 (BioLegend, 302838), CD14 PE/Dazzle™ 594 (BioLegend, 301852), CD16 BV605 (BioLegend, 302040), IgG BV650 (BD Biosciences, 740596), CD3 BUV737 (BD Biosciences, 612750) and Alexa Fluor 488-labeled Wuhan spike (SinoBiological, 40589-V27B-B), and BV421 labeled Omicron Spike (SinoBiological™, 40589-V49H3-B). All antibodies were used as the manufacturer’s instruction and the final concentration of each probe was 0.1 ug/ml. Cells were washed twice in FACS buffer and immediately acquired on a BD FACS Aria III for acquisition and FlowJo for analysis.

Single cell BCR-seq

SARS-CoV-2 spike specific memory B cells gated on singlet CD3-CD14-CD16-CD20+ IgM-IgD-CD27^low/high IgG+ spikes+ were single-cell sorted into individual wells of 96-well plates containing 16 μl of lysis buffer per well using a FACS Aria III. The lysis buffer was composed of 20 U RNAse inhibitor (Invitrogen), 5 mM DTT (Invitrogen), 4 ul 5x RT buffer (Invitrogen), 0.0625 ul Igepal (Sigma), 10 ug/ml Carrier RNA (Applied Biosystems). The 96-well plates went through a quick freeze-thaw cycle, and 0.5 ug Oligo(dT)18 (Thermo Scientific), 0.5 mM dNTP mix (Invitrogen), and 200 U Superscript IV (Invitrogen) was added in a total volume of 4 ul followed by thorough mixing and spinning. The reverse transcription was performed as follows: 10 min at 42 °C, 10 min at 23 °C, 20 min at 50 °C, 5 min at 55 °C, 10 min at 80 °C and finally cooling to 4°C. Ig heavy chain and light chain (kappa/lambda) variable gene fragments were amplified by nested PCR (HotStarTaq DNA Polymerase, QIAGEN) using primer cocktails ^{74, 75} at a concentration of 250 nM per primer. The PCR mix consisted of 2.5 ul 10x PCR buffer, 0.5 ul 10 mM dNTP mix (Invitrogen), 0.5 ul 25 mM MgCl2, 5 ul Q-solution, 1 U HotStarTaq, 0.5 ul 5’ and 3’ primers and 2.5 ul of cDNA. Water was added up to a total volume of 25 ul. The 2nd round PCR products was evaluated on 2% agarose gels, purified using QIAquick spin columns (Qiagen) and sequenced using 2nd round PCR reverse primers. The sequences were analyzed using the online IMGT/HighV-QUEST tool. IGHV and IGLV nucleotide sequences were aligned against their closest germlines and the somatic hypermutation rate was calculated based on the IMGT v-identity output. The average mutation rate was calculated by dividing the sum of all somatic hypermutation rates by the number of sequences used for the analysis in each individual. The Change-O toolkit v.1.0.0. and SHazaM R package were used for B cell clonality analysis ⁷⁶.

T cell stimulation and intracellular cytokine staining assay

Antigen-specific T cell responses were measured using the intracellular cytokine staining assay as previously described ⁷⁷. Live frozen PBMCs were revived, counted, and resuspended at a density of 2 million live cells per ml in complete RPMI (RPMI supplemented with 10% FBS and antibiotics). The cells were rested for 6 h at 37 °C in a CO₂ incubator. At the end of 6 h, the cells were centrifuged, resuspended at a density of 10 million cells per ml in complete RPMI, and 100 μl of cell suspension containing 1 million cells was added to each well of a 96-well round-bottomed tissue culture plate. Each sample was treated with two or three conditions depending on cell numbers: no stimulation or a peptide pool spanning the Spike protein of the ancestral Wu strain or Omicron BA.1 variant (where cell numbers permitted) in the presence of anti-CD28 (1 μg ml⁻¹; clone CD28.2, BD Biosciences) and anti-CD49d (clone 9F10, BD Biosciences), as well as anti-CXCR3 (DataS3). The peptide pools were 15-mer peptides with 10-mer overlaps spanning the entire Spike protein sequence of each variant ⁷⁸. Each peptide pool contained 253 peptides and was resuspended in DMSO at a concentration of 1 mg/ml.

PBMCs were stimulated at a final concentration of 1 μg/ml of each peptide in the final reaction with an equimolar amount of DMSO [0.5% (v/v) in 0.2-ml total reaction volume] as a negative control. The samples were incubated at 37°C in CO2 incubators for 2 hours before the addition of brefeldin A (10 μg ml⁻¹). The cells were incubated for an additional 4 hours. The cells were washed with PBS and stained with Zombie ultraviolet (UV) fixable viability dye (BioLegend). The cells were washed with PBS containing 5% FBS before adding a surface antibody cocktail (Supp Table 5). The cells were stained for 20 min at 4 °C in 100-μl volume. Subsequently, the cells were washed, fixed, and permeabilized with cytofix/cytoperm buffer (BD Biosciences) for 20 min. The permeabilized cells were stained with intracellular cytokine staining antibodies (Supp Table 5) for 20 min at room temperature in 1× perm/wash buffer (BD Biosciences). Cells were then washed twice with perm/wash buffer and once with staining buffer before acquisition using the BD Symphony Flow Cytometer and the associated BD FACS Diva software. All flow cytometry data were analyzed using Flowjo software v.10 (BD Bioscience). DMSO background was subtracted from all samples and the positivity threshold was defined as 3x the median of peptide stimulated samples from healthy control infants.

Quantitation of human IFN-α2a in plasma

Human IFN-α2a was measured using an S-PLEX Human IFN-α2a kit from Mesoscale Discovery (Cat # K151P3S). The assay was performed as per the manufacturer’s instructions. Briefly, the uncoated 96 well plates were washed 3 times with 150 μl per well of 1xMSD wash buffer and coated with 50 μl of coating solution per well with shaking at 700 rpm at room temperature. After 1 hour of coating, the plates were washed, and 25 μl of blocking solution was added to each well, followed by adding 25 μl of neat plasma samples or serially diluted calibrator solution to each plate in the designated wells and incubated at room temperature for 1.5 h with shaking. After calibrator and sample incubation, the plates were washed, and 50 μl of TURBO-BOOST antibody solution was added, and the plates were incubated at room temperature for 1 h. After TURBO-BOOST antibody incubation, the plates were washed and 50 μl of enhance solution was added, and the plates were incubated at room temperature for 30 minutes. After enhance solution incubation, the plates were washed, 50 μl of TURBO-TAG detection solution was added, and the plates were incubated with shaking at 27 °C for 1 hour. After TURBO-TAG dection incubation, the plates were washed gently, 150 μl of MSD-Gold read buffer A was added, and the plates were immediately read using the MSD instrument. The IFN-α concentrations were determined by the calibration curves established by fitting the signals from the calibrators using a 4-parameter logistic model with a 1/Y² weighting in Prism v.9.4.1.

Olink

Cytokines in plasma were measured using Olink multiplex proximity extension assay (PEA) inflammation panel (Olink proteomics: www.olink.com) according to the manufacturer’s instructions as described before ⁷⁹. The PEA is a dual-recognition immunoassay, in which two matched antibodies labelled with unique DNA oligonucleotides simultaneously bind to a target protein in solution. This brings the two antibodies into proximity, allowing their DNA oligonucleotides to hybridize, serving as template for a DNA polymerase-dependent extension step. This creates a double-stranded DNA ‘barcode’ that is unique for the specific antigen and quantitatively proportional to the initial concentration of target protein. The hybridization and extension are immediately followed by PCR amplification and the amplicon is then finally quantified by microfluidic qPCR using Fluidigm BioMark HD system. Normalized Protein eXpression (NPX) values were used for downstream analysis after initial QC filtering. PCA analysis was conducted with the R package “pcaMethods” ⁸⁰.

CyTOF analysis of PBMC samples

1 million peripheral blood mononuclear cells (PBMCs) were fixed using 2% PFA for 30 min at RT and washed. Fixed PBMCs were then resuspended in freezing media (10% DMSO + 90% FBS), transferred to cryovials, and stored at −80°C until read for downstream processing and staining. As previously described ²¹, fixed frozen PBMCs were thawed in a 37°C water bath and gently resuspended using 1 mL CSM (PBS supplemented with 2% BSA, 2mM EDTA and 0.1 % sodium azide) and transferred to 15 mL conical tubes containing 9 mL CSM. Samples were washed twice using CSM and counted. Cells were then permeabilized and barcoded using the Cell-ID^TM 20-Plex Pd Barcoding Kit (Fluidigm, catalog # 201060). Post-barcoding, cells were washed, pooled into one barcode composite, and counted. The pooled composite was then stained for 30 min at RT with a pre-titrated surface antibody cocktail (DataS3). After surface staining, cells were washed twice with CSM and fixed in 4% PFA (freshly prepared paraformaldehyde in PBS) for 10 min at RT. Cells were then washed with CSM and permeabilized with 100% cold MeOH (Sigma), and kept overnight at −80°C. The next day, cells were washed with CSM and counted. Intracellular staining was then performed for 30 min at RT with a pre-titrated intracellular antibody cocktail (DataS3), followed by two CSM washes. Finally, cells were stained with iridium-containing DNA intercalator (Fluidigm) for 20 min at RT, washed first with CSM and then with MilliQ water. The washed cells were resuspended in MilliQ water supplemented with 1x EQ four element calibration beads (Fluidigm) and acquired on Helios mass cytometer (Fluidigm). The raw FCS files were normalized and concatenated using the Fluidigm software. The normalized .fcs files were then processed in FlowJo software v10 (BD Biosciences) for debarcoding. Briefly, the normalized .fcs file was used to gate single cells based on DNA content and event length in FlowJo. The single cells were then reimported and debarcoded using Helios software version 7.0.5189. The debarcoded samples were analyzed using FlowJo or R version 1.2.1335 for downstream analysis and visualization.

CyTOF Analysis

High-dimensional analysis of phospho-CyTOF data was performed using a previously described R-based pipeline ⁸¹. In brief, the raw .fcs files were imported into R, and the data were transformed to normalize marker intensities using arcsinh with a cofactor of 5. For visualization, another transformation was applied that scales the expression of all values between 0 and 1 using percentiles as the boundary. Cell clustering was performed with 4,000 cells randomly selected from each sample using FlowSOM ⁸² and ConsensusClusterPlus ⁸³. The transformed matrix was used as an input for FlowSOM, and cells were separated into 20 clusters. To obtain reproducible results (avoid random start), a seed was set for each clustering. The 20 clusters were manually annotated based on the lineage marker expression and were merged to produce the final clusters. The clusters were visualized in two-dimensional space using UMAP ⁸⁴. In parallel, the data were manually gated to identify 34 immune cell subpopulations that were not well-distinguished in UMAP and used for all quantification purposes.

Single-cell multi-omics experiments

Single-cell ATAC and RNA-seq libraries were prepared using the Chromium Single Cell Multiome ATAC + Gene Expression platform (10X Genomics, Pleasanton, CA). Briefly, cryopreserved PBMCs were thawed and processed for single nuclei multi-omics analysis according to the manufacturer’s instructions (10x Genomics, CG000365 Rev B). Nuclei were obtained by incubating PBMCs for 3.10 minutes in freshly prepared Lysis buffer and washed and resuspended in chilled diluted nuclei buffer (10x Genomics, 2000153). About 9,000 cells were targeted for each experiment. Prepared nuclei were subjected to transposition of open chromatin regions. Next, transposed nuclei, reverse transcription Master Mix, barcoded Gel Beads, and Partitioning Oil were partitioned into single-cell GEMs (Gel Bead in EMulsions) using the 10X Chromium Controller and Next GEM Chip J. Within each GEM, Gel Beads are dissolved and poly-adenlyated (poly-A) mRNA transcripts are captured by uniquely barcoded poly(dt)VN oligos. Simultaneously, accessible chromatin fragments are captured by a separate oligo containing a spacer, unique barcode, and Illumina P5 adaptor sequence. The GEMs were then incubated in a C1000 Touch Thermal Cycler (BioRad) to produce barcodedDNA from the transposed DNA, and full-length cDNA. GEMs are then broken to release and pool single-cell fractions within each sample. Pooled fractions are purified using silane magnetic beads and subjected to PCR amplification to generate sufficient mass for library construction. Next, P7 and a sample index are added to transposed DNA via PCR to generate ATAC libraries. Finally, cDNA is enzymatically fragmented before the addition of P5, P7, i7 and i5 indexes, and TruSeq Read 2 via End Repair, A-tailing, Adaptor Ligation, and PCR resulting in gene expression libraries. Quantitation of gene expression and ATAC libraries was performed using Bioanalyzer High Sensitivity DNA Analysis (Agilent). Libraries were combined into gene expression and ATAC pools and sequenced on an Illumina NovaSeq 6000 system using the read lengths recommended by 10X: 28bp (read 1), 90bp (read 2), 10bp (i7 index), and 10bp (i5 index) for gene expression libraries and 50bp (read 1), 49bp (read 2), 8bp (i7 index) and 24bp (i5 index) for ATAC libraries. CellRanger v.3.1.0 (10xGenomics) was used to demultiplex raw sequencing data and quantify transcript levels against the 10x Genomics GRCh38 reference v.3.0.0.

Single-cell gene expression analysis

The single cell RNA-seq data was processed with Seurat (v4.0.5) ⁸⁵. We removed cells with less than 800 or greater than 6,000 detected genes, less than 1,000 or great than 60,000 mRNA reads, or greater than 15% mitochondrial reads. Normalization and feature selection were performed using sctransform ⁸⁶. Clusters were identified with Seurat SNN graph construction on Harmony ⁸⁷-corrected PCA embeddings followed by Louvain community detection algorithm. After identification of major cell types, the clustering process was repeated on each cell type separately to get refined clusters, e.g., monocyte sub-clustering, and to remove doublets that are not identifiable in the first-round clustering of all the cells. Differentially expressed (DE) genes were identified using Wilcoxon Rank Sum test in Seurat. For comparison between two groups (e.g., group A vs. group B), two modes of DE analyses were performed: 1) DEall: all the samples in group A vs. all the samples in group B; 2) DEeach: each sample in group A vs. all the samples in group B. For each mode of DE analysis, genes with p-value ≤ 0.01 were ranked by log2 fold change and used as input in gene-set enrichment analysis (GSEA) analysis implemented in the fgsea R package ⁸⁸. Enrichment was assessed with gene lists in Blood Transcriptomic Modules ⁸⁹. For each comparison between two groups, enriched gene sets were filtered according to the following criteria: 1) the gene set was enriched in DEall comparison (adjust p value ≤ 0.05) 2) the gene set was enriched in at least two DEeach comparisons (adjust p value ≤ 0.05). Only gene sets satisfying both criteria were kept. Selected gene sets shown in the main figures were manually curated to select gene sets relevant to immunology and often enriched in several cell types across multiple DE comparisons.

The ISG score calculation

The ISG score was calculated as the geometric mean of 33 top differentially expressed genes between acute infection and healthy controls (adjusted p value < 0.05, average log2 fold change >= 0.5) across all cell types in 8 BTM terms that are involved in antiviral interferon response (M165, M75, M150, M127, M67, M68, M111.1, M111.0).

Monocyte single-cell gene expression integration

The single-cell gene expression data of monocytes from the previous adult COVID-19 infections study ⁹⁰ and the Pfizer vaccine study ⁹¹ were integrated with monocytes in the current study using the Seurat integration workflow with reciprocal PCA algorithm. The normalized expression data in the integrated assay was then used to calculate the distance between clusters and the ISG score in the integrated data. The Euclidean distance between clusters was calculated using the average expression of the top 20 marker genes in each cluster. The folder change between C14.1, adults_cov2_c11, and adults_vacc_c8 clusters was calculated using limma ⁹².

Single-cell chromatin accessibility analysis

The single-cell ATAC data was processed with ArchR (v1.0.1) ⁹³. The cell type annotations were transferred from the single-cell RNA-seq analysis. Transcription factor (TF) binding motifs were annotated using JASPAR 2016 transcription factor binding database ⁹⁴. The per-cell motif deviations and scores were calculated using chromVAR ⁹⁵. The difference in motif deviations between two groups was tested using the Wilcoxon Rank Sum test. Adjusted p values were computed using Bonferroni correction. Similar to DE gene analysis, differential accessibility motifs between two groups (e.g., group A and group B) are defined according to the following criteria: 1) the TF motif deviations are significantly different between all cells in group A and group B (adjust p value ≤ 0.001 for comparisons between acute infection and healthy controls, and 0.0001 for comparisons between convalescent and healthy controls). 2) when comparing cells in each sample in group A vs. all cells in group B, the TF motif deviations are significantly different in at least three comparisons between acute infection and healthy controls, and at least two comparisons between convalescent and healthy controls (adjust p value ≤ 0.05). The differential gene score analysis was performed using the “getMarkerFeatures” function in ArchR with Wilcoxon Rank Sum test. Genes with FDR ≤ 0.1 were used as input in gene set over-representation analysis for enrichment of Blood Transcriptomic Modules with hypergeometric test. The tracks of antiviral regions were generated using the “plotBrowserTrack” function in ArchR.

Correlation network

The correlation networks were computed using four data types: 1) transcriptomics, in which the average scores of BTM clusters in each cell type in each sample were used as features; 2) epigenomics, in which the average deviations of TF motif clusters in each cell type in each sample were used as features; 3) proteomics/CyTOF, in which the average level of proteins in each cell type in each sample were used as features; 4) Plasma cytokine levels in monocytes and dendric cells; 5) IFNa2 levels. The BTM score in each cell was calculated as the geometric mean of the expression of all the measured genes in the BTM gene sets. The BTM scores were then aggregated by each sample and each cell type, and hierarchical clustering (Euclidean distance, ward.D2 agglomeration method) was applied to identify BTM clusters. The TF motif clusters were identified in a similar way. Spearman’s rho correlation was used in all the correlation calculations. The correlation networks were identified and visualized using igraph ⁹⁶ implemented in the ggraph package ⁹⁷. Fruchterman-Reingold layout was used.

Bulk RNA-seq

Blood was collected into Tempus Blood RNA tubes (Applied Biosystems) and the RNA was extracted using the MagMAX for Stabilized Blood Tubes RNA Isolation Kit, compatible with Tempus Blood RNA Tubes (ThermoFisher Scientific). RNA quality was assessed using a TapeStation 4200 (Agilent) and then 200 nanograms of total RNA was used as input for cDNA synthesis and library preparation using the KAPA RNA HyperPrep kit with RiboErase (HMR) Globin (Roche) according to the manufacturer’s instructions. Libraries were validated by capillary electrophoresis on a TapeStation 4200 (Agilent), pooled at equimolar concentrations, and sequenced with PE100 reads on an Illumina NovaSeq 6000, yielding ∼60 million reads per sample on average.

Bulk RNA-seq analysis

Alignment was performed using STAR version 2.7.3a ⁹⁸ and transcripts were annotated using a composite genome reference which included GRCh38 Ensembl release 100 and SARS-CoV-2 (GCF_009858895.2,ASM985889v3). Transcript abundance estimates were calculated internal to the STAR aligner using the algorithm of htseq-count ⁹⁹. ENSEMBL IDs were filtered to remove low/non-expressed transcripts (<5 reads in >50% of samples). Gene-level counts were created by averaging counts from all ENSEMBL IDs mapping to the same gene symbol (IDs mapping to multiple symbols were discarded), using the bioMart package. Gene counts were normalized using the estimateSizeFactors function of DESeq2. To evaluate monocyte signatures expression during infection, scores for each signature were computed as the average of all genes within the signature. All analysis was performed in R 4.1.1. if not stated differently. Figures were generated using ggplot2 ¹⁰⁰, and Complexheatmap ¹⁰¹.

Supplementary Figure legends

• Download figure
• Open in new tab

Supp Figure 1.

Humoral response to COVID-19 infection in infants, related to Figure 1 A) Summary of the conducted assays. B,C) Kinetics of developing binding (B) and neutralizing (C) antibody response in infants with COVID-19. D,E) Comparison of binding (D, Non-Omi n = 13, Omi n = 18, Conv n = 30) and neutralizing (E, Non-O n = 13, Omi n = 18) titers between different variants in infants. E) Kinetics of specific autoantibodies which increased during or after COVID-19. Initial statistical comparisons were conducted in a paired fashion with the Wilcoxon signed rank test (Pre vs. Conv, n = 27) and Benjamini-Hochberg correction. Validation tests (M-ctrl vs. Conv) were conducted using Wilcoxon rank sum test (conv n = 30, m-ctrl n = 27).

• Download figure
• Open in new tab

Supp Figure 2.

Prevalence of autoantibodies in infants with COVID-19, related to Figure 1 Heatmap depicting plasma IgG antibodies against the indicated autoantigens and cytokines and chemokines. Prototypes (n=15), adult controls (n=10), adults with acute COVID-19 infection (n=15), infant controls who have not had COVID-19 (n=27), and longitudinal samples from infants who had COVID-19 (pre-infection (n=27), acute infection (n=19), acute omicron infection (n=18), and convalescent samples (n=30)) are shown. Prototypes are positive control samples from patients with known autoimmune disorders. Colors indicate autoantibodies with MFI >5 SD (red) or <5 SD (black) above average for healthy infants. MFIs <3000 were excluded.

• Download figure
• Open in new tab

Supp Figure 3.

Memory B and T cell response, related toFigure 2 A) Gating strategy used for SARS-CoV-2 spike specific IgG+ memory B cell staining and single-cell sorting. Gating was on singlets that were CD20+ and CD3-CD14-IgM-IgD-CD27low/+ IgG+. Sorted cells were Wuhan spike-AlexaFluor 488+ and/or Omicron spike-BV421+. B) The percentage of SARS-CoV-2 spike-specific IgG+ memory B cells in healthy, acute, and convalescent infant individuals. The sample number for each group is indicated in brackets. C) As in (B), the percentage of SARS-CoV-2 spike-specific IgG+ memory B cells in adult individuals with mild, severe, and ICU symptoms and in adult convalescent individuals. The sample number for each group is indicated in brackets. D-F) T cells were stimulated with overlapping peptides against WT (D-F) and Omicron (F) variants. Cytokine production was determined via flow cytometry. D) Box plot showing the fraction of responding T cells at different infection stages. E) Box plot showing the fraction of multifunctional T cells at different infection stages. F) Comparison of the multifunctional T cell response after stimulation with WT and Omicron (Om) peptides. Statistical comparisons were conducted with Wilcoxon rank sum test. Solid line indicates median healthy response; dashed line indicates 3x median healthy response.

• Download figure
• Open in new tab

Supp Figure 4.

Cytokine response to COVID-19 infection in infants, related to Figure 3 A) Comparison of key inflammatory mediators during COVID-19 infection in infants and adults, stratified by severity. Infant infection was overall mild or asymptomatic. B,C) Comparison of plasma IFN⍺2 levels between Non-Omicron (Non-O, n = 19) and Omicron (Omi, n = 22) infected infants (B), and infants (n = 41) and adults (n = 15) (C). Statistical comparisons were conducted with the unpaired t-test (A) and Wilcoxon rank sum test (B,C).

• Download figure
• Open in new tab

Supp Figure 5.

Cellular immune response to COVID-19 infection, related to Figure 4 A) Heatmap showing the expression of CyTOF markers in all manually gated subsets. B) Correlation analysis between the frequency of indicated cell types (y-axis) and plasma IFN⍺2 levels (x-axis). C) Heatmap showing the expression levels of significantly changed markers in healthy and infected samples. Colors indicate the infection stage. D) Boxplots showing the frequency of pDCs as a proportion of total CD45+ cells (left) and the expression of pS6 in pDCs for different infection stages. E) Boxplots showing HLA-DR expression on mDCs for different infection stages. Correlation analyses were conducted using Spearman correlation.

• Download figure
• Open in new tab

Supp Figure 6.

Single-cell multi-omics analysis of immunity to COVID-19 infection in infants, related to Figure 5 A) Full ring plot from Figure 5c. Pairwise comparison of genes from healthy (n = 16) and COVID-19–infected infants at different times during acute infection (D0-5: n = 5, D5-10: n=7, D10+: n=6) was conducted for each cluster. DEGs were analyzed for the enrichment of BTMs. Ring plot shows an abridged representation of enriched pathways in each cluster. Size indicates the number of samples with enrichment; colors indicate the normalized enrichment score. B) UMAP representation of the integrated analysis of monocyte clusters from this study and from adult COVID-19 patients ²¹ and adult subjects immunized with the COVID-19 vaccine ³⁰. Colors indicate the study origin of cells (top), the cell cluster (middle), and the expression of ISGs (bottom). C) DEGs determined between infant C14.1 and adult COVID-19-infection C11 monocyte clusters are plotted and ranked by fold change. C) DEGs determined between infant C14.1 and adult vaccination C8 monocyte clusters are plotted and ranked by fold change. D) Correlation analysis between plasma IFN⍺2 levels and average ISG levels in each cell type. E) Heatmap showing expression of NK cell activation genes enriched in NK cells in Figure 5c (green box).

• Download figure
• Open in new tab

Supp Figure 7.

Single-cell epigenomic analysis of immunity to COVID-19 infection in infants, related to Figure 6 A) Pairwise comparison of TF motif accessibility between convalescent and pre or matched-ctrl samples was conducted for each cluster. Color indicates difference in TF accessibility; non-significant changes (FDR>=0.0001 or changed in less than two subjects) are grey. Size indicates the number of samples with significant change. (pre: n = 9, matched-ctrl: n = 7, Conv: n = 9) B) Box plot showing sample-level accessibility of selected TFs from (A) in CD14+ monocytes. C) Gene tracks showing chromatin accessibility at indicated time points in CD16+ monocytes.