Introduction

Seasonal influenza is estimated to cause up to half a million deaths every year and is associated with a high economic and clinical burden^1–5. Influenza vaccines reduce the disease caused by influenza virus infection although their effectiveness has varied from 10 to 60%^6,7. Seasonal influenza vaccines are composed of representative strains of H1N1, H3N2 and one or two influenza B lineages. The recommended vaccine strains are constantly evaluated and, if needed, updated twice a year, aiming to ensure antigenic match with the predominant circulating strains. The antigenic evolution of human influenza viruses is typically determined in hemagglutination inhibition tests, which captures the antigenic drift of the viral hemagglutinin (HA)^8,9. The neuraminidase (NA) is the second major influenza A and B virus surface antigen. Like HA, NA is subject to antigenic drift, indicating that immune responses against NA exert a selective pressure on influenza A and B viruses¹⁰. Furthermore, antibodies against NA correlate with protection against influenza^11–16. These findings have contributed to the renewed interest in the development of influenza vaccines that by including NA antigen would have improved efficacy and breadth of protection ^17,18.

We previously reported the antigenic evolution of N2 NAs derived from human H3N2 strains that were isolated between 2009 and 2017¹⁹. That work expanded the knowledge of N2 antigenic diversity and identified at least four N2 antigenic groups. Furthermore, we demonstrated that the antigenic assessment of recombinant N2 NA in mice provides similar results to those obtained in ferrets, which is considered the gold standard small animal model to study immunity to and pathogenesis of influenza viruses.

Since its introduction in the human population, the 2009 pandemic H1N1 (H1N1pdm) virus has been subjected to immune selection pressure that has resulted in antigenic drift of HA and NA²⁰. Our aim was to obtain a high-resolution insight in the antigenic relatedness of NAs derived from H1N1pdm viruses that were isolated in the period of 2009-2020. We selected N1 NAs derived from 17 human H1N1 strains that were isolated between 2009 and 2020. The corresponding recombinant NAs were produced and used to immunize mice. NA inhibition (NAI) activity was then determined against a panel of H1N1, H6N1, and an H5N1 virus. NAI data analysis revealed that H1N1 viruses isolated between 2009 and 2020 could be classified into 2 antigenically distinct groups. Furthermore, a machine learning model was generated that allowed the estimation of the antigenic distances between the N1 NAs. Application of this model identified K432E and I321V as the two major classifiers of N1 antigenicity in the tested strains.

Results

Immunogenicity of recombinant N1 NAs

Based on amino acid sequence diversity, and using minimax linkage clustering²¹ implemented in R protoclust package, 17 representative N1 NAs were selected from the GISAID database²² from the 2009-2020 period. The phylogenetic relatedness, based on the amino acid sequences of these selected N1 NAs, together with the NA of 5 other randomly selected H1N1 isolates per year from 2009 to 2020 is shown in Figure 1a. The selected N1s are distributed along all different clades indicating that minimax linkage clustering selected a diverse panel of N1 representatives.

Figure 1.

Next, we generated recombinant N1 NA immunogens based on the 17 selected H1N1pdm viruses. These immunogens consisted of a tetramerizing tetrabrachion zipper fused at the N-terminus of the NA catalytic head and were expressed as secreted proteins by transfected CHO cells, as previously described.¹⁹ The specific activity of the purified N1 NA proteins was determined using MUNANA as a substrate. All 17 selected N1 NAs were active, except the NA derived from A/Ankara/TR-40/2011 (Figure 1b). All 17 N1 NAs could also bind oseltamivir, indicating that the active site in the recombinant NAs was correctly folded (Figure 1c). To assess variability on immunogenicity and antigenicity, two distinct batches (named 946 and 162) of N1 NA derived from A/Stockholm/25/2012 were produced and characterized, with comparable results (Figure 1b and 1c).

Next, groups of 5 mice were immunized with 1µg of N1 NA, adjuvanted with AF03, in a prime boost regimen with 2 weeks interval to generate immune sera (Figure 1d). Two weeks after the last immunization mice were terminally bled, and sera were pooled per group. Analysis of the sera by ELISA using the respective recombinant N1 NA as coated immunogens showed that all groups seroconverted (Figure 1d).

NA inhibition-based antigenic analysis reveals two distinct N1 NA antigenic groups

The NAI activity of the N1 NA mouse immune sera was determined using an enzyme-linked lectin assay (ELLA) against 15 different HxN1 viruses: H6N1 reassortants with N1 derived from A/California/07/2009, A/Michigan/45/2015, A/Hubei Wujiagang/SWL310/2013, A/Cameroon/9766/2017, A/Oman/5532/2017, or A/Wisconsin/588/2019; H1N1 viruses A/Brisbane/02/2018, A/Hawaii/70/2019, A/Wisconsin/588/2019, A/Lviv/N6/2009, A/Stockholm/25/2012, A/Israel/Q-504/2015, A/Hungary/12/2016, and A/Ireland/84630/2018); and NIBRG-14 carrying HA and NA derived from A/Vietnam/1194/2004 (H5N1) and the remaining genes of A/Puerto Rico/8/34. The NAI titers were normalized to the respective ELISA titers and plotted as a heatmap, with columns indicating the influenza virus strains used as a source of NA in ELLA and rows representing the sera raised against the respective recombinant N1 NAs (Figure 2). The serum samples are ordered according to the phylogenetic relatedness of the recombinant N1 NA that was used as an immunogen. Hierarchical clustering of the tested HxN1 strains revealed 3 major N1 NA antigenic groups: a first group of strains from 2009-2013 (AG1), a second group formed by strains from 2015-2020 (AG2), and the H5N1 strain NIBRG-14 which stands apart and was poorly inhibited by the mouse immune sera (Figure 2). The phylogenetic relatedness based on amino acid sequences reveals two major groups of the selected N1 NA immunogens, one composed of strains that circulated from 2009 until 2012, defined here as phylogenetic group 1 (PG1) and a second group comprising N1 NAs from more recent strains, 2015 to 2020, defined here as phylogenetic group 2 (PG2). The immune sera raised against N1 NAs belonging to PG1 displayed stronger inhibition against AG1 and tended to inhibit NA from AG2 viruses less well. Conversely, sera raised against N1 NAs classified in PG2 displayed reduced inhibition of AG1 viruses and stronger inhibition of strains from AG2.

Figure 2.

Although the average NAI of NIBRG-14 (H5N1) by sera generated against PG1 N1 NAs is higher than for anti-PG2 sera, this difference did not reach statistical significance (Figure 2 – figure supplement 1a). Comparison of the normalized NAI IC₅₀ titers of the sera raised against the 2 batches of N1 A/Stockholm/25/2012 (162 and 946) revealed a significantly lower NAI of batch 946 against all but one of the HxN1 viruses (Figure 2 – figure supplement 1b). Despite this difference in immunogenicity the antigenic profile of the sera raised against the two A/Stockholm/25/2012 N1NA batches is highly similar (R =0.86) (Figure 2 – figure supplement 1c).

Visualization of the antigenic differences using multidimensional scaling as proposed by Smith et al, reveals a spatial distribution that is in line with the hierarchical clustering (Figure 2 – figure supplement 2).

A Machine Learning tool to predict N1 NA antigenic distances and identify relevant amino acid substitutions

Next, we calculated the sera-to-sera average antigenic differences and used these differences to generate a random forest model (RF) as detailed in the Methods section (Figure 3a). The model was optimal when using a combination of at least 35 randomly sampled features in 600 decision trees. The N1 antigenic distances predicted by the RF model strongly correlated with the observed values (Figure 3b). Interestingly, the scoring of the individual amino acid substitutions enabled the identification of those substitutions that had the highest impact on the NAI-based antigenic distance. The mean squared error increase (MSEI) reflects the effect of feature permutation, where a larger error increase suggests a greater relevance of the feature (amino acid substitution) to antigenic distance. Another measure of importance is the usage of a given feature as the root of the decision tree: the greater the number of trees that use a given feature at the root, the higher its importance. The resulting RF model predicted that substitutions of N1 NA amino acid residues at positions 431 and 321 stand out on both scores (Figure 3c). Representation of the MSEI values on the surface of N1 NA shows that K432E is in proximity of the catalytic pocket, while I321V is located at the lateral side of N1 NA (Figure 3d).

Figure 3.

Discussion

The protective effect of anti-NA antibody responses against influenza has been known for many decades and repeatedly observed in clinical studies^{14–16,24–26}. Nonetheless, currently used virus-derived influenza vaccines induce a poor response against NA²⁷. Therefore, the addition of NA to influenza vaccines to improve their effectiveness has been proposed²⁸. However, among other reasons, the labile nature of the NA tetramer and limited understanding of NA antigenic drift hamper major efforts to bring this concept to clinical reality²⁹. Previously, using a panel of mouse and ferret sera and recombinant, soluble tetrameric NA immunogens we determined the antigenic landscape of influenza N2 NA of human H3N2 viruses isolated between 2009 and 2017. Here we performed a similar analysis of human H1N1 NAs using immune sera from mice only. We selected a genetically divergent panel of H1N1 pandemic N1 NAs derived from human H1N1 viruses that circulated from 2009 until 2020. These N1 NAs were produced as soluble recombinant proteins, stabilized as a tetramer by a heterologous tetrabrachion domain. This strategy allows the induction of an immune response which is free of anti-HA antibodies, and thus enabled us to capture the NA inhibition without potential interference by such antibodies. The characterization of the recombinant N1 NAs included the determination of enzymatic activity and binding to oseltamivir. The A/Ankara/TR-40/2011 was the only inactive protein in the MUNANA assay, although it showed binding to oseltamivir and induced NAI antibodies, which indicates that, despite non-detectable activity, the protein was correctly folded.

The inclusion of 2 independent groups of mice that were immunized with different preparations (batches 946 and 162) of NA from A/Stockholm/25/2012 revealed a degree of variability in immunogenicity. Batch 946 had lower enzymatic activity, oseltamivir binding ratio and induced slightly lower homologous ELISA titers compared to batch 162. Even after normalization of the NAI IC₅₀ values with the ELISA titers, NAIs were on average 0.66 fold lower in the immune sera raised against batch 946 compared to batch 162 of A/Stockholm/25/2012 N1 NA. This result points out that minor antigenic differences observed in different antigens may be masked by experimental bias, or differences in antigen quality.

NAI determination using the panel of 18 anti-N1 NA sera against 15 HxN1 viruses revealed that pandemic N1 NAs can be classified in at least 2 distinct groups. The antigenic group 1 comprises H1N1 viruses that circulated in the early years after the start of the 2009 H1N1 pandemic and in our panel is represented by strains from 2009 until 2013. The antigenic group 2 is formed by more recent strains, isolated between 2015 and 2020. The antigenic classification agrees well with the phylogenetic relatedness of the stains which is characterized by the substitutions V13I, I34V, L40I, N44S, N200S, V264I, N270K, I314M, I321V, N386K, and K432E (Figure 2).

The calculation of antigenic distances based on sera-to-sera average differences was used to build a RF model that can predict the antigenic distance between two N1 NAs based on their amino acid sequences. This RF model also allows to identify which amino acid substitutions are most likely responsible for the antigenic differences. The substitutions V13I, I34V, L40I, and N44S are located in the N1 stalk, and are absent in the recombinant tetrabrachion stabilized NAs. Therefore, these substitutions were not included in the RF analysis. Among the remaining substitutions, based on the RF model, I321V and K432E stand out as the most important ones (Figure 3c). As those substitutions coincide in our panel, our analysis is not able to specify which of these two substitutions contributes the most, or rather if both are needed to differentiate between the two antigenic N1 NA groups. Interestingly, I321V and K432E are the only substitutions that are never reversed among the two antigenic groups. Using a panel of seven monoclonal antibodies isolated from individuals that had been infected with H1N1pdm, Yasuhara et al., reported the antigenic distinction between the prototype A/California/04/09 strain (clade 1) and more recent isolates such as A/Yokohama/94/2015 (clade 6B.1). The substitution at position N390K at the lateral face together with the substitution K432E, located at the edge of the catalytic site, impacted binding of all seven isolated monoclonal antibodies³⁰. Our results are in line with Yasuhara et al., pointing to a role of K432E, combined with substitution at the lateral surface. In our panel, the individual impact of the substitution K432E or the rather conservative substitution I321V, which is also located on the lateral surface, remains to be addressed.

The substitutions that were reversed in at least one of the strains were ranked lower in importance by RF scores. For example, substitution N200S is present in recent N1s (AG2) but also in A/Jericho/18/2012, a N1 that belongs to AG1. The substitutions V264I and K270N are reversed in A/Hungary/12/2016 and thus seem not to impact the antigenicity of NA from that strain. The substitution I314M is also reversed in A/Israel/Q-504/2015 and, again, is not associated with a major impact on N1 antigenicity. The N386K substitution is present in all strains belonging to PG2 and this substitution results in the loss of an N-glycosylation site. In our panel the N386S mutation occurs in A/Stockholm/25/2012, which also results in loss of the N-glycosylation site, and did not result in higher NAI titers against recent N1s, suggesting that the loss of an N-glycosylation at position N386 does not explain the major antigenic difference between early G1 and the recent G2 H1N1 viruses. Gao et al. compared the antigenicity of NA from vaccine recommended strains and identified the N386K substitution as an important contributor to the antigenic change ³¹. However, this conclusion was based on NAI titers obtained from sera raised against recent strains (GD19, BR18 and MI15) that were tested against CA09 mutants. The presence of an N-glycan at position 386, reduces recognition of CA09 N1 by sera raised against more recent N1 NAs and not the other way around, which may explain the discrepancy with our results. All vaccine recommended strains since 2015 bear the full set of substitutions that characterize the antigenic group 2 (V13I, I34V, L40I, N44S, N200S, V264I, N270K, I314M, I321V, N386K and K432E).

In conclusion, by analyzing the NAI using a diverse panel of anti-H1N1 NA immune sera and HxN1 NA sources, our study identified two distinct antigenic groups of H1N1 NAs from strains that have been isolated between 2009 and 2020. The substitutions I321K and K432E are the critical markers of this antigenic discrimination, although their individual contribution remains to be determined. Despite its limitations, such as the inability of RF to predict the impact of substitutions that are not present in the training data, our study offers a framework and primary dataset that may guide future selection of NA candidates for a next generation of influenza vaccines.

Materials and methods

Design and production of recombinant proteins

The CD5 secretion signal, an amino-terminal His-tag and the thrombin cleavage site followed by tetrabrachion domain, were fused to the N1NA head domain (aa 74-469) coding sequence and cloned under the transcriptional control of the CMV promoter in the pCDNA3.4 plasmid. The recombinant tetNAs were expressed in a mammalian cell culture system as previously described with modifications¹⁹.

Virus generation and amplification

The HxN1 read out panel viruses are 6:2 reassortants on an A/Puerto Rico/8/1934 backbone (PB2, PB1, PA, NP, M, NS). The NA gene from the desired strain was matched with an HA gene matching the H1N1 strain or derived from H6 A/Mallard/Sweden/81/2002. The eight influenza segments were cloned in the bidirectional plasmid pHW2000S, and the resulting plasmids were transfected into HEK293FT cells using a lipid-based transfection reagent³². HEK293FT cells were co-cultured with MDCK ATL cells during transfection, or MDCK ATL cells were added 24hrs post transfection of HEK293FT cells to allow for efficient amplification of Influenza virus. The virus harvest was done 4 to 6 days post transfection, (this initial recovery stage referred to as P0). Harvested P0 virus was further passed on MDCK ATL cells to generate P1 and P2 working virus stocks that were characterized for infectious titer and NA sequence confirmation. Working virus stocks (P2) were passaged on specific pathogen free (SPF) white leghorn chicken eggs to further amplify virus stocks. Harvested allantoic fluid was clarified by centrifugation and characterized for infectious titer, NA gene sequence identity and NA enzymatic activity.

NA enzymatic activity

The enzymatic activity of recombinant tetNAs was determined with the fluorogenic small substrate 4-methylumbelliferyl-N-acetylneuraminic acid (MUNANA) as follows. Ten microliters of 1CmM of 4-MUNANA (Sigma cat # M8639) in 200CmM sodium acetate buffer (pH 6.5) containing 2CmM CaCl₂ and 1% butanol was incubated with 40Cµl of PBS solution containing 30, 15, 7.5, or 3.75Cng of tetNA. Conversion to 4-methylumbelliferone (4-MU) was monitored every 2Cmin for 1Ch using a BMG Fluostar OPTIMA reader (excitation at 365Cnm and emission determined at 450Cnm). A standard curve with 800 to 25Cpmoles/well of 4-methylumbelliferyl (Sigma cat # M1508) was used to extrapolate the molar conversion of MUNANA per microgram of purified recombinant tetNAs.

Tamiflu-NA binding assay on Octet Bio-Layer Interferometry

Binding of NA to Oseltamivir Phosphate-Biotin (Tamiflu-Biotin, GlycoTech Inc.) was measured on the Octet-HTX (ForteBio, Sartorius). All measurements were performed in Sartorius Kinetic Buffer (1xKB buffer containing PBS, 0.02% Tween20, 0.1% BSA, 0.05% sodium azide). The Tamiflu-Biotin was immobilized on streptavidin biosensors (ForteBio, Sartorius) at a concentration of 0.25-0.5 µg/ml in 1xKB buffer. Binding of NA to Tamiflu was initiated by dipping Tamiflu bound biosensors into sample wells containing a 2-fold dilution series of recombinant NA (0.8–200 nM in 1xKB buffer) at 25 °C with a 5 min association step and a 5 min disassociation step. The level of binding response is proportional to the concentration of recombinant NA. The response at the end of the association step was used for determination of kinetics, affinity, and quantitation regarding the NA bound to the biosensor tip.

Mouse experiments

The mouse immunization experiments were conducted according to the Belgian legislation (Belgian Law 14/08/1986 and Belgium Royal Decree 06/04/2010) and European legislation on protection of animals used for scientific purposes (EU directives 2010/63/EU and 86/609/EEC). Experimental protocols were all approved by the Ethics Committee of the Vlaams Instituut voor Biotechnologie (VIB), Ghent University, Faculty of Science (EC2022-056). Female BALB/c mice, aged 6–7 weeks, were purchased from Janvier. The mice were housed in a specified pathogen-free animal house with food and water ad libitum. Animals were immunized twice, with an interval of 14 days, intramuscularly in the right quadriceps. Immunizations containing 1 µg of recombinant protein were performed with a volume of 50 µl (25 μl of antigen in PBS +C25 μl of AF03 per dose) ³³.

ELISA

Anti-tetNA IgG levels in mouse sera were determined by capture ELISA using tetNAs in Pierce nickel-coated plates (cat. # 15442). Recombinant proteins were diluted to 0.5 µg/ml in DPBS (Life Technologies cat. # 14040-182) Then, 50 µl of the coating antigen solution was added to each well and the plate was incubated at room temperature for 1 h on a shaking platform. The plates were then washed three times with PBS-T (Sigma cat. # P3563-10PAK) and blocked for 1 h with 1% BSA. After blocking, the plates were washed once with PBS-T and incubated with a threefold serial dilution, starting at a 1/100 dilution, of serum in DPBS, 0.5% BSA, 0.05% Tween-20 for 2 h at room temperature on a shaking platform. The plates were then washed five times with PBS-T and incubated with a 1:5000 dilution of anti-mouse IgG-HRP (GE Healthcare cat. # NA931-1ml) in DPBS, 0.5% BSA, 0.05% Tween-20. The 3,3’,5,5’-tetramethylbenzidine (TMB) substrate (BD cat. # 555214) was added after three washes with PBS-T and the reaction was stopped after 5 min by addition of 50 µl of 1 M H₂SO₄. The optical density (OD) in each well was determined at 450 nm and, as a reference, 655 nm using an iMark Microplate Absorbance Reader (Bio-Rad). The end point titer was determined for each serum sample by scoring the dilution that resulted in an OD that was equal to or two times higher than the background OD obtained from pre-immune control sera dilution series.

ELLA to determine NAI titers

Fetuin (Sigma cat. # F3385) was diluted into coating buffer (KPL cat. # 50-84-01) to a concentration of 25 µg/ml and 50 µl was added to the wells of Nunc MaxiSorp plates (Thermo Fisher cat. # 44-2404-21), which were incubated overnight at 4°C. The coated plates were then washed three times with PBS-T (Sigma cat. # P3563-10PAK) and incubated overnight with 60 μl of a dilution of H6N1 or H1N1 and 60 µl of twofold serial dilution of serum, starting at a 1/20 dilution, in sample buffer (1× MES VWR cat. # AAJ61979-AP: 20 mM CaCl₂, 1% BSA, 0.5% Tween-20). These dilutions of the H6N1s or H1N1s viruses correspond to the 70% maximum activity of NA from the respective viruses as determined in the ELLA. Fetuin-coated plates were then washed three times with PBS-T and incubated for 1 h with a solution of PNA-HRP (cat. # L6135-1MG, Sigma) at 5 μg/ml in conjugate diluent (MES pH 6.5, 20 mM CaCl₂, 1% BSA). The plates were washed three times with PBS-T, TMB substrate was added, and then the plates were incubated for 5 min before the reaction was stopped by the addition of 50 µl of 1 M H₂SO₄. The optical density was measured at 450 nm and, as a reference, 655 nm in an iMark Microplate Absorbance Reader (Bio-Rad). Half maximum inhibitory concentrations (IC₅₀) values were determined by non-linear regression analysis (GraphPad Prism software).

Random Forest algorithm

The RF algorithm was adapted from Yao et al., 2017 and described previously^19,34. Briefly, the non-conserved positions in the multiple sequence alignment of N1 head domain amino acid sequences of all tested N1s were selected. The feature matrix (X) was created by assessing matches/mismatches on each non-conserved position for each pair of N1 head domain amino acid sequences as:

where X_i represents the non-conserved positions in the multiple sequence alignment, m is the number of non-conserved positions (38 in total), V_i represents the N1 head domain sequence for each tested NA, n is the number of tested sera (18 in total), 0 represents a match, 1 represents a mismatch in the alignment.

The response vector (Y) is composed of the antigenic distances for each pair of tested sera. These antigenic distances were determined based on NAI, similar to the method used to determine HI-based immunogenic distances that was described by Cai et al. The distances between sera were calculated as —, which represent average distances between two serum samples among all tested sources of NA³⁵. We applied the RF function from the randomForest package in R to construct the model. After optimization, we set the bootstrapping (tree) number to 600 and the number of features to 35 for each tree. Accuracy was determined as the RMSE, calculated using the RMSE function from the yardstick package in R. Additionally, predicted antigenic distances of the out-of-bag samples were compared to the observed distances by fitting a linear model and calculating the R2 value. To calculate the importance of each feature (Xi), the metrics were calculated using the measure importance function of the randomForestExplainer package in R.

Antigenic cartography

Antigenic cartography was performed as described by Smith et al., 2004, using the R CRAN package (https://CRAN.R-project.org/package=Racmacs)²³. The antigenic map was optimized 5000 times and the best optimization was kept (lowest stress value). Uncertainty was measured by bootstrapping 1000 repeats with 100 optimizations per repeat.

Additional information

Disclosure of funding sources

This work was funded by Sanofi Vaccines.

Author Contributions

JPPC, TY, AS, LA, MGM, JH, BB, KR, GC, JZ, LH, EC, SR, PD, DT, YC performed experimental work. JPPC, TUV, and XS designed the experiments, analyzed and interpreted the data. JPPC and XS wrote the manuscript and are accountable for accuracy and integrity of the work. All authors proofread and critically reviewed the manuscript.