Single-cell analysis of the epigenomic and transcriptional landscape of innate immunity to seasonal and adjuvanted pandemic influenza vaccination in humans

Lead Contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Bali Pulendran (bpulend{at}stanford.edu).

Materials Availability

This study did not generate new unique reagents.

Data and Code Availability

EpiTOF data is available at: will be available upon publication. Bulk ATAC-seq and RNA-seq data is available here: will be available upon publication. ScATAC-seq and scRNA-seq data are deposited at: will be available upon publication.

H5N1/H5N1+AS03

This study was conducted in Atlanta, GA. Subjects were males and non-pregnant females who met the eligibility criteria as listed on clinicaltrials.gov (NCT01910519). We enrolled a total of 50 healthy adults who were randomized into two groups receiving either the adjuvanted (H5N1+AS03, n=34) or unadjuvanted (H5N1, n=16) GSK avian influenza vaccine. While both vaccines contained split-virion (A/Indonesia/5/2005) inactivated hemagglutinin antigen, the adjuvanted vaccine additionally contained the AS03 adjuvant system (containing DL-alpha-tocopherol and squalene in an oil-in-water emulsion). Subject demographics are listed in (DataS1). Written informed consent was obtained from each subject and protocols were approved by Institutional Review Boards of Emory University.

In-vitro stimulation and intracellular flow cytometry experiments

Samples from healthy subjects were collected at Stanford Blood Center or derived from the before-vaccination time point of a previous vaccination trial (Nakaya et al., 2015). All subjects provided a confidential medical history card and completed informed consent to donate blood for clinical or research uses. We exclude subjects with known diseases, including but not limited to HIV, and hepatitis infections. Purification of buffy coat or LRS chamber from whole blood was performed at Stanford Blood Center to enrich for leukocytes prior to PBMC isolation.

From the vaccination trial (Nakaya et al., 2015), only samples from subjects aged 26 – 41 were selected for this paper. Samples were only selected from the before vaccination time point at day 0. Written informed consent was obtained from each subject with institutional review and approval from the Emory University Institutional Review Board.

Cells, plasma and RNA isolation

Peripheral blood mononuclear cells (PBMCs) and plasma were isolated from fresh blood (CPTs; Vacutainer with Sodium Citrate; BD), following the manufacturer’s protocol. For samples from Stanford Blood Center, PBMCs isolated from whole blood, buffy coat or LRS chamber by Ficoll density gradient centrifugation using Ficoll-Paque PLUS (GE Healthcare, #17-1440-02). PBMCs were frozen in DMSO with 10% FBS and stored at –80C and then transferred on the next day to liquid nitrogen freezers (–196C). Plasma samples from CPTs were stored at –80C. Trizol (Invitrogen) was used to lyse fresh PBMCs (1 mL of Trizol to ∼1.5×10^6 cells) and to protect RNA from degradation. Trizol samples were stored at –80C.

Mass cytometry sample processing, staining, barcoding and data acquisition

Cryopreserved PBMCs were thawed and incubated in RPMI 1640 media (ThermoFisher) containing 10% FBS (ATCC) at 37°C for 1 hour prior to processing. Cisplatin (ENZO Life Sciences) was added to 10 µM final concentration for viability staining for 5 minutes before quenching with CyTOF Buffer (PBS (ThermoFisher) with 1% BSA (Sigma), 2mM EDTA (Fisher), 0.05% sodium azide). Cells were centrifuged at 400g for 8 minutes and stained with lanthanide-labeled antibodies (DataS2) against immunophenotypic markers in CyTOF buffer containing Fc receptor blocker (BioLegend) for 30 minutes at room temperature (RT). Following extracellular marker staining, cells were washed 3 times with CyTOF buffer and fixed in 1.6% PFA (Electron Microscopy Sciences) at 1×106 cells/ml for 15 minutes at RT. Cells were centrifuged at 600g for 5 minutes post-fixation and permeabilized with 1 ml ice-cold methanol (Fisher Scientific) for 20 minutes at 4°C. 4 ml of CyTOF buffer was added to stop permeabilization followed by 2 PBS washes. Mass-tag sample barcoding was performed following the manufacturer’s protocol (Fluidigm). Individual samples were then combined and stained with intracellular antibodies in CyTOF buffer containing Fc receptor blocker (BioLegend) overnight at 4°C. The following day, cells were washed twice in CyTOF buffer and stained with 250 nM 191/193Ir DNA intercalator (Fluidigm) in PBS with 1.6% PFA for 30 minutes at RT. Cells were washed twice with CyTOF buffer and once with double-deionized water (ddH2O) (ThermoFisher) followed by filtering through 35μm strainer to remove aggregates. Cells were resuspended in ddH2O containing four element calibration beads (Fluidigm) and analyzed on CyTOF2 (Fluidigm).

Bulk stimulation experiment

Aliquots of thawed PBMCs from the EpiTOF experiment described above were washed and resuspended in RPMI 1640 (Corning, 10-040-CV) containing 10% FBS (Corning, 35-011-CV) and 1x Antibiotics/Antimycotics (Lonza, 17-602E) [complete media abx] at 4×10^6 cells/mL. 100 µL of cell solution were added to each well of a 96-well round-bottomed tissue culture plate and mixed with 100 µL of either complete media abx (unstim), a cocktail of synthetic TLR ligands mimicking bacterial pathogens (bac: 0.025 µg/mL LPS, 0.3 µg/mL Flagellin, 10 µg/mL Pam3CSK4), or a cocktail of synthetic TLR ligands mimicking viral pathogens (vir: 4 µg/mL R848, 25 µg/mL pI:C). Depending on cell numbers, PBMCs from each sample were stimulated with all 3 conditions in duplicate. After 24h of incubation at 37C and 5% CO2, cells were spun down, supernatant was carefully transferred into new plates, and immediately frozen at -80C until further analysis using Luminex.

Luminex TIV

The Luminex assay was performed by the Human Immune Monitoring Center, Stanford University School of Medicine. Human 62-plex custom Procarta Plex Kits (Thermo Fisher Scientific) were used according to the manufacturer’s recommendations with modifications as follows: Briefly, Antibody-linked magnetic microbeads were added to a 96-well plate along with custom Assay Control microbeads (Assay Chex) by Radix Biosolutions. The plates were washed in a BioTek ELx405 magnetic washer (BioTek Instruments). Neat Cell culture supernatants (25ul) and assay buffer (25ul) were added to the 96 well plate containing the Antibody-coupled magnetic microbeads, and incubated at room temperature for 1 h, followed by overnight incubation at 4°C. Room temperature and 4°C incubation steps were performed on an orbital shaker at 500–600 rpm. Following the overnight incubation, plates were washed in a BioTek ELx405 washer (BioTek Instruments) and then kit-supplied biotinylated detection Ab mix was added and incubated for 60 min at room temperature. Each plate was washed as above, and kit-supplied streptavidin–PE was added. After incubation for 30 min at room temperature, wash was performed as described, and kit Reading Buffer was added to the wells. Each sample was measured in two technical replicates where cell numbers allowed. Plates were read using a FlexMap 3D Instrument (Luminex Corporation). Wells with a bead count <50 were flagged, and data with a bead count <20 were excluded.

Luminex H5N1/H5N1+AS03

This assay was performed by the Human Immune Monitoring Center at Stanford University. A custom 41 plex from EMD Millipore kits was assembled and included: 1. A Pre-mixed 38 plex Milliplex Human Cytokine/Chemokine kit (CAT# HCYTMAG-60K-PX38) 2. ENA78/CXCL5 (CAT# HCYP2MAG-62K-01) 3. IL-22 (CAT# HTH17MAG-14K-01). 4. IL-18 (HIL18MAG-66K). Manufacturer’s recommendations were followed with modifications described. Briefly: neat supernatant samples (25ul) were mixed with antibody-linked magnetic beads in a 96-well plate containing assay buffer, for an overnight incubation at 4°C. Cold and Room temperature incubation steps were performed on an orbital shaker at 500-600 rpm. Plates were washed twice with wash buffer in a BioTek ELx405 washer (BioTek Instruments). Following one-hour incubation at room temperature with biotinylated detection antibody, streptavidin-PE was added for 30 minutes. Plates were washed as above, and PBS was added to wells for reading in the Luminex FlexMap3D Instrument with a lower bound of 50 beads per sample per cytokine. Each sample was measured in duplicate wells where cell numbers allowed. Custom Assay Chex control beads were added to all well (Radix Biosolutions). Wells with a bead count <50 were flagged, and data with a bead count <20 were excluded.

H3K27ac antibody conjugation

α-H3K27ac antibody was labeled using the Lightning-Link Rapid DyLight 488 Antibody Labeling Kit according to manufacturer’s instructions (Novus Biologicals, 322-0010). In brief, 100 µg of antibody was mixed with 10 µL of LL-Rapid modifier reagent and added onto the lyophilized dye. After mixing, solution was incubated at room temperature overnight in the dark. The next morning, 10 µL of LL-Rapid quencher reagent was added.

In-vitro stimulation and intracellular flow cytometry experiments

Cryopreserved PBMCs were thawed, counted, and resuspended in RPMI 1640 (Corning, 10-040-CV) supplemented with 10% FBS (Corning, 35-011-CV) [complete media] at a concentration of 4×10^6 cells/mL. Next, 150µL of cell suspension (6×10^5 cells) was added to each well of a 96-well round-bottomed tissue culture plate and mixed with 50 µL of inhibitor solution containing either Trichostatin A (TSA; CST, 9950S), A-485 (Tocris, 6387), or Cl-Amidine (EMD Millipore, 506282) in complete media. After 2h of incubation at 37C and 5% CO2, the cells were stimulated by adding either LPS (0.025 µg/mL; Invivogen, tlrl-pb5lps) or R848 (4 µg/mL; Enzo Life Sciences, ALX-420-038-M005) to the cultures. After another 2h of incubation, Brefeldin A (10 µg/mL; Sigma Aldrich, B7651-5MG) was added to all cultures and cells were incubated for a final 4h.nAfter a total of 8h of incubation, cells were washed twice with 150 µL PBS (GE Life Sciences, SH30256.LS) and stained for viability using 100 µL of Zombie UV Fixable Viability Dye in PBS (1:1000; Biolegend, 423108). After incubating for 30 minutes at 4C in the dark, cells were washed twice with 150 µL PBS and blocked with 100 µL of PBS supplemented with 5% FBS, EDTA (2 mM; Corning, 46-034-cl), and human IgG (5 mg/mL; Sigma Aldrich, G4386-5G) [blocking buffer] for 15 minutes at 4C in the dark. After incubation, cells were stained for surface markers with 100 µL of antibody cocktail containing α-CD14 BUV805, α-CD3, CD19, CD20 BUV737, α-CD123 BUV395, α-HLA-DR BV785, α-CD16 BV605, α-CD56 PE-CY7, α-CD11c APC-eFluor780 in blocking buffer for 20 minutes at 4C in the dark. Next, cells were washed twice with 150 µL PBS, and fixed in 200 µL eBioscience Foxp3 Fixation/Permeabilization solution (ThermoFisher Scientific, 00-5523-00) for 30 minutes at 4C in the dark. Afterwards, cells were washed twice with 100 µL eBioscience Foxp3 permeabilization buffer and blocked with 100 µL permeabilization buffer containing human IgG (5 mg/mL) overnight at 4C in the dark. Cells were washed and stained for intracellular markers with 25 µL of antibody cocktail containing α-IL-1b Pacific Blue, α-H3K27ac DyLight 488, α-TNFa PE-Dazzle, α-p-c-Jun PE, and α-H3 AF647 in permeabilization buffer containing human IgG (5 mg/mL) for 60 minutes at 4C in the dark. Finally, cells were washed twice with 150 µL of permeabilization buffer, resuspended in 100 µL PBS containing 0.5% FBS and 2 mM EDTA [FACS buffer], and acquired using a BD FACSymphony flow cytometer. Data was analyzed using Flowjo X software (BD). Briefly, cells were identified via FSC/SSC, doublets were discarded via SSC-A/SSC-H and FSC-A/FSC-H gates, and dead cells were discarded as Zombie UV Fixable Viability Dye high. Monocytes were then identified as CD3^-CD19^-CD20^- and CD14⁺.

FACS sorting – bulk ATAC-seq/RNA-seq

Cryopreserved PBMCs were thawed, washed, counted, and resuspended in PBS (GE Life Sciences, SH30256.LS). 5-10×10^6 cells were washed once more with 2 mL of PBS and stained for viability using 500 µL of Zombie UV Fixable Viability Dye in PBS (1:1000; Biolegend, 423108). After incubating for 30 minutes at 4C in the dark, cells were washed with 2 mL of PBS and resuspended in 500 µL blocking buffer. After spinning cells down, supernatant was discarded, and cells were resuspended in 50 µL antibody cocktail containing α-CD3, CD19, CD20 BUV737, α-CD123 BUV395, α-HLA-DR BV785, α-CD14 BV605, α-CD56 BV510, α-CD1c BV421, α-CD327 AF488, α-CD370 PE, α-CD11c APC-eFluor780, α-CD15 AF700, and α-Axl APC in blocking buffer. Cells were stained for 15 minutes at 4C in the dark. Finally, cells were washed with 2 mL of FACS buffer, resuspended in PBS containing 5% FBS at 10-20×10^6 cells/mL, and stored at 4C before sorting on a FACS Aria Fusion (BD). During sort, live innate cells were identified by gating on Viability Dye^- CD3^-CD19^-CD20^- cells. Within this population, CD14⁺ monocytes were identified as CD14⁺, mDCs were identified as CD14^- CD56^-HLA-DR⁺CD16^-CD11c⁺CD123^-, and pDCs were identified as CD14^-CD56^-HLA-DR⁺CD16^-CD11c^-CD123⁺.

Omni ATAC-seq of purified immune cells

Atac was performed on purified innate immune cell subsets immediately after sorting based on the low-input Omni-Atac protocol described before (Corces et al., 2017). In brief, 1,500 – 5,500 cells were washed with ATAC resuspension buffer (10 mM Tris-HCl pH 7.5 [Invitrogen, 15567027], 10 mM NaCl [Invitrogen, AM9760G], 3 mM MgCl2 [Invitrogen, AM9530G], in water [Invitrogen, 10977015]) and supernatant was carefully aspirated, first using a P1000, then a P200 pipette. Next, 10 µL transposition mix (0.5 µL Tn5, 0.1 µL 10% Tween-20, 0.1 µL 1% Digitonin, 3.3 µL PBS, 1 µL water, and 5 µL tagmentation buffer) was added to the pellet and cells were resuspended by pipetting up and down 6 times. Tagmentation buffer was prepared locally by resuspending 20 mM Tris-HCl pH 7.5, 10 mM MgCl2, and 20% Dimethyl Formamide (Sigma Aldrich, D4551-250ML) in water. Cells were incubated at 37C for 30 minutes under constant mixing. After tagmentation, the reaction was cleaned up using the MinElute PCR Purification Kit (Qiagen, 28006) according to manufacturer’s instructions. Cleaned DNA was eluted in 21 µL of elution buffer, stored at -20C, and shipped to Active Motif for sequencing library preparation. At Active Motif, tagmented DNA was amplified with 10 cycles of PCR using customized Nextera PCR Primers 1 and 2 (see Key Resource table), and purified using Agencourt AMPure SPRI beads (Beckman Coulter, A63882). Resulting material was quantified using the KAPA Library Quantification Kit for Illumina platforms (Roche, 07960255001), and sequenced with PE42 sequencing on the NextSeq 500 sequencer (Illumina).

Bulk RNA-seq of purified immune cells

Bulk RNA-seq was performed on purified CD14⁺ monocytes after sorting. In brief, after sorting, 5,500 cells were washed, resuspended in 350 µL chilled Buffer RLT (Qiagen, 79216) supplemented with 1% beta-Mercaptoethanol (Sigma, M3148-25ML), vortexed for 1 minute, and immediately frozen at -80C. RNA was isolated using the RNeasy Micro kit (Qiagen, 74004) with on-column DNase digestion. RNA quality was assessed using an Agilent Bioanalyzer and total RNA was used as input for cDNA synthesis using the Clontech SMART-Seq v4 Ultra Low Input RNA kit (Takara Bio, 634894) according to the manufacturer’s instructions. Amplified cDNA was fragmented and appended with dual-indexed bar codes using the NexteraXT DNA Library Preparation kit (Illumina, FC-131-1096). Libraries were validated by capillary electrophoresis on an Agilent 4200 TapeStation, pooled at equimolar concentrations, and sequenced on an Illumina NovaSeq6000 at 100SR, yielding 20-25 million reads per sample.

FACS sorting – scATAC-seq/RNA-seq

Cryopreserved PBMCs were thawed and innate immune cell subsets were isolated using FACS as described above (FACS sorting – bulk ATAC-seq/RNA-seq). Within the live gated cells, CD14⁺ monocytes were identified as CD14⁺ (fraction A) while a mixture of the remaining monocyte and dendritic cell subsets was identified as CD14^- CD56^-HLA-DR⁺ (fraction B). After sorting and depending on the number of isolated cells, fraction A and B were mixed at a 2:1 ratio to yield a solution of monocytes and dendritic cells enriched for CD14^- cells.

scRNA-seq

FACS-purified cells were resuspended in PBS supplemented with 1% BSA (Miltenyi), and 0.5 U/μL RNase Inhibitor (Sigma Aldrich). About 9,000 cells were targeted for each experiment. Cells were mixed with the reverse transcription mix and subjected to partitioning along with the Chromium gel-beads using the 10X Chromium system to generate the Gel-Bead in Emulsions (GEMs) using the 3’ V3 chemistry (10X Genomics). The RT reaction was conducted in the C1000 touch PCR instrument (BioRad). Barcoded cDNA was extracted from the GEMs by Post-GEM RT-cleanup and amplified for 12 cycles. Amplified cDNA was subjected to 0.6x SPRI beads cleanup (Beckman, B23318). 25% of the amplified cDNA was subjected to enzymatic fragmentation, end-repair, A tailing, adapter ligation and 10X specific sample indexing as per manufacturer’s protocol. Libraries were quantified using Bioanalyzer (Agilent) analysis. Libraries were pooled and sequenced on an NovaSeq 6000 instrument (Illumina) using the recommended sequencing read lengths of 28 bp (Read 1), 8 bp (i7 Index Read), and 91 bp (Read 2).

scATAC-seq

FACS-purified cells were processed for single nuclei ATAC-seq according to the manufacturer’s instructions (10x Genomics, CG000168 Rev D). Briefly, nuclei were obtained by incubating PBMCs for 3.20 minutes in freshly prepared Lysis buffer following manufacturer’s instructions for Low Cell Input Nuclei Isolation (10x Genomics, CG000169 Rev C). Nuclei were washed and resuspended in chilled diluted nuclei buffer (10x Genomics, 2000153). Next, nuclei were subjected to transposition for 1h at 37C on the C1000 touch PCR instrument (BioRad) prior to single nucleus capture on the 10x Chromium instrument. Samples were subjected to post GEM cleanup, sample index PCR, cleanup and library QC prior to sequencing according to the protocol. Samples were pooled, quantified and sequenced on NovaSeq 6000 instrument (Illumina) with at least minimum recommended read depth (25000 read pairs/nucleus).

IFNa SIMOA

Frozen plasma was shipped to Qunaterix and analyzed using the Simoa® IFN-α Advantage Kit (Quanterix, 100860) according to manufacturer’s instructions. In brief, plasma and reagents were thawed at room temperature. Cailbrators, controls, and plasma were transferred to assay plates. Beads were vortexed for 30 seconds and prepared reagents and samples were loaded into a HD-1/HD-X instrument and analyzed with standard settings. All samples were run in duplicate.

Detection of IFNa and IFNg in plasma and cell culture supernatants

Frozen plasma or supernatant was thawed at room temperature and analyzed using the IFNa and IFNg Human ProQuantum Immunoassay Kits according to manufacturer’s instructions. In brief, samples were mixed with assay dilution buffer at a 1:5 or 1:2 ratio and protein standard was serially diluted in assay dilution buffer. Next, Antibody-conjugates A and B were mixed with Antibody-conjugate dilution buffer and added to each well of a 96-well qPCR plate (Bio-Rad, #HSP9601). Next, diluted sample or standard were added to each well and mixtures were incubated for 1h at room temperature in the dark. Finally, Master max and Ligase were mixed and added to each well. QPCR was conducted on a CFX96 Touch Real-Time Detection System (Biorad) using the recommended instrument settings. After measurements were completed, CT values were calculated using a regression model and exported to the ProQuantum Cloud app that accompanied the kit (apps.thermofisher.com/apps/proquantum). ProQuantum Cloud app was then used to construct a standard curve and calculate protein concentrations from CT values.

IP-10 plasma Luminex

Plasma biomarker levels were assayed using a 10-analyte multiplex bead array (fractalkine, IL-12P40, IL-13, IL-1RA, IL-1b, IL-6, IP-10, MCP-1, MIP-1α, TNFβ; Millipore) prepared according to the manufacturer’s recommended protocol and read using a Bio-Plex 200 suspension array reader (Bio-Rad). Data were analyzed using Bio-Plex manager software (Bio-Rad).

Viral infection assay

Dengue virus (DENV-2, Strain Thailand/16681/84) and Zika virus (PRVABC59) were propagated and titrated on Vero cells and stored at -80C until infection. Cryopreserved human PBMCs were thawed, washed, counted, and resuspended in RPMI 1640 (Thermo Fisher, 72400-047) supplemented with 10% FBS (Corning, 35-011-CV), 1mM Sodium pyruvate (Lonza, 13-115E), and 1x Penicillin/Streptomycin (Lonza, 17-602E) at 1.5×10^6 cells/mL. 200 µL of cell solution (3×10^5 cells) was added to each well of a 96-well round-bottomed tissue culture plate and cells were rested in plates for 4h at 37C and 5% CO₂. After resting, PBMCs were infected with DENV-2 or ZIKV at MOI 1. At 0h, 24h, 48h post infection, PBMCs and supernatant were collected for RNA purification and cytokine analysis, respectively. Supernatants were immediately frozen at -20C and stored until analysis. Cells were suspended in RNA lysis buffer and kept at -20C until analysis. RNA was purified using the Purelink RNA kit according to manufacturer recommendations (Thermo Fisher Scientific, #12183052). For viral load detection, quantitative reverse transcription PCR (qRT-PCR) was conducted using Luna universal probe one-step RT-PCR kit (NEB, #E3006) on a CFX96 C1000 Touch Real-Time Detection System with 96-well plates (Bio-Rad, #HSP9601). RNA standards (ATCC, # VR-3229SD, VR-1843DQ) were used to generate standard curves. Viral RNA copies were normalized by cell number. Utilized primers and probes are listed in the Key Resources table.

Detection of IP-10 in culture supernatant

Culture supernatants were thawed at room temperature and analyzed using the IP-10 enzyme-linked immunosorbent assay (R&D Systems, DIP100) according to the manufacturer’s instructions. In brief, samples were thawed at room temperature and mixed with assay dilution buffer at 1:2 ratio. Protein standard was serially diluted in assay dilution buffer. Samples and standards were incubated in plate for 2h at room temperature. Plates were washed and then incubated with human IP-10 conjugate for 2h at room temperature. After wash, substrate solution was added for 30min. Finally, stop solution was added, A450 and A595 were read on a plate reader (Bio-Rad, iMARK). The concentration of IP-10 was determined by the number of A450-A595 based on the standard curve.

EpiTOF analysis

The exported Flowjo data were then normalized following the approach described in (Cheung et al., 2018). In brief, the value of each histone mark was regressed against the total amount of histones, represented by measured values of H3 and H4. For sample level analyses, the values of each histone mark were averaged for each cell type in each sample. Distances of HSC from lymphoid and myeloid epigenetic profiles were obtained by first computing centers of the epigenetic profiles for the two lineages, and then computing Euclidean distances from the centers for each individual HSC. Distances of HSC from epigenetic profiles of specific cell types were similarly obtained by computing Euclidean distances from the centers of the epigenetic profiles for each cell type. Statistical significance of the differences between groups at the sample level was assessed by computing an effect size with Hedges’g formula (Hedges and Olkin, 2014). All p-values were corrected for multiple comparisons with the Benjamini-Hochberg method (Benjamini and Hochberg, 1995). Dimensionality reduction was performed with applying UMAP (McInnes et al., 2020). For single cell analyses, the normalized values were used as input. Correlation between variables was computed using Pearson’s correlation coefficient. All the analyses were performed using the R framework for statistical computing (Version 3.6.3) (R Core Team, 2020).

TIV bulk gene expression analysis Thomas

Processed data and normalized in Bioconductor by RMA (Irizarry et al., 2003), which includes global background adjustment and quantile normalization. Samples from phase1 subjects in the antibiotics and control arm of the study were selected and statistical tests and correlation analyses were performed using MATLAB. Test details and significance cutoffs are reported in figure legends.

Luminex analysis

Statistical analysis was conducted in R (v 4.0.2) (R Core Team, 2020). First, MFI data was log2 transformed and average MFI and CV was calculated from duplicate cultures where available. For samples with CV > 0.25, the duplicate that was closer to the average of all samples of that subject was kept and the other discarded. In case no other sample was available and CV > 0.25, the sample was discarded. Wells without indication of cytokine production were excluded. Statistical tests, correlation analysis, and hierarchical clustering were performed using the R packages stats (v 4.0.2), ggpubr (v 0.4.0) and pheatmap (v 1.0.12). Test details and statistical cutoffs are reported in the figure legends.

Bulk ATAC-seq pre-processing

Analysis of ATAC-seq data was very similar to the analysis of ChIP-Seq data. Reads were aligned using the BWA algorithm (mem mode; default settings; v 0.7.12) (Li and Durbin, 2009). Duplicate reads were removed, only reads mapping as matched pairs and only uniquely mapped reads (mapping quality >= 1) were used for further analysis. Alignments were extended in silico at their 3’-ends to a length of 200 bp and assigned to 32-nt bins along the genome. The resulting histograms (genomic “signal maps”) were stored in bigWig files. Peaks were identified using the MACS algorithm (v 2.1.0) (Zhang et al., 2008) at a cutoff of p-value 1e-7, without control file, and with the –nomodel option. Peaks that were on the ENCODE blacklist of known false ChIP-Seq peaks were removed. Signal maps and peak locations were used as input data to Active Motif’s proprietary analysis program, which creates Excel tables containing detailed information on sample comparison, peak metrics, peak locations and gene annotations. For differential analysis, reads were counted in all merged peak regions (using Subread), and the replicates for each condition were compared using DESeq2 (v 1.24.0) (Love et al., 2014).

Bulk ATAC-seq analysis

Quality control analysis of ATAC-seq data was performed using Rockefeller University workshop on analysis of ATAC-seq data in R and Bioconductor (https://rockefelleruniversity.github.io/RU_ATAC_Workshop.html). Of 57 unique samples processed, 51 passed QC criteria and, on average, we detected more than 42,000 genomic regions and more than 15×10⁶ unique ATAC tags per sample while the average fraction of reads in peaks was larger than 35% (DataS3). Passed samples showed the characteristic fragment length and TSS enrichment distribution (DataS3). DARs were annotated as promoter, distal and trans regulatory peak for a particular gene based on the distance from the middle of the peak to the nearest transcription start site (TSS) using the ChIPpeakAnno package in R (v. 3.24.1). Promoter, distal and trans regulatory peaks were defined as -2000 bp to +500 bp, -10kbp to +10kbp – promoter, and < -10kbp or > +10kbp from TSS, respectively. The hypergeometric distribution-based enrichment analysis was performed to identify the significance of the DARs. Reactome pathways and TF-target relationship using Chip-seq data from ENCODE (both downloaded from https://maayanlab.cloud/chea3/) were used to identify overrepresented pathways and TFs. EnrichmentMap Pipeline Collection (v 1.1.0) (Merico et al., 2010) for CytoScape (v 3.8.2) (Shannon et al., 2003) was used to create the pathway network. Significantly enriched Reactome pathways (p <= 0.05) for each genomic region were used as input. Pathways were clustered and annotated using the AutoAnnotate function within the pipeline. To test for enrichment of TF motifs in DARs, the chromVAR (v 1.8.0) (Schep et al., 2017) and motifmatchr (v 1.8.0) packages were used in R (v 3.6.0) (R Core Team, 2020). In brief, TF motifs were downloaded from the JASPAR2016 core homo sapiens database (Mathelier et al., 2016) and merged regions were annotated for the presence of all TF binding motifs using the matchMotifs (motifmatchr) function with standard settings. Hypergeometric distribution-based enrichment analysis was then performed to identify enrichment of TF motifs in DARs. To determine the relationship between EpiTOF and ATAC-seq data, the Pearson correlation was computed between EpiTOF H3K27ac levels and normalized read counts in each merged peak region. Positively correlated merged peak regions with p-value < = 0.05 were selected for functional annotation. Enrichment analysis was performed as described above.

Bulk RNA-seq of purified immune cells

Alignment was performed using STAR version 2.7.3a (Dobin et al., 2013) and transcripts were annotated using GRCh38 Ensembl release 100. Transcript abundance estimates were calculated internal to the STAR aligner using the algorithm of htseq-count (Anders et al., 2015). DESeq2 version 1.26.0 (Love et al., 2014) was used for differential expression analysis using the Wald test with a paired design formula and using its standard library size normalization.

Analysis of bulk transcriptomics data from previous TIV studies

Processed bulk transcriptomics data from nine independent TIV studies conducted between 2007 and 2012 were obtained from GEO (accessions: GSE47353, GSE59635, GSE29619, GSE74813, GSE59654, GSE59743, GSE74811, GSE29617, GSE74816) (Barrett et al., 2013; Mohanty et al., 2015; Nakaya et al., 2011, 2015; Thakar et al., 2015; Tsang et al., 2014). After removing samples and genes with missing values as well as extraordinary vaccine time points, we selected only samples from subjects matching the same age range as the current study: 18 – 45 years of age. The remaining samples were batch corrected using ComBat from the sva package in R (v 3.36.0) with study as batch, no covariates, and otherwise standard settings. Statistical tests were performed using the R base and ComplexHeatmap (v 2.4.3) packages. Test details and statistical cutoffs are reported in the figure legends.

scATAC analysis

The CellRanger-atac pipeline (v1.1.0) by 10X Genomics was used for alignment (GRCh38 reference genome), de-duplication, and identification of cut sites for each sample. The samples were then combined using the CellRanger-atac aggregation procedure without depth normalization (--normalize=none). The resulting fragment file was read into SnapATAC (Fang et al., 2020). SnapATAC was used to bin the genome (bin size of 5K) and create a cell-by-bin count matrix. Cells were identified as barcodes with at latest 1000 UMIs, and a promotor ratio (defined as: (fragments in promoter regions + 1) / (total fragments + 1)) of at least 0.1, resulting in a total of [state total number of cells in each experiment], as stated in the results section. Bins that mapped to chrY, mitochondrial DNA, or bins that overlap with ENCODE blacklist regions (Amemiya et al., 2019), were removed. The remaining bins were used for dimensionality reduction using Truncated SVD with the irlba R package (Baglama et al., 2019), and the first 50 dimensions were then used for clustering. MACS2 (Zhang et al., 2008) was then used to call peaks within each cluster using recommended parameters for ATACseq data (--nomodel --shift 100 --ext 200 --qval 5e-2 -B --SPMR). The cluster-specific peaks were merged to a single combined set. SnapATAC was then used to map the fragments to the combined peaks set and create a peak-by-cell binary matrix. In the H5N1/H5N1+AS03 dataset, deeply-sequenced libraries were downsampled to an average of 1500 fragments per barcode by randomly removing counts from these samples at a probability p=1500/(mean fragments per cell in the sample). The dimensionality reduction and clustering procedure described above was then repeated on the peak-by-cell matrix. ChromVAR (Schep et al., 2017) was used with default parameters and the JASPAR2016 (Mathelier et al., 2016) motif database to calculate motif accessibility scores and compute differentially accessible motifs in the data. Hotspot was used to identify informative gene modules that explain heterogeneity within the monocyte population (DeTomaso and Yosef, 2020), using the Bernoulli model and the top 2500 regions (ranked by highest autocorrelation z-score) for module calculation. Modules were then identified using the create_modules function, with min_gene_threshold=200. Similar modules were manually identified and merged by taking the average score across modules. Differentially accessible regions were identified using logistic regression with the glm function in R with the design: y ∼ timepoint + donor + log_fragments to control for donor and library size effects. The coefficient corresponding to the time point was then used as the logFC value, and a Wald test was used to get p-values. For numerical stability, we only included peaks that were detected in at least 5% of the cells included in each comparison. All custom scripts for preprocessing, correlation analysis, and differential accessibility analysis are posted in zenodo [doi: 10.5281/zenodo.4446316]. The hypergeometric distribution-based enrichment analysis was performed to identify the significance of the DARs (p <= 0.05 and detected in at least 5% of cells). Reactome pathways database (both downloaded from https://maayanlab.cloud/chea3/) were used to identify overrepresented pathways. Enrichr (Kuleshov et al., 2016) was used to conduct enrichment analysis of genomic regions within Hotspot modules 2, 3. Enrichr was also used to conduct enrichment analysis of DARs containing an IRF1 motif. Briefly, significant DARs (p <= 0.05 and detected in at least 5% of cells) carrying an IRF1 motif, as determined by chromVAR, were selected. Next, gene names with multiple associated DARs were collapsed in case all DARs changed in the same direction or otherwise discared. Subsequently, gene list was submitted to Enrichr for enrichment using the Reactome_2016 database. Similarly, we used Enrichr together with the ChEA_2016 databases to identify TF target genes enriched in genes that were enhanced after booster vaccination with H5N1+AS03 and that overlapped with changes in accessibility at promoter regions.

scRNA analysis

The CellRanger pipeline (v3.1.0) by 10X Genomics was used for alignment (GRCh38 reference genome), demultiplexing, cell-calling, and filtering. The filtered count matrices from each sample were then aggregated using the CellRanger aggregation procedure without depth normalization (--normalize=none). The resulting count matrix was analyzed with scVI (scvi-tools v0.7.1)(Lopez et al., 2018) with default hyperparameters to fit a low-dimensional latent space, using the experiment annotation for each sample as a batch label for batch correction. Visualization, clustering, and exploratory analyses were performed with VISION (v2.1.0)(DeTomaso et al., 2019). Differential expression analysis between time points was performed with edgeR (Robinson et al., 2010) as described in the package documentation, using the exactTest hypothesis testing for each pairwise analysis.

Bulk transcriptomics vax010

Initial data quality was assessed by background level, 3’ labeling bias, and pairwise correlation among samples via the arrayQualityMetrics package in Bioconductor (Kauffmann et al., 2009). CEL files were normalized via RMA (Irizarry et al., 2003), which includes global background adjustment and quantile normalization. Probes mapping to multiple genes were discarded, and the remaining probes were collapsed to gene level by selecting the probe for each gene with the highest mean expression across all subjects. Statistical tests were performed in MATLAB and R.