High-throughput sequencing-based neutralization assay reveals how repeated vaccinations impact titers to recent human H1N1 influenza strains

Our goal was to measure titers against many influenza viruses using the same workflow as a traditional neutralization assay. The fundamental idea behind our approach is to pool many viruses each encoding a different HA, and then read out neutralization of all these viruses simultaneously using next-generation sequencing (Fig. 1). Implementation of this approach required two key technical innovations. First, we incorporated a unique 16-nucleotide barcode sequence within the genomic segment encoding the HA of each strain (13, 14), thereby allowing each HA to be identified by sequencing just this short barcode region (Fig. S1). Second, we developed a spike-in control RNA to enable the conversion of sequencing counts of the viral barcodes to fraction viral infectivity retained upon serum treatment (Fig. 1).

More specifically, we first created a library of influenza viruses encoding barcoded HAs from different natural strains (Fig. 1; Fig. S1). To perform the sequencing-based neutralization assays, we preincubated the library of barcoded viruses with serial dilutions of serum in a 96-well plate (Fig. 1). We then added cells to each well. The non-neutralized viruses infect cells and express high levels of their viral RNA, including the barcoded HA genomic segment. We extracted RNA from the cells and amplified the barcode region of the HA segment using PCR for next-generation sequencing. To quantify the absolute amount of barcoded viral RNA by sequencing, we added a constant amount of a spike-in control containing known barcodes (Fig. S1C), and normalized the viral barcode counts by the counts for the spike-in control in each well (Fig. 1). We quantified the fraction infectivity remaining for each barcoded virus at each serum concentration as the ratio of the normalized counts of that barcode to the normalized counts of that barcode in the no-serum control wells. We then fit neutralization curves to the fraction infectivity values, quantifying the neutralization titer 50% (NT50) as the serum concentration at which 50% of the viral infectivity is neutralized. With this approach, each row or column of the 96-well plate yielded neutralization titers against all the viral variants in the library. Because we designed the libraries to contain multiple distinct barcodes for each different HA, each row or column also yielded multiple replicate titer measurements for each virus. A more detailed description is provided in the methods and a step-by-step protocol is available on protocols.io (DOI: https://doi.org/10.17504/protocols.io.kqdg3xdmpg25/v1). See https://github.com/jbloomlab/seqneut-pipeline for an open-source computational pipeline for analyzing the data.

We aimed to use the sequencing-based assay to measure neutralization titers against the HAs of H1N1 influenza virus strains that circulated in humans between 2020 and 2022. Therefore, we designed a library comprised of the HAs of the five most recent human H1N1 vaccine strains and an additional 31 recent strains that were selected as representatives of subclades within a phylogeny constructed from sequences sampled over a 2-year timeframe, separated from their parent clade by at least one HA1 mutation, and detected at >1% globally (Fig. 2A). The HA ectodomains from the 31 recent strains and the two most recent vaccine strains included in the library vary by 50 amino-acid substitutions at 47 different sites (Fig. 2B), with most of the variable sites located in the globular head of HA, many in regions that are known to be antigenically important (15). Substitutions in all canonical antigenic epitopes are present among different strains within both the 5a.1 clade and 5a.2 clade sequence sets (Fig. S2).

To generate the barcoded library, we cloned the ectodomain of each HA into a unidirectional influenza reverse genetics plasmid (19) that included a 16-nucleotide barcode introduced in a way that does not disrupt viral genome packaging (Fig. S1) (13, 14). We generated plasmids with two to four distinct barcodes for each different HA to provide replicates within the library. Therefore, although our library contained 36 unique HAs, there were 110 different barcodes in the library. We pooled the plasmids encoding replicate barcodes for the HA for each viral strain and then used reverse genetics to separately generate each of the 36 small pools of two to four barcoded viruses encoding each HA. The seven non-HA influenza segments in these viruses were derived from the laboratory-adapted A/WSN/1933 (H1N1) strain (20), a choice that was made both for biosafety reasons and to ensure that all viruses differed only in their HAs and were otherwise identical. The HA genes for all viruses also derived their genome packaging signals, signal peptide, and endodomain from A/WSN/1933, with only the HA ectodomain changed to correspond to the different viral strains (Fig. S1). The use of these HA regions from A/WSN/1933 both simplified construct design and enabled PCR amplification of the HA-embedded viral barcodes from consistent flanking regions, as well as potentially improving viral genome packaging of the HA segment since all other viral genes are derived from HA. Note that although human serum contains antibodies to a variety of influenza proteins (21), typically only antibodies to HA are strongly neutralizing in single-cycle infection assays like the ones used in our experiments (22, 23).

The independently generated stocks of barcoded viruses with each different HA were then pooled so that each HA was represented at approximately equal transcriptional titers. To avoid separately titering dozens of viruses by TCID50 or plaque assays, we used a sequencing-based titering method (Fig. S3). Specifically, we first infected cells with an equal volume mixture of all virus stocks, then used next-generation sequencing of the viral barcodes to determine the relative transcriptional output of each virus stock, and finally re-pooled to equalize these transcriptional titers across HAs such that each strain contributed roughly the same fraction of transcribed barcodes when cells were infected with the re-pooled library in the absence of serum antibodies (Fig. S3). Note that this normalization yields roughly equal titers of transcriptionally active viral particles for each strain; the physical HA content could still differ across strains if different HAs are expressed at different levels or yield differing ratios of physical to transcriptionally active viral particles.

Example neutralization curves from the sequencing-based neutralization assay against the H1N1 vaccine strain from 2021 to 2022 (A/Wisconsin/588/2019) are shown in Fig. 3A for serum samples collected from an individual at 0 and 30 days post-vaccination. Each row of the 96-well plate generates curves like those in Fig. 3A for each of the 36 different viruses (Fig. S4), with each curve generated in replicate because there are multiple barcodes associated with each viral HA (three barcodes in the case of the virus in Fig. 3A). The curves are conceptually identical to those generated by a traditional neutralization assay, with a fraction viral infectivity measurement made for each serial dilution of serum, and a sigmoidal Hill curve fit to the data to estimate the serum dilution at which half of the viral infectivity is retained.

The neutralization curves generated by the sequencing-based assay are highly reproducible. Reproducibility can be assessed at two levels: how similar are measurements for different barcodes associated with the same HA in each row (or column) of a single plate, and how similar are measurements made for the same HA across different plates run on different days. Visual inspection of Fig. 3A shows the consistency of curves both across barcodes within a row of a single plate and across different plates run on different days. A more systematic analysis of additional sera and viruses shows consistently high reproducibility of the NT50s measured both for replicate barcodes within a plate row (Fig. 3B) and the median NT50s measured for the same sera and viruses in different plates run on different days (Fig. 3C).

The neutralization curves and NT50s measured by the sequencing-based neutralization assay are very similar to those measured using a more traditional fluorescence-based one-virus versus one-serum single-cycle neutralization assay (24, 25) (Fig. 4A and B). The sequencing-based neutralization assay NT50s also correlate well with HAI measurements for the 2020–2021 H1N1 vaccine strain, A/Hawaii/70/2019 (Fig. 4C). Note that the HAI titers are not expected to necessarily be identical to the neutralization titers since HAI only assesses the ability of serum antibodies to prevent HA binding to sialic acid on red blood cells, whereas virus neutralization can also be caused by additional mechanisms in addition to blocking of receptor binding (26–28).

We applied the sequencing-based neutralization assay to quantify the magnitude and specificity of the neutralizing H1 antibody response to influenza vaccination in samples from a randomized placebo-controlled trial (29) of Flublok, a recombinant-HA influenza vaccine (30). We selected samples from individuals who were randomly allocated to receive influenza vaccination for the 2021–2022 winter (containing clade 5a.2 A/Wisconsin/588/2019 (H1N1) vaccine strain), stratified by whether they had or had not also received an influenza vaccine in the previous 2020–2021 season (containing clade 5a.1 A/Hawaii/70/2019 (H1N1) vaccine strain) (Table 1; Fig. 5). We selected 15 adult participants in each group (with or without a previous-year vaccination), ranging in age from 20 to 45 years at enrollment. For this study specifically, we chose individuals within each of the groups who had an HAI titer of <10 against the A/Hawaii/70/2019 (H1N1) vaccine strain prior to vaccination in year 1 (2020–2021); this criterion means that the study focuses on individuals with relatively low initial titers. None of the participants were infected with influenza during the study time frame because there was minimal influenza circulation at the study site in Hong Kong between the initial 2020–2021 blood draw and the 2021–2022 blood draws due to COVID-19 mitigation measures (29, 31, 32). We analyzed sera collected immediately before (day 0) and approximately 30 days after the participants were given the 2021–2022 vaccine (Fig. 5). In addition, for a subset of participants, we also analyzed sera collected approximately 182 days after vaccination (Table 1; Fig. 5).

For each study participant, we generated neutralization landscapes that quantified the neutralization titers against all 36 H1 viral strains by the serum collected at each timepoint (see Fig. 6A for example landscapes, and Fig. S5 and S6 for landscapes for all participants). These neutralization landscapes showed a variety of patterns. For instance, vaccination of participant D10066 induced a broad response against all viral strains (Fig. 6). Vaccination of participant D10011 induced a large increase in neutralizing titers against 5a.1 clade viruses, but a smaller increase in titer to most 5a.2 clade viruses other than the vaccine strain A/Wisconsin/588/2019 and a closely related strain (Fig. 6). Participant D10366, who had been vaccinated the previous season, had high pre-vaccination titers against 5a.1 clade viruses but no neutralizing response to vaccination in 2021–2022 (Fig. 6). A variety of other patterns can be seen by examining the neutralization landscapes for all participants in Fig. S5 and S6.

Across the entire study cohort, participants who had been vaccinated the previous year had higher initial titers against recent viruses prior to receipt of the 2021–2022 vaccine compared to those who had not been vaccinated in the previous year (Fig. 7A; Fig. S6). The pre-vaccination titers of these previous-year vaccinees were higher against 5a.1 clade than 5a.2 clade viruses, probably because the H1N1 component of the previous year (2020–2021) vaccine was a 5a.1 clade strain (A/Hawaii/70/2019). The group that received the placebo in the 2020–2021 season had low pre-vaccination titers against all recent strains included in the library (Fig. 7A; Fig. S5), an expected result since study participants were chosen based on low initial baseline titers against the 2020–2021 vaccine strain.

Following vaccination in 2021–2022, most individuals who had not been vaccinated in the previous year responded with large increases in neutralizing antibody titers against all strains (Fig. 7; Fig. S5). By contrast, only some individuals who had also been vaccinated in the previous year showed an appreciable response to vaccination in 2021–2022, and the magnitude of this response was generally smaller than in the group with a previous-year vaccination (Fig. 7; Fig. S5 and S6). Specifically, only 8 of 15 participants in the previous-year vaccination group had at least a fourfold increase in titer against the vaccine strain at 30 days post-vaccination, while 14 of 15 individuals in the group without a previous-year vaccination had at least a fourfold increase (Fig. S5 and S6). Both groups showed some decrease in titer from day 30 to day 182, consistent with other studies showing waning of the antibody response to influenza vaccination within a few months (33, 34).

However, despite the much stronger response to vaccination in the first-time versus previous-year vaccinees, the absolute neutralization titers after vaccination are not dramatically different between the two groups (Fig. 8; Fig. S7). At 30 days post-vaccination, the first-time vaccinees do have higher titers than the previous-year vaccinees (as assessed by the median titers across the groups) to nearly all strains, with some of the differences being statistically significant—but there is substantial overlap between the two groups, and the differences in absolute titers are much more modest than the differences in fold change (Fig. 8; Fig. S7). Many of the strains for which the first-time vaccinees have higher absolute titers than the previous-year vaccinees are in the 5a.2 clade, which is the same clade as the 2021–2022 vaccine strain, A/Wisconsin/588/2019 (Fig. 8; Fig. S7). However, the modest differences in absolute neutralization titers between the two groups at 30 days post-vaccination are nearly gone by 182 days post-vaccination (Fig. 8). Overall, these results show that while the titer rise caused by vaccination is much larger in participants who have not been vaccinated the previous year compared to those who received a previous-year vaccine, the actual titers post-vaccination are fairly similar, with first-time vaccinees having slightly higher titers at 30 days, but not 182 days post-vaccination.

The fact that we generate full neutralization landscapes for each study participant against many viral strains enables the identification of specific viruses with reduced neutralization by specific participants. For example, several individuals have markedly lower post-vaccination neutralization titers against a small subclade of the 5a.1 clade (represented by viruses A/Belgium/H0017/2022 and A/England/220200318/2022) than against other viral strains. A subset of individuals with limited neutralization against this particular subclade is shown in Fig. 9. These more poorly neutralized viral strains carry HA mutations P137S and G155E (in H1 numbering), both of which are in known antigenic regions of the H1 HA (the Ca and Sa epitopes, respectively) (15). Only one other strain in the library has either of these mutations: P137S is found in the 5a.2 strain A/Norway/25089/2022, but that strain does not have reduced neutralization by any of these participants (Fig. 9). This fact suggests that the reduced neutralization of the subclade containing A/Belgium/H0017/2022 and A/England/220200318/2022 could be due to G155E, although it also could be due to the combined action of several mutations in a fashion that depends on viral genetic background.

Inspection of Fig. S5 and S6 also reveals other instances where specific participants have enhanced or reduced neutralization of specific viral strains. Overall, these observations emphasize how neutralization titers are both participant- and virus-specific.

Influenza is a segmented, negative-sense RNA virus. Proper packaging of the viral RNA (vRNA) segments into the virion relies on segment-specific RNA sequences called “packaging signals” which are present at the 3′ and 5′ ends of each segment and include both coding and noncoding regions (52, 53) (Fig. S1A). To insert a unique barcode into the HA genomic segment without substantially disrupting vRNA packaging, we designed a chimeric construct with a duplicated packaging signal at the 5′ end of the negative-sense HA vRNA (Fig. S1B) (13, 14, 52, 54–57), see plasmid maps in https://github.com/jbloomlab/flu_seqneut_DRIVE_2021-22_repeat_vax/tree/main/plasmids. As all non-HA genes used to generate the barcoded influenza viruses are from the laboratory-adapted A/WSN/1933 (H1N1) strain, we used terminal packaging signals from this strain to better facilitate the incorporation of HAs from various recent H1N1 strains into the virions (Fig. S1B) (55–57). We chose A/WSN/1933 (H1N1) for practical and biosafety reasons. Specifically, this laboratory-adapted strain produces high titer stocks by reverse genetics, and its extensive laboratory adaptation is thought to have reduced its risk to humans. At the 3′ end of the negative-sense HA vRNA, the native 3′ packaging signal [defined as the 32-nucleotide noncoding region and first 67 nucleotides of the coding sequence of HA (52)] from A/WSN/1933 HA was used (Fig. S1B). We then incorporated the ectodomain from the HA strain of interest, spanning amino acids 23–520. The transmembrane domain and cytoplasmic tail of the HA, spanning the last 46 amino acids, are replaced with the homologous 46 amino-acid region from the WSN strain followed by a pair of consecutive stop codons (144 nucleotides in total). We incorporated synonymous mutations into the last 105 nucleotides of this region encoding the C-terminal portion of the HA to remove the packaging signals without altering the encoded protein (Fig. S1B). This synonymous re-coding reduces homology between the duplicated region within the coding sequence and the terminal packaging sequence used for incorporating the segment into the virion, which may be important for reducing intra-segment recombination. Immediately following the consecutive stop codons, we incorporated a random 16-nucleotide barcode followed by the Illumina read 1 priming sequence. We then included the last 150 nucleotides of the negative-sense 5′ HA vRNA from the WSN strain (this is the full WSN HA 5′ packaging signal). A stop codon was introduced into what would have been the coding sequence of this final region to improve barcode retention: if recombination between the partially homologous regions immediately upstream and downstream of the barcodes occur, this stop codon should truncate the HA protein, therefore, yield a non-functional virus (Fig. S1B, and https://github.com/jbloomlab/flu_seqneut_DRIVE_2021-22_repeat_vax/tree/main/plasmids).

To select strains for inclusion in the barcoded influenza library, we first identified emerging clades using a Nextstrain (16, 17) build for human H1N1pdm sequences from 7 October 2022 (https://github.com/jbloomlab/flu_seqneut_DRIVE_2021-22_repeat_vax/tree/main/librarydesign). This was done using a method available through Augur on Nextstrain (similar to https://github.com/nextstrain/augur/blob/master/scripts/identify_emerging_clades.py). Clades were selected based on the following criteria: they must contain at least one HA1 amino acid mutation compared to the parent clade and have descendants that were represented at a global frequency of greater than 1%. This yielded a set of 33 clades; however, two clades were excluded. One of the clades was eliminated as all members of that clade overlapped with an included subclade, and the other was excluded as all sequences were from 2020 or before, indicating that this clade was likely no longer circulating. Next, we manually selected a sequence from each of these clades that contained the clade-defining mutations and were close to the node of the clade such that the selected strains did not include additional mutations which were not shared with all members of the clade. This yielded 31 strains that represented the diversity of viruses that were circulating in the fall of 2022. This list of strains is at https://github.com/jbloomlab/flu_seqneut_DRIVE_2021-22_repeat_vax/tree/main/librarydesign. In addition to the selected recent strains, the following vaccine strains were also included in the library: A/California/07/2009, A/Michigan/45/2015, A/Brisbane/02/2018, A/Hawaii/70/2019, and A/Wisconsin/588/2019.

The ectodomains of the hemagglutinin proteins from selected influenza virus strains were ordered as synthetic DNA fragments from Twist Biosciences and cloned with a barcoded fragment encoding the last 46 amino acids of WSN HA as three-segment assembly reaction into a construct containing the 3′ non-coding region of the packaging signal and signal peptide from A/WSN/1933 influenza HA and the Read 1 Illumina sequence and the full 5′ packaging signal from A/WSN/1933 virus using Hifi Assembly Mastermix (NEB). The backbone for this cloning reaction was a pHH21 plasmid (19) which had been modified to contain the packaging sequences of WSN HA, with restriction enzyme cut sites flanking the 3′ packaging signal, specifically, a BamHI site, which results in cleavage after nucleotide 67 of the viral HA genomic segment, and an XbaI immediately preceding the read 1 sequence to allow for incorporation of unique barcodes into each plasmid. The full plasmid map for the vector used for this reaction is made available as 2805_WSNHAflank_GFP_NterminalWSNHA_duppac-single-stop.gb in https://github.com/jbloomlab/flu_seqneut_DRIVE_2021-22_repeat_vax/tree/main/plasmids. Each strain was cloned with a minimum of two barcodes/strain, and these were sequence-verified with whole plasmid sequencing through Primordium Labs. See https://github.com/jbloomlab/flu_seqneut_DRIVE_2021-22_repeat_vax/tree/main/plasmids for the Genbank maps for all plasmids.

To generate actual viruses from these plasmids, the unidirectional reverse genetics plasmids encoding all barcoded variants of the same HA sequence were pooled at equal concentrations and used to generate influenza viruses containing unique barcodes by reverse genetics, with the WSN pHW18* series of plasmids for all seven non-HA viral genes (20), and pHAGE2-EF1aInt-TMPRSS2-IRES-mCherry-W (58). Briefly, 5e5 293T cells and 5e4 MDCK-SIAT1-TMPRSS2 cells (58) were plated in a 6-well dish in D10 media (DMEM supplemented with 10% heat-inactivated fetal bovine serum, 2 mM L-glutamine, 100 U per mL penicillin, and 100 µg per mL streptomycin). Approximately 24 hours after plating cells, a master mix containing 240 ng/well of each of the reverse genetics plasmids for WSN virus and 80 ng of the TMPRSS2 containing plasmid was prepared in 100 μL DMEM with 3 μL BioT reagent. This was incubated for 15 minutes at room temperature and then added dropwise to the plate. Approximately 20 hours post-transfection, media was removed, and cells were washed once with phosphate-buffered saline (PBS) and 2 mL Influenza Growth Media (Opti-MEM supplemented with 0.1% heat-inactivated FBS, 0.3% bovine serum albumin, 100 µg per mL of calcium chloride, 100 U per mL penicillin, and 100 µg per mL streptomycin) was added to cells. At ~65 hours post-transfection, cell supernatants containing the barcoded influenza viruses were collected and centrifuged for 4 min at 845 × g to remove cell debris. Aliquots of clarified viral supernatant were then frozen at −80°C for storage. To expand viruses to high titer, all virus pools generated by reverse genetics were passaged once in MDCK-SIAT1-TMPRSS2 cells (58). For this, MDCK-SIAT1-TMPRSS2 cells were seeded at a density of 4e5 cells/well in a six-well dish in D10, and at 4 hours after seeding were washed with PBS and 2 mL of Influenza Growth Medium was added to each well. Each well was then inoculated with 200 μL of the virus supernatant from reverse genetics and the virus was allowed to grow in the cells for ~40 hours. Cell supernatants were then collected and clarified by centrifugation at 845 × g for 4 minutes, and aliquots of clarified viral supernatants were stored at −80°C for use in the sequencing-based neutralization assay.

To generate a pooled library where each strain was represented at approximately equal transcriptional titers, we used a sequencing-based titering method. First, an equal volume mix of all of the passaged barcoded virus strains that were to be incorporated into the virus library was made and titered. The TCID50 for this initial pool was ~35,900 TCID50/uL on MDCK-SIAT1-TMPRSS2 cells. Serial dilutions of this pool were used to inoculate a single well of a 96-well plate containing 50,000 MDCK-SIAT1 cells, and after 16 hours the viral barcode frequencies were determined by deep sequencing as described in a later section. We determined the relative fraction of all barcode counts that corresponded to each strain (Fig. S3A, left panel). We then took the fraction of the library that each strain should represent in an ideal equal-representation pool (the reciprocal of the number of strains, in this case, 1/36) and divided this number by the fraction in the initial equal-volume pool to find the ratio at which that strain should be added. We multiplied this ratio by 200 μL to determine the volume to add for each strain to generate the final pool, with sufficient volume to run ~100 neutralization assay plates. The resulting pool was titered by TCID50 and found to be at ~39,600 TCID50/uL. We sequenced cells infected with the transcriptional-titer balanced pool and found that all strains were within approximately twofold representation of each other (Fig. S3A, right panel). We also examined the relative representation of each barcode for each strain (Fig. S3B). There is variability in barcode representation for each strain (Fig. S3B); however, as the replicate barcodes for each strain were rescued by reverse genetics as a pool, these ratios cannot be manipulated at this stage. Independent reverse-genetics rescues could be performed for each barcoded strain to be able to balance barcodes, rather than strains. That approach would substantially increase the number of rescues needed to prepare a barcoded library; however, by improving the balance of barcodes in the library, the total library concentration needed in the final assay to obtain reasonable coverage of each barcode would likely be lower. Therefore, better normalization at this step could facilitate larger library construction. The yield of the virus library obtained with the approach described above depends on the titers of the viruses included in the library. If viruses are generated at titers between 1,000 TCID50/uL–10,000 TCID50/uL, with 2 mL of each strain per strain included in the library, and a similar library size (36 strains), one should have sufficient volume to run between 26 and 260 sequencing-based neutralization assay plates, or sufficient volume to perform measurements against the library for ~200–2,000 serum samples. However, the relevant titer is the transcriptional titer of the library (how many virions can infect cells to transcribe the HA barcode), not TCID50/uL. This transcriptional titer should be measured experimentally for each new library as described in a later section. The necessity of this rebalancing depends on how different the titers are between the different strains one wishes to include in the library. For this library, only a small number of strains were poorly represented in the equal volume pool; however, if there is significant variability in growth of library strains, this rebalancing would be essential to ensure that each strain is adequately represented in the final pool.

Our laboratory has previously confirmed that barcodes introduced into the influenza HA gene segment using a broadly similar strategy to that in this paper are stably retained by the virus (13, 14). However, to confirm that the barcode region was genetically stable within the exact HA constructs in the current study, we analyzed the HA sequences of the virus library generated by pooling the passaged viral stocks. For this analysis, RNA was isolated from the viral supernatant of the pooled library using the RNeasy PLUS Mini Kit (Qiagen, 74104). Lysis buffer (RLT) was mixed with the viral supernatant sample and 70% ethanol was added, the sample was then applied to the RNAeasy purification column and RNA was isolated according to the protocol supplied by the manufacturer. The RNA for each sample was eluted in 22 μL of RNase-free water. Reverse transcription was performed using 5 μL of the extracted RNA according to the manufacturer’s protocol for the SuperScript III First-Strand Synthesis SuperMix kit (Thermo Scientific), with a gene-specific primer full-length_HA_Fwd: 5′-AGCAAAAGCAGGGGAAAATAAAAACAACC-3′. Following cDNA synthesis, PCR using full-length_HA_Fwd primer and full-length_HA_Rev: 5′- AGTAGAAACAAGGGTGTTTTTCCTTATATTTCTGAAATCC-3′ was performed to amplify the HA segment. A plasmid containing barcoded HA was used as a control to verify the expected size. We examined both a single virus strain after initial reverse genetics and blind passage and found that the full-length barcoded construct was maintained. In the pooled library, some deletions are observed in the HA segment (Fig. S8A). We Sanger sequenced all bands observed in the gel and verified that all of them included the barcode region, the largest band is the intact, full-length HA segment. Significant diversity in nucleotide identity was observed within the barcode region of this segment, as expected. The next largest band contains a deletion within the HA ectodomain. This segment also had diversity barcodes associated with it, indicating that this deletion likely occurred in multiple barcoded HAs. Finally, the smallest segment is a deletion product missing much of the HA protein sequence, but containing only a single barcode, indicating that it likely originated in one of the passages from a single strain. This barcode is not associated with any of the sequence-confirmed HAs included in the library, so no curves were generated for this sequence. The presence of some amount of internal deletions in influenza gene segments including HA is common (59) and does not pose any problem for our experimental approach as long as the barcode region is retained (which is the case here). To confirm that internal deletions were not rampant in any other influenza genomic segments (which would suggest high amounts of defective interfering particles), we performed cDNA synthesis and PCR as described above using primers that annealed to all eight segments of the influenza genome as previously described by others (60). All segments appear to be present in these stocks, suggesting that they are not dominated by defective interfering particles (Fig. S8B).

Our approach used an RNA spike-in control to convert sequencing counts to an absolute fraction of viral infectivity remaining (Fig. 1). To generate a spike-in control single-stranded RNA that contained known barcodes and resembled the vRNA of barcoded influenza viruses, we first cloned the protein sequence for green fluorescent protein in place of the HA coding sequence in the plasmid backbone we used for the barcoded HAs (see schematic in Fig. S1C and the full plasmid map 3808_WSNHAflank_GFP_WSNtermini_definedbarcodes_forSpikeinControl, available in https://github.com/jbloomlab/flu_seqneut_DRIVE_2021-22_repeat_vax/tree/main/plasmids), and then sequence-verified 10 plasmids containing unique barcodes with whole plasmid sequencing. We then generated a linear PCR template from the pooled set of barcoded plasmids using primers complementary to the termini of the vRNA encoded in the plasmid, full-length_HA_Fwd: 5′-AGCAAAAGCAGGGGAAAATAAAAACAACC-3′ and full-length_HA_Rev: 5′- AGTAGAAACAAGGGTGTTTTTCCTTATATTTCTGAAATCC-3′, respectively. The PCR product was run on a 1% agarose gel at 85 V for 35 minutes and the gel is extracted using the Nucleospin Gel Extraction Kit (Takara). We then used PCR with KOD Polymerase 2× Mastermix (Millipore Sigma) to append a T7 promoter sequence, using the following primers: Add_T7_Rev: 5′-TTACGATAATACGACTCACTATAGGGAGTAGAAACAAGGGTGTTTTTCCTTATATTTCTG-3′ and full-length_HA_Fwd: 5′-AGCAAAAGCAGGGGAAAATAAAAACAACC-3′, with 500 ng of the purified linear template. Then we used a T7 Ribomax Express (NEB) in vitro RNA synthesis kit to generate single-stranded negative-sense RNA that resembles the vRNA molecule. To remove DNA template after in vitro transcription, 0.5 μL of RQ DNAse was added to the reaction and incubated for another 15 minutes at 37°C. The Monarch RNA Cleanup Kit (NEB) was used to purify the resulting RNA. The product was quantified by Nanodrop, diluted to 200 pM, aliquoted for use in each experiment, and stored at −80°C. We verified that the RNA spike-in stock was free of DNA template by performing cDNA synthesis using the iScript Select cDNA Synthesis Kit, with different concentrations of template (0.12–12 ng) with and without the reverse-transcriptase enzyme. PCR amplification was run with the Round 1 PCR primers and the product was run on a 1% agarose gel at 85 V for 1 hour. No product was observed in the no-template and no-enzyme controls, indicating that spike-in was free of contaminating DNA.

The sequencing-based neutralization assays require that the number of sequencing counts for each HA viral barcode normalized by the spike-in RNA counts is directly proportional to the number of virions encoding that barcode that infected cells (Fig. 1). This condition will be met as long as the multiplicity of infection (MOI) is not too high since at low to modest MOIs an increased number of virions infecting cells will produce more barcoded viral RNA, but at very high MOIs the capacity of infected cells to produce viral RNA becomes saturated. To identify the amount of virus library that can be added to a well while maintaining linearity between the amount of input library and the resulting spike-in normalized HA counts, we performed a dilution series where wells were infected with different amounts of virus library and then sequenced as described above using 50,000 MDCK-SIAT1 cells/well with serial dilutions of the virus library at known TCID50s. As described above, we added a constant amount of the RNA spike-in control (50 μL of 2 pM RNA) in lysis buffer to each well containing cells infected with the virus pool. We then sequenced the barcodes in the well and identified the highest amount of virus library we could add and still be in the range where a twofold decrease in input virus library led to a twofold decrease in the counts of viral barcodes relative to spike-in control barcodes. We found that adding 25,000 TCID50/well (or a volume of ~0.6 μL/well of library) and 50,000 MDCK-SIAT1 cells per well was appropriate for this particular library; note however that the TCID50 to cell ratio could vary for other libraries depending on the ratio of TCID50 to transcriptionally active virions. Our library has a total of 110 barcoded variants (the 36 different HAs each have 2–4 barcodes), corresponding to ~220 TCID50 per well when using the amount of virus chosen for our experiments. These values present some anticipated limits to the library size that might be able to be characterized by the sequencing-based neutralization assays. At a coverage of at least 100 TCID50/well, and a maximum inoculum of 25,000 TCID50/well, with three barcodes per variant, we anticipate that up to 80 variants per well could be characterized with this method. A smaller number of barcodes per variant might be employed to increase the library size. The experiments described in this subsection should be repeated for each new library, as different stocks of viruses may have different numbers of transcriptionally active particles.

To verify that competition among different viruses via serum absorption of antibodies was not substantially affecting the measured neutralization titers, we performed the following experiment. First, we tested whether repeating the experiments using the virus library at half concentration altered the NT50s. As shown in Fig. S9 (top panel), the measured NT50s were not dependent on the virus concentration used. In addition, we prepared two additional pools of virus library, one lacking a 5a.1 clade virus (A/South Africa/R16462/2021) and one lacking a 5a.2 clade virus (A/India-Pune-Nivcov2221170/2022). We then tested whether the titers measured against all other viruses in these one-virus dropout pools were similar to those measured for the full library pool. As shown in Fig. S9 (middle and bottom panel), the NT50s measured with the one-virus dropout pools were similar to those measured for the full virus pool.

All viral infections for barcode sequencing were performed using MDCK-SIAT1 cells. We used these cells rather than the MDCK-SIAT1-TMPRSS2 cells to limit secondary viral replication (TMPRSS2 cleaves HA in producing cells to activate it for infection) (61–63). The cells were plated with influenza virus library (at indicated concentrations) with 50,000 cells per well in 96-well plates in IGM in a total volume of 150 μL. At 16 hours post-infection, cells were washed once with 150 μL PBS, and then lysed by a 5-minute incubation with 50 μL iScript RT-qPCR Sample Preparation Reagent (BioRad) (64), containing the RNA spike-in control added at a final concentration of 2 pM,. We use 1 μL of this lysate in a 10 μL cDNA synthesis reaction, using the iScript cDNA Select Synthesis kit according to the manufacturer’s instructions for synthesis with a gene-specific primer that binds immediately upstream of the barcode region (cDNA_Fwd 5′ CTCCCTACAATGTCGGATTTGTATTTAATAG-3′). The resulting cDNA product was then amplified and the Read 2 sequence was added with PCR, using a forward primer that binds immediately upstream of the barcode and appends the Read 2 sequence (Rnd1_Fwd 5′-GTGACTGGAGTTCAGACGTGTGCTCTTCCGATCTCTCCCTACAATGTCGGATTTGTATTTAATAG-3′) and a reverse primer that binds to the terminal sequence of the vRNA (Rnd1_Rev 5′-AGTAGAAACAAGGGTGTTTTTCCTTATATTTCTGAAATCC-3′), in a 50 μL reaction, with 5 μL of the cDNA reaction as template, according to the manufacturer’s instructions, using KOD Polymerase Hot Start 2× Mastermix (Sigma). This amplification was performed with 20 cycles total, using an annealing temperature of 56°C, and elongation time of 20 seconds, ending with a final elongation step at 70°C for 2 minutes. A second round of PCR with unique dual indexing primers was performed to uniquely label each sample, using primers designed as follows UDI_i7 5′- CAAGCAGAAGACGGCATACGAGATnnnnnnnnnnGTGACTGGAGTTCAGACGTGTGCTCTTCCGATCT-3 and UDI_i5 5′-AATGATACGGCGACCACCGAGATCTACACnnnnnnnnnnACACTCTTTCCCTACACGACGCTCTTCCGATCT-3′, using a 20-cycle protocol with KOD Polymerase Hot Start 2× Mastermix (Millipore Sigma), 20 seconds at 95°C, 10 seconds at 66°C, and 20 seconds at 70°C, and a final 2-minute elongation at 70°C. After round 2 indexing PCR, samples were pooled at equal volume, and run on a 1% agarose gel at 85V for 40 minutes. The band was extracted and purified with a Nucleospin Gel Extraction Kit (Takara). The gel-purified DNA was then further purified with Ampure XP beads (Beckman Coulter) by adding 2× volume of DNA binding beads to the sample, incubating 5–10 minutes, then applying to a magnet and incubating again for 10 minutes, removing the supernatant, washing the pellet 2× with 80% ethanol, drying for approximately 10 minutes, taking care not to overdry the pellet, and then resuspending in elution buffer (5 mM Tris/HCl, pH 8.5).

The indexed DNA was quantified with the Qubit 1X dsDNA High Sensitivity Kit (Thermo Scientific), diluted to 4 nM, and submitted for Illumina Sequencing on the NextSeq as either P1 or P2 depending on whether one or two plates were sequenced, to give an average number of reads per well of approximately 500,000–1,000,000 reads per well. This read depth is larger than necessary, but it allows for variability in sample loading in the absence of more rigorous, concentration-dependent pooling, which may introduce additional errors given the number of samples pooled per run. We use a 50 bp read length, and custom UDI indexes based on the sequences from Twist Biosciences. We identified an issue in the sequencing of one of the spike-in control barcodes that began with “GG.” Low read counts were detected for this barcode in a consistent well on each plate, though the affected well was different between plates that used either index set A or index set B, indicating that the index sequence may play a role in whether this occurs. The affected barcode was excluded from the analysis. To eliminate this issue, barcodes that start with “GG” should be avoided in library design, and indices starting with “GG” should be also avoided (as recommended by Illumina).

A step-by-step protocol is available on protocols.io (https://doi.org/10.17504/protocols.io.kqdg3xdmpg25/v1). Briefly, the setup for a sequencing-based neutralization assay resembles other 96-well plate neutralization assays (Fig. 1). Prior to performing the assays, serum was treated with receptor-destroying enzyme (RDE) and heat inactivated. RDE treatment was done to remove sialic acids in the serum which may interfere with the assay. This protocol was performed as previously described (39). To set up sequencing-based neutralization assays, virus and RDE-treated serum aliquots were thawed at room temperature. Then serial dilutions of serum were made in Influenza Growth Medium in a 96-well plate for a final total volume of 50 μL at 2 x of the intended serum concentration to be incubated with the virus. For this study, we performed threefold dilutions of serum with a starting concentration of 1:20 to capture a large range of potential NT50 values. Pooled virus library was diluted to an appropriate concentration, as determined by transcriptional titering, in Influenza Growth Medium (for this library, 70 μL of the stock library in 5,430 μL of IGM), and 50 μL of diluted virus library was added to the plate (Fig. 1). Virus-serum mixtures were incubated at 37°C with 5% CO2 for 1 hour. Following incubation of virus and serum, 50 μL of MDCK-SIAT1 cells at 1e6/mL in Influenza Growth Media was added to the plate, and the plate was then incubated at 37°C with 5% CO2 (Day 1, Fig. 1). At approximately 16 hours post-infection (Day 2, Fig. 1), supernatant was removed, infected cells were washed once with 150 μL PBS to remove any residual virus. Cells were lysed with 50 μL of iScript RT-qPCR Sample Preparation Reagent (BioRad) per well containing the RNA spike-in control at 2 pM final concentration (Day 2, Fig. 1), cDNA for the barcoded region of the vRNA was synthesized and PCR was performed to prepare samples for Illumina sequencing as described in detail above.

To illustrate the potential this new sequencing-based neutralization assay, we selected human serum samples from a randomized vaccine trial of repeated annual vaccination of Flublok done in Hong Kong in 2020–2025 (DRIVE-1 study) for testing. The DRIVE-1 study (ClinicalTrials.gov: NCT04576377) is a randomized placebo-controlled trial of repeated annual influenza vaccination in adults 18–45 years of age. We enrolled participants from the community if they met eligibility criteria and did not meet exclusion criteria. Individuals were eligible to participate if they were between 18 and 45 years of age, capable of providing informed consent and intending to reside in Hong Kong for at least the next 2 years. Potential participants were excluded if they had received influenza vaccination in the preceding 24 months if they were included in a priority group for influenza vaccination such as being a healthcare worker, if they had a diagnosed immunosuppressive condition or were taking immunosuppressive medication, or if they had severe allergies or bleeding conditions that contraindicated intramuscular influenza vaccination. We enrolled 447 participants between 23 October 2020 and 11 March 2021, and randomly allocated them to five separate groups, with annual receipt of vaccination (V) or saline placebo (P) as follows: V-V-V-V-V, P-V-V-V-V, P-P-V-V-V, P-P-P-V-V, or P-P-P-P-V. For the present study, we selected samples provided by 15 participants from years 1 and 2 of the V-V-(V-V-V) and P-V-(V-V-V) groups, that is, participants randomly allocated to receive vaccination in both years or placebo in year 1 (2020–2021) and vaccination in year 2 (2021–2022). The vaccine used was Flublok, a quadrivalent recombinant-HA vaccine produced in insect cells, including 45 mcg HA of each of the four included strains as recommended by the World Health Organization for the northern hemisphere in the 2020–2021 and 2021–2022 (65, 66). The relevant H1N1 vaccine strains in 2020–2021 and 2021–2022 seasons were A/Hawaii/70/2019 (H1N1) (clade 5a.1) and A/Wisconsin/588/2019 (H1N1) (clade 5a.2). We collected 9 mL sera immediately before administration of vaccine/placebo, and again at scheduled follow-up visits after approximately 30 days and 182 days each year. The study protocol was approved by the Institutional Review Boards of the University of Hong Kong (ref: UW19-551) and of the University of Chicago Biological Sciences Division (ref: IRB20-0217). Written informed consent was obtained from all participants. For the present analysis, we selected day 0, 30, and 182 sera collected from study year 2 in a subset of participants who had a pre-year 1 vaccination HAI titer <10 against the A/Hawaii/70/2019 (H1N1) vaccine strain, and who were assigned to either receive the placebo or the Flublok Vaccine (A/Hawaii/70/2019) in study year 1 (2020–2021 season) but all received the Flublok Vaccine (A/Wisconsin/588/2019) in year 2 (2021–22 season), that is, two groups of participants who were vaccinated in 2021–2022 and who had either received prior year vaccination in 2020–2021 (“2xVax”) or not (“1xVax”).

To analyze the data from the sequencing-based neutralization assay, we developed the modular seqneut-pipeline (https://github.com/jbloomlab/seqneut-pipeline), which we have made publicly available so that it can be easily used by others. This pipeline, in turn, makes use of the neutcurve Python package (https://github.com/jbloomlab/neutcurve) that was previously developed by some of us to automate the process of fitting neutralization curves. All analyses performed for this paper, which make use of seqneut-pipeline and neutcurve are available https://github.com/jbloomlab/flu_seqneut_DRIVE_2021-22_repeat_vax. We have also made available at https://jbloomlab.github.io/flu_seqneut_DRIVE_2021-22_repeat_vax/ interactive HTML versions of key plots and analysis notebooks to enable easier exploration of the data.

The initial processing of sequencing-based neutralization assay data is done by plate, as controls are included on each plate. First, we process Illumina barcode sequencing using the parser https://jbloomlab.github.io/dms_variants/dms_variants.illuminabarcodeparser.html to determine the counts of each barcoded strain in the well of the 96-well plate in each selection condition. We then determine the average counts per barcode in each well and filter wells on the plate with insufficient barcode counts. In addition, we examine the fraction of reads from the RNA spike-in barcodes and filter any wells with too few of these barcodes. Next, check that the fraction of counts from each barcode is consistent across all no-serum control wells and exclude any barcodes with highly variable representation in multiple control wells. See the config.yml in the main project GitHub repo (https://github.com/jbloomlab/flu_seqneut_DRIVE_2021-22_repeat_vax) for the exact values used in these quality control steps.

We then calculate a fraction infectivity (F_b,s) for each viral barcode (b) in each sample (s) as follows: F_b,s = [c_b,s/Σ_n c_n,s]/[median_s0(c_b,s0/Σ_n c_n,s0)] where c_b,s is the count of barcode b in sample s, Σ_n c_n,s is the sum of the counts for all neutralization standard barcodes (n) in sample s, c_b,s0 is the counts of viral barcode b in the no-serum sample (s0), and median_s0(c_b,s0/Σ_n c_n,s0) is the median taken across all no-serum control samples for the counts of the viral barcode divided by the total counts of all neutralization standard barcodes.

There is some variability in the amount of vRNA each influenza-infected cell produces (13, 59). In general, this variability is averaged out by the fact that each barcoded virus infects hundreds of cells in our assay setup; nonetheless, rarely, a single barcode may be over-represented in a particular well of the plate. We can detect this issue by examining fraction infectivity for each barcode and excluding measurements where a given barcode in a single well yields an extremely high fraction infectivity (such as a value >5). After filtering is completed, we clip all fraction infectivity values at an upper limit of one and then fit a Hill curve for each barcoded virus (Fig. 1), dropping any curves with a very poor coefficient of determination (r²). Again, see the config.yml for the exact filtering criteria used in this study.

Next, we calculate a neutralization titer for each barcoded strain. An NT50 can be calculated either by determining the midpoint of the fraction infectivity curve or directly calculating the value at which the fraction infectivity would be 50%. NT50s obtained by these two methods are generally very well correlated for each viral barcode; however, for curves with an upper plateau that is less than one (which sometimes happens if there is a slightly lower representation of the RNA spike-in control barcode in the serum dilution wells than the no-serum control wells), then midpoint measurement may be a better representation of neutralizing titer. For this reason, we used the midpoint to calculate NT50s for this study.

Finally, we report the NT50 for each viral strain by as the median of the NT50s measured for each barcode corresponding to that strain. In instances in which there is a specific viral barcode with an NT50 that is an outlier with respect to the median, we manually examined the curves for these barcodes and excluded them from the NT50 calculations if they were aberrant for an identifiable reason. At least two barcodes per strain must pass all filters for an NT50 to be reported in this study. See https://jbloomlab.github.io/flu_seqneut_DRIVE_2021-22_repeat_vax/ to visually inspect all neutralization curves generated for each plate and serum in this study.

A traditional one-virus versus one-serum neutralization assay was used for the validations in Fig. 4. These assays were performed as previously described (24, 25, 67) (see https://github.com/jbloomlab/flu_PB1flank-GFP_neut_assay) with the modifications described here. Briefly, viruses for neutralization assays were generated by reverse genetics as described above, using a single barcoded HA per strain. In place of the plasmid encoding the PB1 segment of A/WSN/1933, a plasmid (pHH-PB1flank-eGFP) which encodes GFP within the flanking sequences from the PB1 segment to allow for packaging into virions was used (24). These viruses were produced and passaged in cells constitutively expressing the PB1 protein (68). For reverse genetics, 5e5 293T-CMV-PB1 and 5e4 MDCK-SIAT1-CMV-PB1-TMPRSS2 cells were seeded in D10 media in a six-well plate (67). Approximately 20 hours after seeding, cells were transfected with 250 ng of each reverse genetics plasmid using BioT transfection reagent according to the manufacturer’s instructions. At 20 hours post-transfection, D10 media was removed, and cells were washed once with PBS and 2 mL NAM (Neutralization Assay Media, consisting of Medium-199 supplemented with 0.01% heat-inactivated FBS, 0.3% BSA, 100 U/mL penicillin, 100 µg/mL streptomycin, 100 µg/mL calcium chloride, and 25 mM HEPES) was added to the well. At approximately 48 hours post-transfection, the transfection supernatants were harvested, clarified by centrifugation at 845 × g for 4 min, aliquoted, and frozen at −80°C. Viruses were passaged at low MOI (0.02) in MDCK-SIAT1-CMV-PB1-TMPRSS2 cells, as previously described. We titered the passaged virus stocks by TCID50 in 1e5 MDCK-SIAT1-CMV-PB1-TMPRSS2 (67). In addition, a test was performed to determine the amount of virus to add to each well such that the fluorescence signal decreases twofold with a twofold dilution of virus stock. Once the appropriate concentration of virus was determined, RDE-treated serum was serially diluted in a 96-well plate in NAM and, the virus was added. The viruses and serum were incubated at 37°C for 1 hour to allow antibodies in the serum to bind to GFP-encoding influenza viruses containing the HA segment of interest and then 5e4 MDCK SIAT1-CMV-PB1 cells were added to each well. At approximately 20 hours post-infection, fluorescence was measured using a Tecan M1000 plate reader. Fraction infectivity was calculated as compared to a no-serum well. Neutralization curves were then fit to the data using the neutcurve package (https://jbloomlab.github.io/neutcurve).

The HAI titers reported in Fig. 4C were measured using the protocol described in reference (69). Briefly, serum samples were treated with receptor-destroying enzyme for 18 hours and heat-inactivated at 56°C for 30 minutes. Serum samples were serially diluted twofold with a starting dilution of 1:10 and incubated for 30 minutes with four agglutinating doses (4AD) of test antigens. Titers were read as the reciprocal serum dilutions that inhibited agglutination of 0.5% turkey red blood cells.