Epistatic interactions between neuraminidase mutations facilitated the emergence of the oseltamivir-resistant H1N1 influenza viruses.

Oseltamivir-resistant H1N1 influenza viruses carrying the H275Y neuraminidase mutation predominated worldwide during the 2007-2009 seasons. Although several neuraminidase substitutions were found to be necessary to counteract the adverse effects of H275Y, the order and impact of evolutionary events involved remain elusive. Here we reconstruct H1N1 neuraminidase phylogeny during 1999-2009, estimate the timing and order of crucial amino acid changes and evaluate their impact on the biological outcome of the H275Y mutation. Of the 12 neuraminidase substitutions that occurred during 1999-2009, 5 (chronologically, V234M, R222Q, K329E, D344N, H275Y and D354G) are necessary for maintaining full neuraminidase function in the presence of the H275Y mutation by altering protein accumulation or enzyme affinity/activity. The sequential emergence and cumulative effects of these mutations clearly illustrate a role for epistasis in shaping the emergence and subsequent evolution of a drug-resistant virus population, which can be useful in understanding emergence of novel viral phenotypes of influenza.

The evolutionary course of influenza A viruses is shaped by interplays among mutation, reassortment, and natural selection[1]. Influenza A viruses, like all RNA viruses, have a high mutation rate[2]; whether a mutation can spread at a population level (epidemiologic fitness) is dependent upon its impact on viral biologic fitness (replication fitness within a host and transmission fitness between hosts)[3]. Understanding the impact of mutations and interaction of mutations on viral fitness is, therefore, critical for a mechanistic understanding of viral phenotype emergence.

The influenza virus neuraminidase (NA) inhibitors oseltamivir and zanamivir are the options currently approved in the U.S. for immediate control of influenza virus infection. Their clinical use, however, provides a selection force to drive emergence of resistance within treated individuals. Before 2007, resistant viruses were detected only infrequently during NA inhibitor treatment[4-7] and very rarely during surveillance[7-9], suggesting that those drug-driven resistant viruses had little epidemiologic fitness. However, during the 2007–2009 influenza seasons, oseltamivir-resistant H1N1 viruses surged from <1% to >90% prevalence worldwide[10-12]. Such spread of resistance at population level was not attributed to oseltamivir use in individuals, but to global transmission of the resistant viruses carrying the NA H275Y mutation[13, 14], suggesting these H275Y-mutant viruses had acquired advantageous epidemiologic fitness. A mechanistic understanding of such drug-independent resistance spread would give us insights to the adaptability and evolution of drug-resistant influenza viruses.

Recent studies have advanced our understanding of the biological properties of the H275Y-mutant viruses related to their different epidemiologic fitness outcomes. Genetically, the NA genes of most H275Y-mutant viruses were closely associated with the genetic 2B clade (represented by A/Brisbane/59/2007 [BR07]) of H1N1 viruses but not with the other three clades (clade 1, represented by A/New Caledonia/1999 [NC99]; clade 2A, represented by A/Solomon Island/23/2006 [SI06]; and clade 2C)[15-19]. This clade-specific resistance distribution suggested a link between biologic fitness and genetic context of the H275Y-mutants. Indeed, phenotypically, NC99-like H275Y-mutants manifested greater biological cost relative to their respective wild-type counterparts than did BR07-like mutants, as measured by growth in cells, mice, and ferrets[20-23] and by their NA affinity[15, 16, 19], and cell surface accumulation[24].

Several mutations have been identified elsewhere in the NA that can counteract the adverse effects of the H275Y mutation. It has been found that the D344N[16, 25], R222Q, and V234M[24] NA substitutions can counteract the reduced NA affinity and surface accumulation caused by the H275Y mutation; therefore, these mutations are “permissive” for the H275Y mutation. Another study confirmed that changing the permissive substitutions to the non-permissive substitutions (Q222R, M234V) compromised the replication fitness of a clade 2B H275Y-mutant virus in vitro and in ferrets[26]. While illuminating, the identification of these NA permissive mutations has not provided a full understanding of the evolutionary path and molecular process involved in the fitness changes of the H275Y-mutant viruses.

Here we reconstruct the molecular evolutionary path of the NA protein of seasonal H1N1 viruses from the NC99 to BR07 genetic lineage, during 1999–2009. We then evaluated the biological outcomes of the H275Y mutation in different NA genetic contexts at different stages of the path. We further investigate the chronological order and nature of NA mutations and their impact on the phenotypic outcome of the H275Y mutation in vitro and in vivo. These biologic fitness assessments include NA functional activity (accumulation and substrate affinity), replication in cells (plaque morphology), and viral growth and transmissibility in ferrets. We find that multiple mutations, having either permissive or compensatory epistatic interactions with the H275Y, sequentially alter the outcomes of the H275Y mutation in NA functionality and virus biological fitness. Our findings show biologic evidence for epistasis in which the phenotypic outcome of a mutation depends on the absence or presence of others. These results support a role for epistasis in shaping the evolution of influenza viruses, and in the emergence and spread of viral phenotype.

Results

NA protein sequence evolution from 1999 to 2009

The NA sequences of seasonal H1N1 viruses isolated during 1999–2009 showed a major genetic divergence, from the NC99 (clade 1) to BR07 (clade 2B) lineage[15, 16, 19]. We first sought to identify changes in the NA protein sequence during this lineage transition that would explain the different H275Y mutation outcomes in the NC99 and BR07 lineages. Comparison of the consensus NA gene sequences revealed ten amino acid substitutions (H45N, K78E, E214G, R222Q, V234M, G249K, T287I, K329E, D344N, G354D, excluding H275Y) differentiating the lineages (Supplementary Table 1) and one substitution (D354G) differentiating wild-type and H275Y-mutant NA sequences within the BR07-like lineage. None of the substitutions resided at the NA active site (Supplementary Fig. 1).

We reconstructed the N1 NA protein’s evolutionary path, focusing on the period 1999–2009, and estimated the time and order of occurrence of the amino acid substitutions by using a tip-calibrated relaxed molecular clock phylogenetic tree reconstruction method[27] (Fig. 1). The NC99-like NA population (the tree trunk in Fig. 1) acquired most of the mutations in succession during 1999–2007, rather than through a single selective sweep, evolving into the BR07-like NA population in early 2007; multiple mutations were seemingly acquired at a few time-points. Several intermediate NA populations (clustered away from the trunk Fig. 1) that had acquired some of the mutations co-circulated for some time but were extinguished before 2007. The H275Y mutation occurred only sporadically in those populations, in which it ended as an NA evolutionary terminus (the tip of a branch Fig. 1). In mid-2007 (2007.485, Bayesian credible intervals, 2007.302–2007.698), soon after the BR07-like NA population was established, the H275Y was fixed into the population. Only one mutation (D354G) was acquired after the H275Y, thereby completing the NA protein evolution such that oseltamivir resistance dominated the population during 2008–2009.

Figure 1

Temporally-structured maximum clade credibility phylogenetic tree of the NA gene of seasonal H1N1 viruses, showing the timeline of amino acid substitutions

Ancestral state reconstruction of tested amino acids is annotated on the tree backbone to indicate when specific mutations were fixed in the population. The trunk color change indicates a >95% posterior probability of amino acid fixation. Purple bars on tree nodes indicate the Bayesian credible intervals of ancestral divergence time estimates (>50% posterior probability). Isolates with H275Y mutation prior to H275Y fixation into the population are indicated with red branches.

Different outcome of the H275Y mutation in NA functionality during the NC99-to-BR07 transition

The NA protein sequence changes seen in our phylogenetic calculations prompted us to evaluate the possible accompanying phenotypic changes in NA functionality. To this end, we expressed the NA proteins (using equal quantities of NA-containing plasmids to partially mimic NA expression by equal quantities of infectious virus particles) and assessed their total activity and accumulation. To quantify the different NA proteins, we adapted a method for uniformly epitope-tagging recombinant proteins[24]. After comparing untagged, C-terminal–tagged, and N-terminal–tagged recombinant NABR/59 proteins, we chose the C-terminal epitope tag for use in all NA proteins, as it did not significantly affect NA accumulation or activity (Supplementary Fig. 2). We compared lineage-representative NA proteins as well as intermediate proteins that had partial complements of the 10 substitutions (Supplementary Table 2).

We first evaluated NA enzyme activity (Fig. 2a). Across the two lineages, the BR07-like wild-type NABR/59 had enzyme activity ~3–4 times that of two NC99-like wild-type proteins (NAGA/17 and NAMEM/13), while the intermediate NAMS/03 (it had acquired the R222Q and D344N substitutions but not V234M yet and therefore it had atypical intermediate genotype), had activity near that of NABR/59. Within the NC99-like lineage, the three pairs of wild-type and H275Y-mutant NA proteins differed only at position 275; the H275Y-mutants (NAGA/20 (H275Y) and NAMEM/13-H275Y) had activity of only ~50% that of their wild-type counterparts, while the activity of the intermediate mutant NAMS/03 (H275Y) was reduced least (~20%). Within the BR07-like lineage, the activity of H275Y-mutant proteins did not significantly differ from that of wild-type NABR/59, although the two representative mutants (NANJ/15 (H275Y) and NANY/1692 (H275Y)) had slightly higher activity and two intermediates (NABR/59-H275Y and NANY/3467 (H275Y)) that had not acquired the D354G or D344N mutations showed slightly lower activity.

Figure 2

NA protein phenotypes of seasonal H1N1 influenza viruses of the NC99 and BR07 lineages

(a) Neuraminidase enzyme activity of total NA proteins expressed by equal quantities (0.5 μg) of the respective plasmids. Rep: representative NA protein of the indicated lineage; Int: intermediate NA protein of the indicated lineage, fully or partially conforming to the 10 consensus lineage residues. The residue at the NA protein 275 position (wild-type [wt] or H275Y) is indicated in parentheses. NA-H275Y is shown when the H275Y was introduced by mutagenesis. Data were normalized to NABR07 protein. Dotted and dashed lines indicate the activity level of lineage-representative NAGA/17 and NABR/59 proteins, respectively. (b) Representative western blot showing total accumulation of NA protein variants expressed by equal quantities (0.5 μg) of the respective plasmids. Black and grey indicate wt and H275Y-mutant NA proteins, respectively. Panels show detection by anti-HA-tag (upper), anti-β actin (middle), and anti-GFP (lower). (c) Representative flow cytometric plots of cell surface accumulation of NA protein variants expressed by equal quantities (0.5 μg) of respective plasmids. The pairs differed only at position 275, with the exception of the last panel, the BR07-like NAs. (d) Km values of the NA proteins. The dotted and dashed lines indicate the value of lineage-representative NAGA/17 and NABR/59, respectively. All graphs show mean ± s.e.m. from three to five times independent experiments. * p<0.05, two-tailed t test, versus NABR/59; † p<0.05, two-tailed t test, versus its counterpart wt NA at left side.

We next examined whether the observed difference in the NA activities reflected differences in their protein amount. Across the two lineages, two BR07-like wild-type proteins (NABR/59 and NACA/27) showed higher accumulation than three NC99-like wild-type proteins (NAGA/17, NAMEM/13 and NAMS/03) (Fig. 2b, Supplementary Fig. 3a). Within the NC99-like lineage, accumulation of the three H275Y-mutant proteins was substantially lower than that of their wild-type counterparts, but accumulation of the intermediate NAMS/03-H275Y was reduced least. Within the BR07-like lineage, accumulation of the five H275Y-mutant NA proteins was similar to that of the wild-type NABR/59. Cell surface accumulation of the NA proteins showed the same trends of change as observed in total NA accumulation (Fig. 2c).

As the protein amount differed among the NAs, we next examined each NA-specific activity by standardizing each NA’s activity relative to its protein level, which reflected each NA’s intrinsic conversion ability (Supplementary Fig. 3b). The differences in the standardized NA-specific activities showed similar patterns to those in the total NA activity. Across the two lineages, the standardized activity of wild-type NABR/59 was ~2–3 times of that of the two representative NC99-like wild-type NAs (NAGA/17, NAMEM/13). Within each lineage, all the H275Y NAs showed slightly lower standardized activity than their wild-type counterparts, except that two representative BR07-like H275Y NAs (NANJ/15(H275Y) and NANY/1692 (H275Y)) showed higher standardized activity than the NABR/59. The distinct levels of NA-specific activity between the two lineages reflected that the NAs had varying intrinsic conversion ability in addition to their different expression level.

Enzyme Km is an intrinsic enzyme property and its value is inversely correlated with the enzyme substrate-binding affinity and enzyme conversion ability: the higher Km value indicates lower binding affinity and reaction rate, leading to lower activity. We next examined whether the different levels of standardized NA-specific activity reflect different Km values. Across the two lineages, the Km values of two BR07-like wild-type proteins (NABR/59 and NACA/27) were only ~25%–30% those of two NC99-like wild-type proteins; the intermediate NAMS/03 protein had an intermediate Km (Fig. 2d). Within lineages, the H275Y-mutant proteins had consistently higher (~2 fold) Km values than their wild-type counterparts; only the intermediate NANY/3467 (H275Y), which had not acquired the D344N, had an unusually high Km within the BR07 lineage. The lineage-associated Km levels corroborated previous reports that BR07-like viruses had greater NA affinity than the NC99-like viruses[15, 19, 28].

Taken together, these results revealed two NA phenotypic transformations during the NC99-to-BR07 NA evolution. First, the accumulation and substrate affinity of the wild-type BR07-like NA proteins increased greatly, resulting in inherently higher functionality. Second, the outcome of the H275Y mutation in net NA function differed across lineages: the H275Y-generated NA defects were severe in the NC99-like NA proteins but minimal or absent in the BR07-like NA proteins. These results indicate that certain mutations acquired by the BR07-like NA proteins contributed to the observed NA phenotypic transformations.

Different outcome of the H275Y mutation in virus transmissibility during the NC99-to-BR07 transition

The observed NA phenotypic transformations during NC99-to-BR07 NA evolution prompted us to assess the transmissibility of representative viruses isolated at different stages of the evolution path, an essential factor for virus spread. To this end, we chose a ferret model to recapitulate virus human transmission. Although so far it is recognized as the best model to recapitulate human transmission of influenza virus, ferret model has limitations such as outbred model and random factors in sneezing activity etc, so we chose to assess multiple wild-type and H275Y-mutant viruses of each lineages, looking for common observations for each lineage. We assessed the growth and transmissibility of selected viruses, via direct contact and respiratory droplets in influenza-naïve ferrets. Transmission was indicated by nasal virus shedding and seroconversion (Fig. 3, Supplementary Table 4).

Figure 3

Transmissibility of seasonal H1N1 influenza viruses of the NC99 and BR07 lineages in ferrets

(a) Virus titers in nasal wash samples from individual donor ferrets (black), direct-contact (DC) ferrets (blue), and respiratory-droplet (RD) ferrets (red). Arrows indicate the day (day 1 p.i for donors) when recipient ferrets were housed with donor ferrets. (b) Virus titers in DC ferret nasal wash samples at day 1 post-contact. Each dot represents an individual ferret. * p<0.05, two-tailed t test, versus MEM/13. (c) Virus titers in RD ferret nasal wash samples on the first day of detection (the day post-contact), which is indicated above each column.

In the NC99 lineage, both the wild-type (MEM/13/06) and H275Y-mutant (GA/20/06 and MS/03/02) viruses caused productive infection in 3 of 3 inoculated donor ferrets, but the two H275Y mutants were shed at much lower titers on day 1 post-inoculation (p.i.), (~3 and 1.5 log10TCID50 lower, respectively) (Fig. 3a). At day 1 post-contact (p.c.), only the wild-type virus had been transmitted to 3 of 3 recipients by direct contact, and neither H275Y-mutant had been transmitted to all recipients, and their titers were ~3 log10TCID50 lower (Fig. 3b). The wild-type virus was transmitted via respiratory droplets to 2 of 3 ferrets at days 3 and 5 p.c., while the both H275Y mutants were transmitted to only 1 of 3 ferrets, at days 7 and 9 p.c. (Fig. 3c). Therefore, the transmission of the two NC99-like H275Y-mutants was commonly slower and less efficient than the NC99-like wild-type virus, suggesting that the transmissibility of the NC99-like H275Y-mutants was impaired, consistent with their severe NA defects observed.

In the BR07 lineage, the wild-type (BR/59/07) and three H275Y-mutant viruses (HW/28/07, NY/3467/09, and NY/1692/09) caused productive infection in all inoculated donors (3 of 3). Virus shedding between wild-type and H275Y-mutant virus only differed in the intermediate NY/3467/09 virus (~1.5 log10TCID50 lower in the mutant virus at day 1 p.i) (Fig. 3a). Each virus had been transmitted to 3 of 3 recipients by direct contact at day 1 p.c. (Fig. 3b), showing similar transmission kinetics, although the intermediate HW/28/07 virus was transmitted later to one recipient. Each virus was transmitted to 2 of 3 recipients via respiratory droplets, with similar kinetics (Fig. 3c), except that the representative H275Y-mutant NY/1692/09 virus was transmitted earlier (day 3 p.c.) than the wild-type BR/59/07 (day 7 p.c.) to one recipient. Therefore, the transmission efficiency of the three BR07-like H275Y-mutants was commonly similar to the BR07-like wild-type virus, showing that the transmissibility of BR07-like H275Y mutant viruses is minimally impaired or not impaired, consistent with their minimal or absent NA defects observed.

We next generated two H275Y-mutant reverse-genetics (rg) viruses, rgBR/59/07-NAH275Y and rgBR/59/07-NAH275Y+D354G, to represent the evolutionary acquisition of the H275Y and D354G mutations within the BR07 lineage. Both H275Y-mutants productively infected 3 of 3 inoculated ferrets (Fig. 3a). They were transmitted to 3 of 3 recipients via direct contact and to 2 of 3 via respiratory droplets, but with different kinetics. At day 1 p.c., only the NAH275Y virus had not been transmitted to all direct-contact ferrets (were transmitted only 2 of 3), showing slightly slower transmission than that of the wild-type and NAH275Y+D354G viruses (Fig. 3b); the NAH275Y+D354G virus was transmitted to 2 of 3 recipients via respiratory droplets at days 3 and 5 p.c., showing slightly faster transmission than did the wild-type and NAH275Y viruses (transmission at days 3 and 7 p.c.) (Fig. 3c). These results were consistent with the observations from the three BR07-like H275Y-mutant viruses and further confirmed that the H275Y mutation minimally reduced the transmissibility of the BR/59/07 virus and suggested that addition of D354G restored the transmissibility of the H275Y-mutant virus to the wild-type level. Another report[29] found that D354G not only increased the transmissibility of a resistant virus in guinea pigs but also to a level higher than that of the wild-type virus, suggesting that the extent of D354G increasing resistant virus transmissibility may depend on other viral protein context.

In summary, the growth and transmissibility of NC99-like and BR07-like wild-type viruses was indistinguishable despite their different NA phenotypes. However, the growth and transmissibility outcome of the H275Y-mutant variants differed greatly between the two lineages, as did the extent of H275Y-generated NA defects, suggesting that the NA defects was directly correlated with virus transmissibility.

Interaction between mutations alters the phenotypic outcome of the H275Y mutation in the BR07-like NA

We next sought to identify which of the 11 amino acid substitutions that occurred during the NC99-to-BR07 NA evolution account for the observed NA phenotypic transformations. To this end, we examined the impact of each individual mutation on NA functionality with or without the H275Y mutation.

We firstly adopted a strategy using plaque size to evaluate NA-mediated virus release during the replication cycle. A plaque is formed from a single virus particle via multiple cycles of replication, and release of the progeny viruses at each cycle is mediated by NA functionality. Thus NA mutation that could alter NA activity sufficiently would directly affect the speed of virus release and consequently determine the plaque size. We first validated the feasibility of this strategy by confirming that plaque size paralleled NA enzyme activity in two pairs of reverse-genetics viruses, a BR07-like pair (NABR/59 and NABR/59-H275Y) and an NC99-like pair (NAGA/17 and NAGA/20 (H275Y)) (each pair differs only at position 275) (Supplementary Fig. 4a, b). We next generated rgBR/59/07 viruses whose NA proteins carried the respective substitutions N45H+E78K, G214E, Q222R, M234V, K249G, I287T, E329K, N344D, and D354G (representing a change from BR07-like to NC99-like), in the presence and absence of H275Y. Four substitutions (Q222R, M234V, E329K, and N344D) in NABR/59 significantly reduced plaque size, with and/or without the H275Y (Fig. 4a, Supplementary Fig. 5a). These four positions and the final substitution site (354) were selected for further evaluation; the remaining five substitutions were unlikely to reduce virus growth appreciably and were not further evaluated.

Figure 4

Plaque size, NA activity, NA accumulation, and substrate affinity of the BR/59/07 NA proteins with different amino acid substitutions

(a) Plaque size of viruses carrying ten different NA amino acid substitutions, in the absence (white) or presence (gray) of the H275Y mutation. The ten amino acid residues were replaced with the corresponding NC99-like consensus residues. (N45H+E78K were substituted together, as both reside in the stalk region.) Data represent mean ± s.e.m. diameter of ~10 randomly selected plaques. Dashed line indicates mean plaque diameter of wild-type NABR/59. (b) Enzyme activity of total NABR07 and NABR07-H275Y variant proteins with 5 different amino acid substitutions, expressed by equal quantities (0.5 μg) of the respective plasmids. Dotted and dashed lines indicate mean activity of NABR/59and NABR/59-H275Y, respectively. (c) Representative western blot showing total accumulation of protein variants with substitutions; the variants were expressed by equal quantities (0.5 μg) of the respective plasmids. Black and grey indicate wt and H275Y-mutant NA, respectively. Panels show detection by anti-HA-tag (upper), anti-β actin (center), and anti-GFP (lower). (d) NA Km values of rgBR/59/07 viruses carrying the respective NA amino acid substitutions in the absence and presence of the H275Y mutation. Dotted and dashed lines indicate the values of NABR/59and NABR/59-H275Y, respectively. All graphs represent mean ± s.e.m. of three to five times independent experiments. * p<0.05, two-tailed t test, versus NABR/59; † p<0.05, two-tailed t test, versus NABR/59-H275Y.

The four substitutions (Q222R, M234V, E329K, and N344D) variably reduced enzyme activity in NABR/59 or NABR/59-H275Y (Fig. 4b), consistent with the reduced plaque size observed. Conversely, while not notably increasing plaque size, the D354G substitution significantly increased enzyme activity in both NABR/59 and NABR/59-H275Y; it increased the enzyme activity NABR/59-H275Y to a level equivalent or slightly higher than that of NABR/59. Two substitutions (Q222R, M234V) in both NABR/59 and NABR/59-H275Y greatly reduced the NA protein accumulation (Fig. 4c, Supplementary Fig. 5b) and the level of accumulation resulted by each mutation with H275Y was much lower than that by each mutation alone (such as H275Y-Q222R vs Q222R alone). Three substitutions (Q222R, E329K, N344D) in both NABR/59 and NABR/59-H275Y significantly increased the Km value ~2 folds compared with Km of NABR/59 and NABR/59-H275Y respectively (Fig. 4d), and the Km values resulted by each mutation with H275Y was much higher than that by each mutation alone (such as H275Y-Q222R vs Q222R alone). The three Km-changing mutations differently reduced NA-specific activity with Q222R causing biggest reduction (Supplementary Fig. 5c). The D354G substitution had no appreciable effect on protein accumulation or Km values but significantly increased the NA-specific activity (Supplementary Fig. 5c). The 354 position resides at the NA surface but faces toward the virion’s membrane (Supplementary Fig. 1) and therefore is structurally unlikely to affect substrate interaction directly. Thus D354G increased enzyme conversion ability by unidentified mechanisms other than Km level. In contrast, the three affinity-changing mutations (Q222R, E329K, and N344D) are located at the NA surface and face away from the virion’s membrane and tetramer interface (Supplementary Fig. 1), structurally allowing substrate interaction.

Together, our results identified two classes of mutation interaction in NABR/59 that sequentially altered the outcome of the H275Y mutation in net NA function. First, acquisition of R222Q, V234M, K329E, and D344N mutations by NABR/59 each increased NA activity and thus minimized subsequent H275Y-generated defects. Second, acquisition of the D354G mutation by NABR/59–H275Y increased NA activity to wild-type protein levels, resulting in unaltered net NA functionality. Thus, the higher levels of enzyme activity (Fig. 2a, Supplementary Fig. 3b) of NANJ/15 (H275Y) and NANY/1692 (H275Y) reflected their further acquisition of the D354G after the H275Y addition and demonstrated the role of D354G in increasing H275Y-NA activity.

Reconstruction of NC99-to-BR07 NA phenotypic evolution in vitro

We next tested whether the molecular determinants identified above were sufficient to explain the NA phenotypic transformations that accompanied the NC99-to-BR07 transition. The V234M, R222Q, K329E and D344N mutations (representing a change from NC99-like to BR07- like) were introduced into the NAGA/17 and NAGA/20 (H275Y) (already possessing the D354G mutation) NC99-like proteins, which differed only at position 275. All substitutions, except D344N, significantly increased the functional activity of NAGA/17 and, to a lesser extent, of NAGA/20 (H275Y) (Supplementary Fig. 6a). Enzyme activity paralleled the plaque size of the mutant viruses (Supplementary Fig. 6b). In NAGA/17, two substitutions (R222Q and V234M)greatly increased NA protein accumulation (Fig. 5a), and three (R222Q, K329E and D344N) variably reduced the Km value (Fig. 5b); all of these changes occurred to a lesser extent in NAGA/20 (H275Y). These findings confirmed that these four substitutions play a role in increasing overall NA function and in reducing the H275Y-generated defects, although no single substitution could alone fully transform the NA phenotype and fully offset the H275Y-generated defects. Thus, the intermediate level of NA accumulation observed in NAMS/03, NAMS/03 (H275Y) (which had acquired R222Q and D344N), and unusually high Km value in NANY/3467 (H275Y) (which had not acquire the D344N) (Fig. 2) can be explained by partial acquisition of these phenotypic determinants.

Figure 5

Total protein accumulation and Km values of the NC99-like NA proteins with different amino acid substitutions

(a) Representative western blot showing total accumulation of the NC99-like NA proteins (NAGA/17and NAGA/20 (H275Y)) with the indicated single substitutions; the variant NAs were expressed by equal quantities (0.5 μg) of the respective plasmids. Black and grey indicate wt and H275Y-mutant NA proteins, respectively. (b) NA Km values of rgBR/59/07 viruses carrying NAGA/17 or NAGA/20 (H275Y) with the indicated single substitutions. Dotted and dashed lines indicate Km values of NAGA/17 and NAGA/20 (H275Y), respectively. (c) Representative western blot showing total accumulation of NAGA/17 and NAGA/20 (H275Y) proteins with sequentially added substitutions (sub), the order is the same in next panel; NAs were expressed by equal quantities (0.5 μg) of respective plasmids. (d) NA Km values of rgBR/59/07 viruses carrying the NAGA/17 and NAGA/20 (H275Y) with sequentially added substitutions. Dotted and dashed lines indicate Km values of NAGA/17 and NAGA/20 (H275Y), respectively. All graphs show mean ± s.e.m. of three to five times independent experiments. *p<0.05, two-tailed t test, versus NAGA/17; † p<0.05, two-tailed t test, versus NAGA/20 (H275Y).

We then evaluated whether sequential addition of all four phenotypic determinants were sufficient to fully reconstruct the BR07-like NA phenotype. The four substitutions were consecutively added to the NC99-like NAGA/17 and NAGA/20 (H275Y) proteins in the order indicated by the phylogenetic reconstruction of ancestral states (V234M-R222Q-K329E-D344N or V234M-R222Q-D344N-K329E) (Fig. 1). The sequential addition of the four substitutions progressively increased the functional activity of NAGA/17 and NAGA/20 (H275Y) to a final level equivalent to that of NABR/59 and NABR/59-H275Y, respectively (Supplementary Fig. 7a). The increased NA activity closely paralleled increased plaque size (Supplementary Fig. 7b). Addition of V234M and R222Q, but not the two subsequent mutations, to NAGA/17 and NAGA/20 (H275Y) increased NA protein accumulation to a level equivalent to that of NABR/59 and NABR/59-H275Y, respectively; (Fig. 5c). Addition of V234M did not appreciably reduce Km, but subsequent addition of R222Q, K329E, and D344N progressively reduced the Km values of NAGA/17 and NAGA/20 (H275Y) (Fig. 5d) to that of NABR/59 and NABR/59-H275Y, respectively.

In summary, sequential addition of the four identified molecular determinants into the two NC99-like NA proteins (already possessing the D354G mutation) fully reconstructed the BR07-like NA phenotype. The impact of the interactions of the V234M and R222Q with the H275Y on total NA accumulation corroborated and expanded previous findings about surface NA accumulation [24]. The two subsequent substitutions, K329E and D344N, increased NA affinity and further reduced H275Y-generated NA functional losses.

The outcome of the H275Y mutation in other N1-subtype NA proteins

As the NA H275Y mutation is associated with oseltamivir resistance mainly in the N1 subtype, we evaluated the functional outcomes of the H275Y mutation in other N1 proteins, including those of two human 2009 pandemic H1N1 (pdmH1N1) (NACA/04 and NADM/524) (Fig. 6a, b) and two human H5N1 (NAVN/1203 and NAHK/213) viruses (Fig. 6c, d). The H275Y mutation caused loss of neuraminidase activity by 45–70% in these four NA proteins and greatly diminished their surface accumulation (Fig. 6). The extent of H275Y-generated defects observed in these four N1 proteins showed that that the H275Y mutation remained deleterious for NA functionality in these NA genetic contexts, suggesting that additional mutations would be required to restore NA functionality, in order for these four viruses to carry H275Y mutation to spread to any degree. This results demonstrates evaluation of the biology cost by the H275Y mutation at protein level can an informative analytic tool for assessing other emerging resistant viruses in other future studies.

Figure 6

Defects in enzyme activity and protein accumulation caused by the H275Y mutation in other N1-subtype NA proteins

Each panel shows a pair of NA protein variants differing only in the absence or presence of the H275Y mutation, including two pairs of pdmH1N1 NA proteins (a, b) and two pairs of hpH5N1 NA proteins (c, d). Enzyme activity of total accumulated NA protein is shown at left in each panel; Representative flow cytometry plots showing surface accumulation are at right. NA proteins were expressed by equal quantities (0.5 μg) of respective plasmids. Data are the mean ± s.e.m. of at four determinations. * p<0.05, two-tailed t test, versus its counterpart wild-type NA.

Discussion

The findings presented in this study detail the molecular mechanisms and temporal process by which the biologic fitness of oseltamtivir-resistant H1N1 viruses was altered by epistatic interactions between successive NA mutations. We discovered that the NA molecular changes (in the temporal order V234M-R222Q-K329E-D344N-H275Y-D354G) which occurred during NA genetic evolution during 1999~2009 gradually reduced and finally mitigated H275Y-generated defects on NA functionality. The genetic context-dependent fitness outcome of the H275Y-mutant viruses is characteristic of epistasis effects. Therefore, the sequential and cumulative epistatic interactions in NA contributed substantially to the evolutionary course of the H275Y mutation-containing H1N1 viruses.

Epistatic interactions occurs when the phenotypic outcome of a mutation is conditional to the presence or absence of other mutations in the genome[30] and it can have either positive (alleviating) or negative (aggravating) effects[30]. Epistasis in influenza A viruses has been estimated by modeling sequence evolution of human H3N2 and H1N1 viruses[31-33], and was suggested to be prevalent in HA and NA proteins[31]. There have, however, been few biological observations of epistatic interactions at the protein or virus phenotypic level. In one example, two advantageous immune-escape NP mutations that would have otherwise been destabilizing were acquired through interaction with other stabilizing mutations[34]. Our findings provide direct biological evidence of epistatic interactions in determining phenotypic outcome of a drug-resistant mutation during NA evolution. The interaction sites were concordant with a previous statistical analysis based on NA sequence alone[31], supporting a role for epistasis in shaping influenza virus evolution.

In the context of current knowledge, our results suggest a four-stage scenario of molecular evolution of the NA protein in H1N1 viruses during 1999–2009: 1) During 1999–2006, incidental introduction of the NA H275Y mutation into the predominant circulating NC99-like viruses severely impaired NA functionality and consequently virus fitness and ability to spread. These viruses became extinct. 2) Meanwhile, during 1999–2006, the NC99-like viruses gradually evolved under host pressures, including immunity. The NA protein sequentially acquired twelve substitutions, of which R222Q, K329E, G249K, and D344N have been shown to cause NA antigenic drift[35] while V234M, R222Q, K329E, and D344N greatly enhanced NA accumulation and/or substrate affinity. 3) In 2007, a new BR07 lineage, which had undergone hemagglutinin (HA) and NA antigenic drift, emerged and became prevalent. The introduction of the NA H275Y mutation into BR07-like viruses minimally reduced NA function due to its increased surface accumulation and/or substrate affinity as compared to the NC09-like NA. 4) After the H275Y substitution was acquired in mid-2007, final addition of the NA D354G fully compensated for the remaining H275Y-generated functional defects. Therefore, these BR07-like, oseltamivir-resistant viruses had no net loss of NA function, had biological fitness equal at least to that of the BR07-like wild-type viruses, and spread readily during the 2008–2009 influenza season.

Our findings show the presence of two types of epistatic interactions that collectively altered the H275Y phenotype. Four permissive epistatic mutations (V234M, R222Q, K329E, and D344N) were selected before addition of H275Y, increasing overall NA functionality fist and providing a permissive context for potential function-decreasing mutation. The D354G mutation, which was added to the BR07-like NA population immediately after H275Y acquisition, was very likely selected due to its ability to mitigate remaining H275Y-associated functional defects in a typical epistatic compensatory fashion. The combined phenotypic effects of both the permissive and compensatory mutations were required to maintain the net NA function of the H275Y-mutant protein.

The emergence of the V234M and R222Q substitutions before 1999 suggests that the selection forces for these substitutions had existed and was not drug pressure. The extensive overlap of the four permissive NA mutations (V234M, R222Q, K329E, and D344N) and the four reported antigenic-drift NA mutations (R222Q, G249K, K329E, and D344N)[35] suggests that host anti-NA immunity may have been the primary driver of the overlapping mutations. All four antigenic-drift mutations are located at the surface of the NA globular head, opposite the tetramer interface and viral membrane (Supplementary Fig. 1), potentially allowing interaction to both antibody and substrate. Our screening did not find the G249K substitution able to affect virus growth sufficiently to alter plaque size, but a small effect on enzyme affinity is possible. Thus, selection of antigenic-drift NA mutations likely inadvertently provided a permissive context for the H275Y mutation and altered its phenotypic outcome. The antigenic drift-driven NA mutation spread can be a plausible explanation in addition to previously proposed genetic hitchhiking process[36]. The non-overlapping mutation V234M, not located at the NA surface, did not affect enzyme affinity and was not a likely antigenic drift mutation. However, it preceded a series of NA affinity-changing and antigenic-drift mutations, showing that the first altered NA phenotype was expression level and the second was functionality, suggesting that protein expression level could affect introduction of function-changing mutations during protein evolution.

Our mechanistic reconstruction of the evolutionary course of the H1N1 viruses has direct implications for assessing the potential fitness, evolution and spread of other H275Y-mutant viruses of N1-subtype. In the case of 2009 pdmH1N1 viruses, the severe defects on NA functionality generated by the H275Y mutation is consistent with the reduced transmissibility of H275Y mutants experimentally[37, 38] and the limited detection of such viruses during surveillance (1%~3.5% detection rate)[39-41]. However, experimental transmission and community spread of H275Y-mutant pdmH1N1 viruses has been reported[40, 42–44], including a very recent alarming community spread in Sapporo, Japan[45], suggesting that potential permissive mutation might have been emerged in certain cluster of variants. Three studies using computational predictions or in vitro and in vivo experiments have suggested several potential permissive NA mutations for the H275Y in pdmH1N1 viruses[46-48], which are entirely different from the one identified in seasonal H1N1 viruses. Given that NA genes of the seasonal H1N1 and pdmH1N1 were phylogenetically distant and had extensive amino acid dissimilarities, it is not surprising that different H275Y-permissive mutations would be required for pdmH1N1 viruses, which is also very likely the case for other N1-subtype viruses. Our study suggest that the evolution course of pdmH1N1 NA toward the potential resistance spread should be monitored longitudinally by assessing the extent of H275Y-generated NA defects, especially in community-transmitted resistant isolates. The NA defects evaluated by a spectrum of assays from protein expression, enzyme activity to virus transmission level would be most informative for assessing dominance potential of resistant viruses. The occurrence of smaller or no NA defects in all assay levels should worth a high alert.

Similarly, the deleterious outcome of the H275Y mutation on the NA protein from H5N1 viruses is consistent with the rare case of resistant viruses detected in the field (only 3 cases where drug treatment or prophylaxis was used)[49, 50]. Because the highly-pathogenic H5N1 (hpH5N1) viruses carry unique virulence determinants[51], H275Y-generated NA defects may not alter their overall pathogenicity. However, as hpH5N1 infection is a substantial public health concern, it would be advisable to compare the biological fitness of wild-type vs. emerging resistant viruses by identifying H275Y-generated NA defects. Of special value in our study is the use of NA-containing plasmids to assess defects in NA protein expression and activity, as described here, averting the need to engineer infectious hpH5N1 virus.

Our findings of the H275Y mutation in H1N1 viruses may help in assessing the biological outcome of other NA inhibitor–resistance markers. NA inhibitor–resistance mutations are virus subtype- and drug-dependent and may be NA catalytic residues (such as NA 118, 276, 292, etc., N2 numbering) or framework residues (such as NA 119, 222, 274, 294, etc., N2 numbering)[7, 52]. The H275Y (274 in N2 numbering) is a framework residue mutation[52]. To date, surveillance and clinical observation have shown that R292K, E119V, N294S, and I222V mutant H3N2 viruses[4, 7–9], I222V and I222R mutant pandemic H1N1 viruses[40, 53, 54], and R292K mutant H7N9 viruses [55] have emerged only from drug-treated individuals and have not been able to spread at the population level. Their compromised epidemiologic fitness indicates that the biological outcome of these mutations on the proteins and viruses remains detrimental. Among these, resistance mutations at the catalytic residues presumably would cause much more profound loss of NA functionality than framework mutations, and could consequently completely eliminate the virus biological fitness. However, the smaller loss of function by framework mutation might be restored more readily, like the H275Y. Another two framework mutations I222V and I222R in the pandemic H1N1 viruses have shown less loss of NA function[56, 57], and its NA evolution course may therefore merit more attention.

Interaction of the H275Y mutation with permissive and compensatory NA mutations explains the restoration of the NA function and the fitness of the resistant viruses to the wild-type level, but NA function alone does not explain why the resistant viruses outtransmitted over the wild-type viruses. Therefore, the roles of other viral proteins in conferring the greater fitness of resistant viruses should be examined. For example, HA-binding to receptors also plays important role in influenza virus transmission. The NA and HA genes of H1N1 viruses showed the same genetic clustering and antigenic lineage patterns[16-19], indicating their co-evolution. The HAs of the resistant viruses exhibited reduced reactivity to antibody inhibition[58]. HA mutations can promote replication of BR07-like H275Y-mutant viruses in cells[59] and drive NA mutations that alter antigenicity and NA-inhibitor susceptibility[60]. All of these findings suggest the occurrence of phenotypic interactions between the HA and NA proteins (probably to promote functional balance) during the evolution of H275Y-mutant viruses. Comparison of the HA sequences of the BR07-like viruses used in the ferret experiment, revealed four HA mutations (G185V, N186D, A189T and H192R) shared by the H275Y mutant viruses (Supplementary Table 5). These four positions are located in or near to the 190 helix which is one of the three structure elements enclosing the base of the HA binding site[61]; although these positions are not directly involved in receptor binding, whether and how they can affect the receptor binding and virus transmissibility would be of interest for further study. In addition, intra-subtype, inter-clade segment reassortment events observed in the BR07-like resistant viruses[18, 62] suggest that a combination of internal protein interactions plays a role in enhancing overall virus fitness.

In conclusion, we have shown how the fitness cost of a drug-resistance NA mutation in human seasonal influenza virus can be reduced or eliminated during NA evolution, allowing the resistant strain to spread. As influenza viruses continue to evolve under immune pressure and NA inhibitors remain the primary therapeutic option, we suggest that NA genetics, antigenic drift, and resistance mutation–generated NA functional outcome should be collectively and continually monitored to assess the biologic fitness of emerging resistant viruses and their potential evolutionary course toward resistance spread.

Methods

Cells and viruses

Pandemic DM/524/09 and DM/528/09 viruses were provided by Statens Serum Institute, Copenhagen, Denmark. NY/3467/09 and NY/1692/09 viruses were provided by the Wadsworth Center of the New York State Department of Health. Wild-type and NA H275Y-mutant MS/03/02 viruses were provided by the Neuraminidase Inhibitor Susceptibility Network. The remaining H1N1 viruses, including CA/04/09, BR/59/07, HW/28/07, NJ/15/07, GA/17/06, GA/20/06, TX/30/07 and MEM/13/06, were provided by the U.S Centers for Disease Control and Prevention. The NA gene segments of VN/1203/04 and HK/213/03 H5N1 viruses were generated into plasmids previously[63, 64], and no infectious H5N1 viruses were used in the present study. All virus stocks used underwent an initial limited number (1~2) of passages in MDCK cells upon reception, to maintain their original properties. The 50% tissue culture infectious dose (TCID50) or plaque forming unit (pfu) was used to measure virus infectivity in MDCK cells as previously described[37].

Sequence analysis

Viral RNA was isolated by using the RNeasy Mini kit (Qiagen). Samples were reverse-transcribed and analyzed by PCR, using segment-specific primers as described previously[65]. Sequencing was performed by the Hartwell Center for Bioinformatics and Biotechnology at St. Jude Children’s Research Hospital. DNA sequences were completed and edited by using the Lasergene software package (DNASTAR).

Phylogenetic estimation of timing and sequence of mutation events

All the published NA gene sequences of seasonal H1N1 viruses were downloaded from NCBI (Influenza Virus Resource). Initial visual inspection of multiple NA sequence alignments containing taxa with and without the H275Y oseltamivir resistance mutation identified 11 amino acid mutations associated with the NA lineage divergence. In order to estimate the time of most recent common ancestor that possessed the mutations of interest, we used all available full-length H1N1 NA sequences that also had exact date of isolation in the accession, and generated three datasets comprised of randomly selected taxa plus the past vaccine strains, reference sequences and six available pre-2007/2008 H275Y variants (ntax=413 for each dataset). The SRD06 approximate codon model of evolution with an HKY85+gamma nucleotide substitution model for codon position 1+2 and codon partition 3 was used to estimate the phylogenetic tree. The ancestral state reconstruction of all tested mutations was jointly estimated with the phylogeny. The NA evolutionary history and the ancestral states reconstruction were estimated using a tip date calibrated Bayesian relaxed-clock phylogenetic tree reconstruction method[27]. We assigned the individual mutations as a tip observed state[66]. These mutations were indicators for possible adaptive mutations. We assessed which mutations were fixed in the population, the sequence of fixation and when this fixation occurred. Phylogenetic analysis was carried out in BEAST (v1.7.3).

Plasmid construction and site-directed mutagenesis

The full-length cDNA (8 segments) of BR/59/07 virus and the NA segments of GA/17/06 and GA/20/06 viruses were cloned into the pHW2000 plasmid vector to generate reverse-genetics viruses, as described previously[67]. Mutations of interest were introduced into the corresponding NA genes by using QuickChange site-directed mutagenesis (Stratagene) and confirmed by Sanger sequencing. The full-length protein-coding sequences of the NA genes of interest were cloned into the pCAGGS plasmid vector for expression of recombinant NA proteins in BHK cells, as described previously[68]. The HA epitope tag (YPYDVPDYA) was added to the N or C terminus of the NA protein by incorporating the coding sequences into the primers.

Generation of reverse-genetics influenza viruses

Reverse-genetics viruses were rescued by transfecting a 293T/MDCK cell co-culture with eight PHW2000 plasmids containing the eight virus segments, using TransIT LT-1 (Panvera), as described previously[67]. Rescued viruses contained variant NA segments in the uniform backbone from other seven segments of the BR/59/07 virus. The rescued supernatants from co- culture were first titrated by plaque assay and then one passage of the supernatant in MDCK cells (MOI = 0.001~0.01) was used to prepare stocks of all reverse-genetics viruses.

Generation and detection of recombinant NA proteins

BHK cells were transfected with the respective NA gene-containing pCAGSS plasmids by using lipofectamine 2000 (Invitrogen), as described previously[68], and NA protein expression was measured by NA enzyme activity assay, western blot, or flow cytometry. For comparison of different proteins, BHK cells were transfected with equal quantities of the respective pCAGGS plasmids (0.5 μg unless indicated otherwise) in a 12-well plate, in parallel. For western blot, 0.05 μg of pCMV-GFP plasmid was co-transfected with the respective pCAGGS plasmid. Cells were harvested in approximately 36 h post-transfection. For the subsequent NA enzyme activity assay, cells were lysed with non-denaturing NA enzyme buffer containing 0.5% Triton X-100. For western blot, cells were lysed with denaturing SDS-containing RIPA buffer. All samples were processed in parallel. For protein expression level determination, anti-HA epitope monoclonal antibodies (Sigma, 1:20,000) were used to detect the epitope-tagged recombinant proteins in western blot or flow cytometry. In western blot, β actin was detected (anti-β actin, Santa Cruz, 1:20,000) as loading control, and GFP (detected by anti-GFP, Santa Cruz, 1:2,000) was used as transfection efficiency control. Fluorescent dye-labeled (IRDye680 or IRDye800) anti-mouse or –rabbit secondary antibodies (LI-COR, 1:20,000) were used to detect the respective primary antibodies. The blot signals were visualized and analyzed in Odyssey imaging system (LI-COR).

NA enzyme activity and kinetics

The NA activity assay was based on a fluorometric assay using the substrate 2′-(4-methylumbelliferyl)-α-D-N-acetylneuraminic acid (MUNANA) (Sigma-Aldrich) in enzyme buffer containing 33 mM 2-(N-morpholino) ethanesulfonic acid hydrate (Sigma-Aldrich), 4 mM CaCl2 at pH 6.5 at 37°C for 30 mins[69]. Hydrolysis of 1 mol MUNANA by NA liberated 1 mol fluorescent product 4-methylumbelliferone[69] and the fluorescence of liberated product was read in a Synergy 2 multi-mode microplate reader (BioTek), using excitation and emission wavelengths of 360 and 460 nm, respectively. Purified 4-methylumbelliferone (Sigma-Aldrich) was used as a standard for estimation of product generation and substrate consumption in the assay (e.g. Supplementary Fig. 8a). For NA enzymatic activity measurement, all the harvested samples were serially diluted and were assayed in parallel with 100 μM MUNANA, and then NA quantity-activity curves were generated. The dilution at which all NAs were within their linear range of product generation and approximately consumed less than 20% of substrate was adopted for activity measurement (e.g. Supplementary Fig. 8b, c). Then the activity was normalized to a reference NA as the % relative activity (e.g. Fig. 2a). The enzyme activity of all expressed NA proteins was measured by the same process.

For NA kinetics assay, all viruses were standardized to a dose of 106.0 PFU/ml, a dose at which all viruses were determined to consume less than 10% of 100 μM MUNANA substrate in NA activity assays (e.g. Supplementary Fig. 8a). We then measured NA enzyme kinetics in the enzyme buffer using serial concentrations of MUNANA (final concentration, 800-6.25 μM, 2-fold serial dilution) in a total volume of 50 μl. The fluorescence of released 4-methylumbelliferone was measured at 37°C every 60 sec for 40 mins. Km was calculated by fitting the data to the Michaelis-Menten equations by nonlinear regression in Prism5 software (GraphPad).

Transmission experiments in ferrets

All experiments were conducted in an ABSL2+ laboratory under applicable laws and guidelines and after approval from the St. Jude Children’s Research Hospital Animal Care and Use Committee. Four- to 5-month-old ferrets were obtained from Triple F farm; all were tested sero-negative for contemporary influenza A H1N1 and H3N2 viruses and influenza B viruses. Ferrets were initially housed in isolators and monitored for 3–5 days to establish baseline body temperature and overall health. Donor and contact ferrets were housed separately. Transmission experiments (1 donor + 1 direct-contact recipient + 1 respiratory-droplet recipient) were conducted in triplicate for each virus. Three donor ferrets were inoculated intranasally with 105 TCID50 of virus in 1.0 ml sterile PBS. At day 1 p.i., each donor was co-housed with 1 naïve direct-contact ferret. One additional naïve ferret was placed in an adjacent cage separated by double-layered (~5 cm apart) perforated dividers to assess respiratory-droplet transmission; air flow was directed from the donor cage toward the recipient cage using a Borazine gun (Zero Toys). Nasal wash samples were collected at days 1 and 2 p.i and then every other day for 14 days by flushing both nostrils with 1.0 ml PBS. TCID50 titers of nasal washes were determined in MDCK cells. Serum samples were collected 3 weeks after inoculation, treated with receptor-destroying enzyme, heat-inactivated at 56 °C for 30 min, and tested against the homologous virus by HI assay with 0.5% packed turkey red blood cells.

Statistical analysis

All statistical analyses were performed using student t test in Prism5 software. Data are presented as mean ± s.e.m. The level of significance was determined as P<0.05.