Research Article

Antigenic evolution of human influenza H3N2 neuraminidase is constrained by charge balancing

Department of Biochemistry, University of Illinois at Urbana-Champaign, United States
Department of Physics, University of Washington, United States
Max Planck Institute for Dynamics and Self-Organization, Germany
Fred Hutchinson Cancer Research Center, United States
Center for Biophysics and Quantitative Biology, University of Illinois at Urbana-Champaign, United States
Carl R. Woese Institute for Genomic Biology, University of Illinois at Urbana-Champaign, United States
Carle Illinois College of Medicine, University of Illinois at Urbana-Champaign, United States

Dec 8, 2021

https://doi.org/10.7554/eLife.72516

Version of Record

Accepted for publication after peer review and revision.

Download
Cite
Share
CommentOpen annotations (there are currently 0 annotations on this page).

Version of Record published: December 17, 2021 (This version)
Accepted Manuscript published: December 8, 2021 (Go to version)
Accepted: December 7, 2021
Received: July 27, 2021
Preprint posted: July 12, 2021 (Go to version)

1. Of interest
Native American genetic ancestry and pigmentation allele contributions to skin color in a caribbean population

Khai C Ang, Victor A Canfield ... Keith C Cheng

Research Article Jun 9, 2023
Further reading

Article
Figures and data
Abstract
Editor's evaluation
Introduction
Results
Discussion
Materials and methods
Data availability
References
Decision letter
Author response
Article and author information
Metrics

Abstract

As one of the main influenza antigens, neuraminidase (NA) in H3N2 virus has evolved extensively for more than 50 years due to continuous immune pressure. While NA has recently emerged as an effective vaccine target, biophysical constraints on the antigenic evolution of NA remain largely elusive. Here, we apply combinatorial mutagenesis and next-generation sequencing to characterize the local fitness landscape in an antigenic region of NA in six different human H3N2 strains that were isolated around 10 years apart. The local fitness landscape correlates well among strains and the pairwise epistasis is highly conserved. Our analysis further demonstrates that local net charge governs the pairwise epistasis in this antigenic region. In addition, we show that residue coevolution in this antigenic region is correlated with the pairwise epistasis between charge states. Overall, this study demonstrates the importance of quantifying epistasis and the underlying biophysical constraint for building a model of influenza evolution.

Editor's evaluation

This paper presents a systematic analysis of the fitness landscape of the influenza virus protein neuraminidase (NA). The authors generate 864 different combinations of amino acids at seven positions in six genetic backgrounds sampled 10 years apart and measure the fitness of the resulting virus. This fitness landscape is characterized by strong epistatic interactions, including a strong tendency to maintain the local charge of the protein. Such systematic characterizations of important proteins of viral pathogens are crucial to develop principled models to understand and predict their evolution.

https://doi.org/10.7554/eLife.72516.sa0

Introduction

There are two major antigens on the surface of influenza virus, hemagglutinin (HA) and neuraminidase (NA). Although influenza vaccine development has traditionally focused on HA, NA has emerged as an effective vaccine target in the past few years because recent studies have shown that NA immunity has a significant role in protection against influenza infection (Monto et al., 2015; Weiss et al., 2020; Couch et al., 2013; Memoli et al., 2016; Krammer et al., 2018). Influenza NA has an N-terminal transmembrane domain, a stalk domain, and a C-terminal head domain. The head domain of NA functions as an enzyme to cleave the host receptor (i.e., sialylated glycan), which is essential for virus release. Most NA antibodies target the surface loop regions that surround the highly conserved catalytic active site (Gulati et al., 2002; Malby et al., 1994; Venkatramani et al., 2006; Colman et al., 1983; Air, 2012; Zhu et al., 2019; Tulip et al., 1992). Due to the need to constantly escape from herd immunity (also known as antigenic drift), both HA and NA of human influenza virus have evolved extensively (Kilbourne et al., 1990; Sandbulte et al., 2011; Westgeest et al., 2015). For example, since influenza H3N2 virus entered the human population in 1968, its HA and NA have accumulated more than 83 and 73 amino acid mutations, respectively (Figure 1—figure supplement 1), which accounted for ~15% of their protein sequences. However, the evolution of NA is much less well characterized as compared to HA.

To understand how the evolutionary trajectories of NA are being shaped, it is important to characterize the underlying biophysical constraints that govern the fitness of individual amino acid mutations and epistatic interactions between mutations (Starr and Thornton, 2016; Echave and Wilke, 2017; Gong et al., 2013). Epistasis is a phenomenon in which the fitness effect of a mutation is dependent on the presence or absence of other mutations. Since epistasis can lead to differential fitness effects of a given mutation on different genetic backgrounds, it can restrict evolutionary trajectories or open up a new functional sequence space that would otherwise be inaccessible (Wu et al., 2016). As a result, epistasis is a primary challenge for predicting evolution (Miton and Tokuriki, 2016; Luksza and Lässig, 2014). Nevertheless, epistasis is pervasive in natural evolution in general (Breen et al., 2012) and has been shown to influence the evolution of influenza virus (Gong et al., 2013; Wu et al., 2016; Lyons and Lauring, 2018; Wu et al., 2018; Kryazhimskiy et al., 2011; Wu et al., 2020; Koel et al., 2019; Bloom et al., 2010; Abed et al., 2011; Wu et al., 2013; Duan et al., 2014). While epistasis is critical for the emergence of oseltamivir-resistant mutants of influenza NA (Bloom et al., 2010; Abed et al., 2011; Wu et al., 2013; Duan et al., 2014), the role of epistasis and the underlying biophysical constraints on NA antigenic evolution remains largely unclear.

Deep mutational scanning combines saturation mutagenesis and next-generation sequencing to determine the phenotypic effects of numerous mutations in a highly parallel manner (Lee et al., 2018; Narayanan and Procko, 2021; Fowler and Fields, 2014). Deep mutational scanning has been employed to measure the replication fitness effect of all possible single amino acid mutations across different influenza proteins (Lee et al., 2018; Thyagarajan and Bloom, 2014; Doud et al., 2015; Soh et al., 2019; Hom et al., 2019), which in turn can help to model the natural evolution of human influenza virus (Lee et al., 2018; Thyagarajan and Bloom, 2014). However, most deep mutational scanning experiments lack the power to systematically measure the fitness of variants with two or more mutations, hence epistasis. More recently, combinatorial mutagenesis was used in conjunction with next-generation sequencing to examine the local fitness landscape and to identify epistasis in influenza HA (Wu et al., 2018; Wu et al., 2020; Wu et al., 2017). Unlike saturation mutagenesis in conventional deep mutational scanning, combinatorial mutagenesis enables us to measure the fitness of high-order mutants. By dissecting the mechanistic basis of epistasis, these studies have provided important insight into the biophysical constraints on HA antigenic evolution (Wu et al., 2018; Wu et al., 2020).

In this study, we aim to understand how epistatic effects influence NA antigenic evolution and investigate the underlying biophysical constraints. Specifically, we coupled combinatorial mutagenesis and next-generation sequencing to characterize the local fitness landscape of an NA antigenic region in six different human H3N2 strains that were isolated ~10 years apart. Our results indicate that the local fitness landscape of this NA antigenic region is highly correlated across six different genetic backgrounds. In-depth analyses further demonstrate that local net charge balancing is a biophysical constraint that governs the epistasis within this NA antigenic region. Lastly, we show that epistasis is correlated with residue coevolution in naturally circulating influenza strains.

Results

Local fitness landscape of an antigenic region on NA

When the structure of N2 NA was first reported in 1983 (Varghese et al., 1983), seven regions (I–VII) were proposed to be targeted by antibodies (Colman et al., 1983). Within regions I–III, mutations at residues 329, 344, 368, and 370 have been shown to escape monoclonal antibodies (Gulati et al., 2002; Varghese et al., 1988; Air et al., 1985). These four residues, along with residues 328 (region I), 367 (region III), and 369 (region III), form a cluster of seven residues that are very close in space (Figure 1a and b). Both residues 328 and 367 are under positive selection in human H3N2 NA (Westgeest et al., 2012), whereas coevolution of residues 367 and 369 in human H3N2 NA has created an N-glycosylation site (NXT) at residue 367 during 2010 (Figure 1c). These observations indicate that residues 328, 367, and 369 also participated in the antigenic drift of human H3N2 NA. By focusing on this seven-residue antigenic region, this study aimed to dissect the biophysical constraints on NA antigenic evolution.

Figure 1 with 3 supplements see all

Download asset Open asset

Natural evolution of an antigenic region in human H3N2 neuraminidase (NA).

(a) The enzyme active site and the antigenic region of interest are highlighted as red and blue, respectively, on one protomer of the NA tetramer. The sialic acid within the active site is shown in yellow. (b) Seven residues in the antigenic region of interest are highlighted in blue spheres. (c) Natural occurrence frequencies of the amino acid variants that have a natural occurrence at >50% in any given year at the residues of interest are shown (Figure 1—source data 1 and Figure 1—source data 2). Of note, only those amino acid variants with a natural occurrence at >80% in any given year were included in our mutant libraries. Therefore, although D368 is shown in this plot, it was not included in our mutant libraries since it only reached a maximum occurrence of 67%.

Figure 1—source data 1 Human H3N2 neuraminidase (NA) protein sequences.: https://cdn.elifesciences.org/articles/72516/elife-72516-fig1-data1-v2.zip
Download elife-72516-fig1-data1-v2.zip
Figure 1—source data 2 Natural occurrence frequencies of the major amino acid variants at the residues of interest in different years.: https://cdn.elifesciences.org/articles/72516/elife-72516-fig1-data2-v2.zip
Download elife-72516-fig1-data2-v2.zip

We first compiled a list of amino acid variants at NA residues 328, 329, 344, 367, 368, 369, and 370 that reached an occurrence frequency of >80% in any given year during the natural evolution of human H3N2 viruses since 1968 (Figure 1c). This list includes two amino acid variants at residue 328 (Asn and Lys), three at residue 329 (Asn, Ser, and Asp), three at residue 344 (Lys, Glu, and Arg), three at residue 367 (Asn, Ser, and Gly), two at residue 368 (Lys and Glu), four at residue 369 (Lys, Glu, Asp, and Thr), and two at residue 370 (Leu and Ser). Together, there are 2 × 3 × 3 × 3 × 2 × 4 × 2 = 864 possible amino acid combinations (also called haplotypes) across these seven residues, although only 53 of them have been observed in naturally circulating human H3N2 strains (Figure 1—figure supplement 2). Using PCR primers that carry degenerate nucleotides (see Materials and methods and Supplementary file 1), these 864 variants were introduced into the NA of six different strains (genetic backgrounds) from 1968 to 2019, namely, A/Hong Kong/1/1968 (HK68), A/Bangkok/1/1979 (Bk79), A/Beijing/353/1989 (Bei89), A/Moscow/10/1999 (Mos99), A/Victoria/361/2011 (Vic11), and A/Hong Kong/2671/2019 (HK19). All these six strains were historical vaccine strains and isolated approximately 10 years apart.

To measure the virus replication fitness of all 864 variants in the six different genetic backgrounds of interest (i.e., HK68, Bk79, Bei89, Mos99, Vic11, and HK19), we employed a high-throughput experimental approach that coupled combinatorial mutagenesis and next-generation sequencing. Unlike conventional deep mutational scanning, which studies all possible single amino acid mutations across a protein or domain of interest, our approach analyzes all possible combinations of a subset of mutations. Subsequently, six different local fitness landscapes, each with 864 variants, were obtained. Fitness value of each variant, which was defined in log scale (see Materials and methods), was normalized to the corresponding wild type (WT), such that the WT fitness value was 0, whereas positive and negative fitness values represented beneficial and deleterious variants, respectively. Two biological replicates of each experiment were performed. Overall, a high correlation was observed between the replicates (Pearson correlation = 0.81–0.89) except for HK68 and Bk79 (0.39 and 0.63, respectively), mostly due to the high measurement noise for low fitness variants (Figure 1—figure supplement 3).

Dynamics of the local fitness landscape across genetic backgrounds

To examine whether the local fitness landscape of the NA antigenic region of interest changes over time, fitness measurements of the 864 variants were compared across different genetic backgrounds (Figure 2a and b). Interestingly, while most variants were strongly deleterious (i.e., had a fitness of <-1) in HK68, variants with a fitness of <-1 were rare in other genetic backgrounds (Figure 2a). Notably, the difference in variant fitness distribution across genetic backgrounds was not due to the difference in their WT replication fitness since the viral titers from a rescue experiment were similar among HK68 WT, Bei89 WT, and Mos99 WT (Figure 2—figure supplement 1). In contrast, the fitness of individual variants correlated well across genetic backgrounds (Pearson correlation = 0.48–0.79, Figure 2b), despite their differences in variant absolute fitness values (Figure 2a). These results demonstrate that although epistasis exists between the seven-residue antigenic region and the rest of the NA sequence, such epistasis is largely variant-nonspecific. Consequently, the topology of the local fitness landscape is highly conserved across genetic backgrounds. This observation is very different from a similar study on a major antigenic site of HA, where the topology of the local fitness landscape differed dramatically among genetic backgrounds (Wu et al., 2020).

Figure 2 with 1 supplement see all

Download asset Open asset

Comparing the local fitness landscapes of the neuraminidase (NA) antigenic region in six human H3N2 strains.

(a) Variant fitness distributions in different genetic backgrounds are shown as a sina plot (Figure 2—source data 1). Each data point represents one variant. A total of 864 data points are present for each row (each genetic background). (b) Correlations of variant fitness distribution among different genetic backgrounds are shown, with each data point representing one variant. Pearson correlation coefficients (R) are indicated. (**a, b**) Data points corresponding to the wild type (WT) sequences of HK68, Bk79, Bei89, Mos99, and HK19 are colored as indicated. Of note, the WT sequence of Vic11 contains a naturally rare variant T329. Therefore, the WT sequence of Vic11 was not included in our mutant libraries. (c) Naturally occurring variants were grouped by the year of isolation, and their mean fitness in different genetic backgrounds is shown (Figure 2—source data 2). Different genetic backgrounds are represented by different lines, which are color coded as indicated on the right. Error bars represent the standard error of mean. This analysis included 66,562 NA sequences from human H3N2 strains that were isolated between 1968 and 2020 (Figure 1—source data 1).

Figure 2—source data 1 Fitness value of each variant across six different genetic background.: https://cdn.elifesciences.org/articles/72516/elife-72516-fig2-data1-v2.xlsx
Download elife-72516-fig2-data1-v2.xlsx
Figure 2—source data 2 The yearly mean fitness of naturally occurring variants in different genetic backgrounds.: https://cdn.elifesciences.org/articles/72516/elife-72516-fig2-data2-v2.zip
Download elife-72516-fig2-data2-v2.zip

We further examined the fitness of naturally occurring variants in different genetic backgrounds (HK68, Bk79, Bei89, Mos99, Vic11, and HK19) (Figure 2c). Most natural variants have a fitness between –0.5 and 0.5 on the background of Bk79, Bei89, Mos99, Vic11, and HK19. In contrast, most natural variants from 1980s onward are strongly deleterious with fitness <-1 on the background of HK68, indicating many natural variants would not have emerged if the genetic background had not evolved. In other words, the emergence of natural variants in the antigenic region of interest is contingent on the evolution of the rest of the NA sequence. This result is consistent with the observation that HK68 has a much lower mutational tolerance at this antigenic region (Figure 2a).

Conservation of pairwise epistasis across genetic backgrounds

To understand the biophysical constraints on NA antigenic evolution, the local fitness landscapes were decomposed into additive fitness effects of individual amino acid variants and pairwise epistasis between amino acid variants (Starr and Thornton, 2016). Additive fitness describes the independent contributions of each amino acid variant to fitness, whereas pairwise epistasis describes the nonadditive interactions between amino acid variants. For each genetic background, additive fitness and pairwise epistasis were inferred from the variant fitness data using an established statistical learning model (Tareen et al., 2020) (see Materials and methods). The model was evaluated using repeated k-fold cross-validation and hyperparameters were chosen by maximizing the R² of model prediction and the Pearson correlation coefficient of model parameters (Figure 3—figure supplement 1). While the correlations between additive fitness contributions varied hugely across genetic backgrounds (Pearson correlation = 0.12–0.88, Figure 3a and c, Figure 3—figure supplement 2), we observed generally strong correlations between pairwise epistatic effects across the six different genetic backgrounds (Pearson correlation = 0.69–0.86, Figure 3b and c, Figure 3—figure supplement 3). Of note, the variation of additive fitness contributions across genetic backgrounds does not seem to strictly depend on the similarity between genetic backgrounds since the correlation between the additive fitness contributions of HK68 and HK19 (Pearson correlation = 0.46) is much higher than that of HK68 and Bk79 (Pearson correlation = 0.12), which have a shorter time separation. This observation points at a possibility for complex epistatic interactions between the antigenic region of interest and other regions on NA. Overall, these results indicate that pairwise epistasis, but not additive fitness, is highly conserved at the NA antigenic region of interest across genetic backgrounds. Consistently, an ‘additive-only’ model, without accounting for epistasis, shows a poor fit to the fitness landscape data compared to the model above with epistasis (Figure 3—figure supplement 4a), despite the fact that the inferred additive fitness effects correlate well between the two models (Figure 3—figure supplement 4b).

Figure 3 with 5 supplements see all

Download asset Open asset

Inference of additive fitness and pairwise epistasis.

(a) Parameters for additive fitness in different genetic backgrounds are shown. (b) Pairwise epistatic effects show strong correlations among six different genetic backgrounds, and therefore, Bk79 is shown as a representative. The identity of a given double amino acid variant is represented by the labels on the x- and y-axes. The same plots for other genetic backgrounds are shown in Figure 3—figure supplement 5. (**a, b**) Positive additive fitness and pairwise epistasis are in red. Negative additive fitness and pairwise epistasis are in blue. The magnitude is proportional to the size of the circle. (c) Correlation matrices of additive fitness and pairwise epistasis among six genetic backgrounds are shown as a heatmap. See Figure 3—figure supplements 4 and 5 for the related scatter plots.

Local net charge imposes a biophysical constraint on NA antigenic evolution

Next, we investigated the biophysical constraints that have led to such conserved patterns of epistatic interactions in the NA antigenic region of interest. We noticed that amino acid variants with opposite charges usually exhibited positive epistasis, whereas amino acid variants with the same charge usually exhibited negative epistasis (Figures 3b and 4a, Figure 3—figure supplement 5, Figure 4—figure supplement 1; see Materials and methods). In contrast, pairwise interactions that involved a neutral amino acid variant did not show any bias towards positive or negative epistasis. These results suggest that balancing of charged amino acid variants imposes a key biophysical constraint on the evolution of this NA antigenic region.

Figure 4 with 4 supplements see all

Download asset Open asset

The importance of local net charge in the neuraminidase (NA) antigenic region.

(a) Pairwise epistasis in genetic background Bk79 is shown and classified based on amino acid charges. +chg represents positively charged amino acids (K/R), -chg represents negatively charged amino acids (D/E), and neutral represents the remaining amino acids. The same plots for other genetic backgrounds are shown in Figure 4—figure supplement 1. (b) NA structure from H3N2 A/Memphis/31/98 (PDB 2AEP), which has K328, E344, and K369, was analyzed (Gulati et al., 2002). A similar NA structure from H3N2 A/Tanzania/205/2010, which has K328 and E344, is shown in Figure 4—figure supplement 2. (c) The relationship between variant fitness and net charge in genetic background Bk79 is shown. A smooth curve was fitted by loess and shown in teal. The same plots for other genetic backgrounds are shown in Figure 4—figure supplement 4. (**d–f**) Relationship between the side-chain–side-chain distances (Figure 4—source data 1) and epistasis for (d) variant pairs with opposite charges, (e) variant pairs with the same charge, and (f) variant pairs that involve a neutral amino acid in genetic background Bk79 is shown. The same plots for other genetic backgrounds are shown in Figure 4—figure supplement 3. Pearson correlation coefficient (R) is indicated.

Figure 4—source data 1 Pairwise side-chain–side-chain distances within neuraminidase (NA) antigenic region.: https://cdn.elifesciences.org/articles/72516/elife-72516-fig4-data1-v2.zip
Download elife-72516-fig4-data1-v2.zip

To further probe the mechanism of epistasis in this NA antigenic region, we analyzed a published crystal structure of human H3N2 NA that has K328, E344, and K369 (Figure 4b; Venkatramani et al., 2006). Both variant pairs E344/K369 and K328/E344 exhibited positive epistasis across all genetic backgrounds, whereas K328/K369 exhibited negative epistasis (Figure 3—figure supplement 5). In fact, E344 and K369 had the strongest positive epistasis in five of the six genetic backgrounds of interest (except Bei89). Our structural analysis showed that the distance between the side-chain carboxylate oxygen of E344 and the side-chain amine nitrogen of K369 is 6.6 Å, whereas the distance between the side-chain carboxylate oxygen of E344 and the side-chain amine nitrogen of K328 is 5.5 Å (Figure 4b, Figure 4—figure supplement 2). At these distances, salt bridges cannot be formed and the electrostatic attraction force is negligible (Kumar and Nussinov, 2002; Shashikala et al., 2019). Similarly, the distance between the side-chain amine nitrogen atoms of K328 and K369 is 11.6 Å (Figure 4b), which is too far for any significant electrostatic repulsion force (Shashikala et al., 2019). As a result, direct side-chain–side-chain interaction via electrostatic attraction or repulsion is unlikely to be a major determinant for epistasis in this NA antigenic region. Consistently, pairwise epistasis between two amino acid variants does not correlate with their side-chain–side-chain distances (Figure 4d–f, Figure 4—figure supplement 3), further substantiating that direct side-chain–side-chain interaction is not a determinant for epistasis here.

We then analyzed the relationship between variant fitness and local net charge as a molecular phenotype. Here, local net charge was defined as the sum of charges at the seven residues of interest (residues 328, 329, 344, 367, 368, 369, and 370), where positively charged amino acids (K and R) were +1 and negatively charged residues (D and E) –1. In all genetic backgrounds, variants tended to have a higher fitness when the local net charge was around –1, while variants with a more positive or negative local net charge usually had a lower fitness (Figure 4c, Figure 4—figure supplement 4). Overall, our results demonstrate that the local net charge is a key molecular phenotype with a biophysical function that is under balancing selection. This biophysical phenotype imposes a strong constraint on the evolution of the NA antigenic region, reflected in the conserved epistatic interactions between amino acid variants of this region.

Impact of local net charge and epistasis on NA evolution

Next, we aimed to understand whether the local net charge at the NA antigenic region of interest influences its evolution in circulating human H3N2 strains. We retrieved 66,562 human H3N2 NA sequences spanning from 1968 to 2020 from the Global Initiative for Sharing Avian Influenza Data (GISAID) (Shu and McCauley, 2017). Most natural variants in the antigenic region of interest had a local net charge between –1 and +1 (Figure 5a). Natural variants with a local net charge of −3,–2, + 2, or +3 could also be observed but were rare. This observation suggests that the natural evolution of this NA antigenic region is constrained by balancing the local net charge.

Figure 5 with 6 supplements see all

Download asset Open asset

Correlation between epistasis and coevolution of charge states in the neuraminidase (NA) antigenic region.

(a) The natural evolution of frequency of haplotypes by local net charge at the NA antigenic region of interest is shown (Figure 5—source data 1). Key mutation events that led to changes in local net charge are indicated. (b) Pairwise epistasis of charge states among different residues. Amino acids were classified based on charges. (+) represents positively charged amino acids (K/R). (-) represents negatively charged amino acids (D/E). (n) represents the remaining amino acids. (c) Relationship between coevolution score (Figure 5—source data 2) and pairwise epistasis in genetic background Bk79 is shown. Pearson correlation coefficient (R) is indicated. The same plots for other genetic backgrounds are shown in Figure 5—figure supplement 4.

Figure 5—source data 1 The local net charge at the neuraminidase (NA) antigenic region of interest in different strains from 1968 to 2020.: https://cdn.elifesciences.org/articles/72516/elife-72516-fig5-data1-v2.zip
Download elife-72516-fig5-data1-v2.zip
Figure 5—source data 2 Pairwise coevolution score for pairs of charge states that emerge or disappear within 5 years.: https://cdn.elifesciences.org/articles/72516/elife-72516-fig5-data2-v2.csv
Download elife-72516-fig5-data2-v2.csv

We further tested whether epistasis is correlated with the evolution of local net charge at the NA antigenic region of interest. A coarse-grained analysis that only considered three charged states, namely, positive, negative, and neutral, for each residue was performed (see Materials and methods, Figure 5—figure supplement 1). In this analysis, amino acids were classified into (+) as positively charged (amino acids K/R), (-) as negatively charged (amino acids D/E), and (n) as neutral representing the remaining amino acids. We then computed the pairwise epistasis between different charge classes in two loci by averaging over all the corresponding amino acids in each class (Figure 5b, Figure 5—figure supplement 2). For example, the epistasis between positive charge at residue 344 and negative charge at residue 369 (i.e., +344/–369) is the averaged epistasis over K344/D369, R344/D369, K344/E369, and R344/E369. In addition, we evaluated a coevolution score between a pair of charged states at two different residues (Materials and methods, Figure 5—figure supplement 3). In our definition, two charge states that emerged or disappeared shortly one after another would have a positive coevolution score. In contrast, emergence of a charge state in one residue followed by disappearance of a charge state in another residue would result in a negative coevolution score. Our analysis only included pairs of charge states that emerged or disappeared within 5 years from each other. Since residues 367 and 370 were dominated by neutral charge since 1968, they were not included on our coevolution analysis (Figure 5—figure supplement 1).

When we compared the coevolution score with the pairwise epistasis (Figure 5c, Figure 5—figure supplement 4), high correlation was observed (Pearson correlation = 0.79–0.92). Specifically, pairs with opposite charges usually have positive coevolution scores and positive epistasis, whereas pairs with the same charge usually have negative coevolution scores and negative epistasis. We also attempted to use the fitness data to predict the identity of amino acid mutations along the evolutionary trajectory. Nonetheless, the realized evolutionary trajectories in nature are seldom the most probable ones among an ensemble of trajectories that were constructed from the fitness data (Figure 5—figure supplement 5). As a result, while our analysis demonstrates that a biophysically grounded epistatic landscape is correlated with the natural coevolution of residues in this NA antigenic region, additional considerations are needed for the prediction of the exact amino acid mutations along the evolutionary trajectory.

Discussion

Understanding the biophysical constraints is a key to model the molecular evolution of proteins (Echave and Wilke, 2017; Luksza and Lässig, 2014; Nourmohammad et al., 2013; Otwinowski et al., 2018; Mustonen et al., 2008; Wylie and Shakhnovich, 2011; Lässig et al., 2017). Through a systematic analysis of pairwise epistasis, this study shows that the local net charge is a major biophysical molecular phenotype that constrains the evolution of an antigenic region in human influenza H3N2 NA. An important feature for this antigenic region of interest is that pairwise epistasis is highly conserved across diverse genetic backgrounds, which in turn is correlated with residue coevolution in naturally circulating human influenza H3N2 virus. Although the homologous regions in the H1N1 NA and influenza B NA do not evolve as extensively as human influenza H3N2 NA (Figure 5—figure supplement 6), residue 329 coevolves with 344 in seasonal influenza H1N1 NA, indicating that similar local net charge balancing may also influence the evolution of H1N1 NA. In fact, earlier theoretical studies have indicated that charge balancing is one of the strongest signatures of correlated evolution across different protein families (Neher, 1994; Magliery and Regan, 2004; Callahan et al., 2011). Our study here further provided empirical evidence for this phenomenon.

A key finding of this study is that the optimal local net charge at the antigenic region of interest is slightly negative, and an increase or decrease in local net charge is deleterious. This characteristic suggests that the local net charge may be the key molecular phenotype under stabilizing and balancing selection, imposing a biophysical constraint on evolution of the NA antigenic region. This local net charge may have several nonexclusive functional roles. First, since the host cell membrane and sialylated glycan receptor are both negatively charged, charge distribution on the virus surface, including the antigenic region of interest, may affect the kinetics of host membrane attachment and virus release. Second, the antigenic region of interest is proximal to the catalytic active site (Figure 1a), the local net charge may affect the catalytic efficiency of NA. Lastly, the local net charge may influence the protein stability of NA (de Graff et al., 2016; Raghunathan et al., 2013). Understanding the detailed molecular mechanisms of biophysical constraints, albeit beyond the scope of this study, will likely further enhance the ability to model evolution.

One interesting observation in this study is that the ancestral strain HK68 has a much lower mutational tolerance at the antigenic region of interest than the subsequent strains, although the topology of the local fitness landscapes is largely conserved across strains. This result suggests that other biophysical features of NA are evolving over time and that they epistatically interact with the antigenic region of interest in a residue nonspecific manner. Nonspecific epistasis is often related to protein stability (Starr and Thornton, 2016), which is best described by the threshold robustness model (Bloom et al., 2006). Under the threshold robustness model, slightly destabilizing mutations may have neutral fitness when the protein has an excess stability margin but are deleterious when the protein is marginally stable. Threshold robustness model can be used to explain the differential variant fitness distribution among the six genetic backgrounds in our study. For example, it is possible that Mos99 NA has an excess stability margin such that many variants are nearly neutral despite being slightly destabilizing. In contrast, the same slightly destabilizing variants are highly deleterious in HK68 NA because it is marginally stable. When comparing the variant fitness in different genetic backgrounds, some nonlinearity can be observed (Figure 2b), which is a feature of the threshold robustness model (Otwinowski et al., 2018; Bloom Jesse et al., 2005). Nevertheless, additional studies are needed to confirm whether the difference in variant fitness distribution among genetic backgrounds is due to protein stability or other biophysical factors.

Predicting the evolution of human influenza virus is a challenging task, yet important for seasonal influenza vaccine development. An accurate predictive model of human influenza evolution likely requires an integration of epitope information (Luksza and Lässig, 2014), antigenic data (Neher et al., 2016), mutant fitness measurement (Lee et al., 2018), and experimental selection for antibody escape variants (Li et al., 2016). This work further suggests that a biophysical epistatic model of antigenic fitness landscape can also be instrumental in modeling the evolution of human influenza virus. As the knowledge about the evolutionary biology of influenza virus accumulates, a unifying model that can accurately predict emerging mutation may one day be built.

Materials and methods

Key resources table

Reagent type (species) or resource	Designation	Source or reference	Identifiers
Cell line (Homo sapiens)	MDCK-SIAT1	Sigma-Aldrich	Cat#: 05071502-1VL
Cell line (H. sapiens)	HEK 293T	Scripps Research	N/A
Strain, strain background (Escherichia coli)	MegaX DH10B T1R	Thermo Fisher Scientific	Cat#: C640003
Recombinant DNA reagent	A/WSN/33 (H1N1) eight-plasmid reverse genetics system	Neumann et al., 1999	N/A
Commercial assay or kit	Lipofectamine 2000	Thermo Fisher Scientific	Cat#: 11668019
Commercial assay or kit	KOD Hot Start DNA Polymerase	EMD Millipore	Cat#: 71086-3
Software, algorithm	Python	Python Software Foundation	RRID:SCR_008394
Software, algorithm	R	R Core Team	RRID:SCR_001905

Cell lines

Request a detailed protocol

HEK 293T (human embryonic kidney) cells and MDCK-SIAT1 (Madin–Darby canine kidney) cells were used in this study. The identification of the cell lines was confirmed morphologically. Cells were maintained in a humidified 37°C, 5% CO₂ incubator and cultured in Dulbecco’s modified Eagle’s medium (DMEM) (Life Technologies), supplemented with 10% fetal bovine serum (FBS) (VWR), and 1% penicillin-streptomycin (Life Technologies). Cells were tested monthly for mycoplasma contamination. Mycoplasma contamination was not detected.

Recombinant influenza virus

Request a detailed protocol

All H3N2 viruses generated in this study were based on the influenza A/WSN/33 (H1N1) eight-plasmid reverse genetics system (Neumann et al., 1999). Chimeric 6:2 reassortments were employed with the HA ectodomains from the H3N2 A/Hong Kong/1/1968 (HK68) and the entire NA coding region from the strains of interest (Wu et al., 2017). For HA, the ectodomain was from HK68, whereas the noncoding region, N-terminal secretion signal, C-terminal transmembrane domain, and cytoplasmic tail were from A/WSN/33. For NA, the entire coding region was from the strains of interest, whereas the noncoding region of NA was from A/WSN/33. H3N2 strains of interest in this study were as follows with GISAID (Shu and McCauley, 2017) accession numbers in parentheses: A/Hong Kong/1/1968 (EPI_ISL_245769), A/Bangkok/1/1979 (EPI_ISL_122020), A/Beijing/353/1989 (EPI_ISL_123212), A/Moscow/10/1999 (EPI_ISL_127595), A/Victoria/361/2011 (EPI_ISL_158723), and A/Hong Kong/2671/2019 (EPI_ISL_882915). For virus rescue using the eight-plasmid reverse genetics system, transfection was performed in HEK 293T/MDCK-SIAT1 cells (Sigma-Aldrich, Cat#: 05071502-1VL) that were co-cultured (ratio of 6:1) at 60% confluence using Lipofectamine 2000 (Life Technologies) according to the manufacturer’s instructions. At 24 hr post-transfection, cells were washed twice with phosphate-buffered saline (PBS) and cell culture medium was replaced with OPTI-MEM medium supplemented with 0.8 μg mL^–1 TPCK-trypsin. Virus was harvested at 72 hr post-transfection. For measuring virus titer by TCID₅₀ assay, MDCK-SIAT1 cells were washed twice with PBS prior to the addition of virus, and OPTI-MEM medium was supplemented with 0.8 μg mL^–1 TPCK-trypsin.

Mutant library construction

Request a detailed protocol

For each mutant library, insert and vector fragments were generated by PCR using PrimeSTAR Max DNA Polymerase (Takara) according to the manufacturer’s instructions, with WT NA-encoding plasmid (pHW2000-NA) as templates. Primers for the insert contained the combinatorial mutations of interest and are shown in Supplementary file 1. For insert fragment, two rounds of PCR were performed. Forward primers (P2 set) and reverse primers (P3 set) were mixed at the indicated molar ratio and used for the first round PCR (Supplementary file 1). Products from the first round PCR were then purified using Monarch DNA Gel Extraction Kit (New England Biolabs) and used as the templates for the second round insert PCR. Forward primers (P1 set) and reverse primers (P3 set) were mixed at the indicated molar ratio and used for the second round PCR (Supplementary file 1). Of note, the same reverse primers (P3 set) were used for both rounds PCR. Primers for the vector PCR are also shown in Supplementary file 1. The final PCR products of the inserts and vectors were purified by PureLink PCR purification kit (Thermo Fisher Scientific), digested by DpnI and BsmBI (New England Biolabs), and ligated using T4 DNA ligase (New England Biolabs). The ligated product was transformed into MegaX DH10B T1R cells (Thermo Fisher Scientific). At least 1 million colonies were collected for each mutant library. Plasmid mutant libraries were purified from the bacteria colonies using Plasmid Midi Kit (QIAGEN).

High-throughput fitness measurement of NA mutants

Request a detailed protocol

Each plasmid mutant library was rescued as described above (see section ‘Recombinant influenza virus’) in a T75 flask, titered by TCID₅₀ assay using MDCK-SIAT1 cells then stored at –80°C until use. To passage the virus mutant libraries, MDCK-SIAT1 cells in s T75 flask were washed twice with PBS and then infected with a multiplicity of infection (MOI) of 0.02 in OPTI-MEM medium containing 0.8 µg mL^–1 TPCK-trypsin. Infected cells were washed twice with PBS at 2 hr post-infection, then fresh OPTI-MEM medium containing 0.8 µg mL^–1 TPCK-trypsin was added to the cells. At 24 hr post-infection, supernatant containing the virus was collected. Viral RNA was extracted using QIAamp Viral RNA Mini Kit (QIAGEN). Purified viral RNA was reverse transcribed to cDNA using Superscript III reverse transcriptase (Thermo Fisher Scientific). The adapter sequence for Illumina sequencing was added to the plasmid or cDNA from the post-passaging virus mutant libraries by PCR using sequencing library preparation primers listed in Supplementary file 1. An additional PCR was performed to add the rest of the adapter sequence and index to the amplicon using primers: 5′-AAT GAT ACG GCG ACC ACC GAG ATC TAC ACT CTT TCC CTA CAC GAC GCT-3′ and 5′-CAA GCA GAA GAC GGC ATA CGA GAT XXX XXX GTG ACT GGA GTT CAG ACG TGT GCT-3′. Positions annotated by an X represent the nucleotides for the index sequence. The final PCR products were purified by PureLink PCR purification kit (Thermo Fisher Scientific) and submitted for next-generation sequencing using Illumina MiSeq PE250.

Sequencing data analysis

Request a detailed protocol

Sequencing data were obtained in FASTQ format and analyzed using custom Python scripts. Briefly, sequences were parsed by SeqIO module in BioPython (Cock et al., 2009). After trimming the primer sequences, both forward and reverse-complement of the reverse reads were translated into protein sequences. A paired-end read was then filtered and removed if the protein sequences of the forward and reverse-complement of the reverse reads did not match. Subsequently, amino acids at the residues of interest were extracted. The number of reads corresponding to each of the 864 variants was then counted. The unnormalized fitness of each variant $i$ in each replicate was estimated as follows:

f_{i} = {l o g}_{10} \frac{({O u t p u t}_{C o u n t}_{i} + 1) / ({i n p u t}_{C o u n t}_{i} + 1)}{({O u t p u t}_{C o u n t}_{W T} + 1) / ({i n p u t}_{C o u n t}_{W T} + 1)}

where the ${O u t p u t}_{C o u n t}_{i}$ represents the number of reads corresponding to variant $i$ in the post-passaging virus mutant library, and the ${i n p u t}_{C o u n t}_{i}$ represents the number of reads corresponding to variant $i$ in the plasmid mutant library. A pseudocount of 1 was added to the counts to avoid division by zero. Of note, the WT sequence of Vic11 contains a naturally rare variant T329. As a result, the WT sequence of Vic11 was not included in our mutant library design. However, due to incomplete DpnI digestion of the vector during mutant library construction, the WT sequence of Vic11 was present in the Vic11 mutant library and could be detected in the next-generation sequencing data.

The final fitness value for each mutant is

f_{i} = {l o g}_{10} (\frac{({O u t p u t}_{C o u n t}_{i, r e p 1} + 1) / ({i n p u t}_{C o u n t}_{i, r e p 1} + 1)}{({O u t p u t}_{C o u n t}_{W T, r e p 1} + 1) / ({i n p u t}_{C o u n t}_{W T, r e p 1} + 1)} + \frac{({O u t p u t}_{C o u n t}_{i, r e p 2} + 1) / ({i n p u t}_{C o u n t}_{i, r e p 2} + 1)}{({O u t p u t}_{C o u n t}_{W T, r e p 2} + 1) / ({i n p u t}_{C o u n t}_{W T, r e p 2} + 1)})

where rep1 and rep2 represent replicate 1 and replicate 2, respectively. The final fitness value of each variant is listed in Figure 2—source data 1.

Total charge of residues of interest

Request a detailed protocol

Positively charged amino acids (K/R) were assigned with a charge of 1, negatively charged amino acids (D/E) were –1, and neutral amino acids were 0. The net charge of a given variant is defined as the algebraic sum of the charges of the seven residues of interest. For example, the net charge of KNEGKKL, which has three positively charged amino acids and one negatively charged amino acid, is 2.

Modeling fitness and decomposition of interactions

Request a detailed protocol

We modeled variant fitness as a nonlinear function of the sum of the additive and the pairwise effects by MAVE-NN (Tareen et al., 2020). Statistical learning model was trained to estimate the latent phenotype $ϕ$ and prediction fitness $\hat{y}$ .

The sum of the additive and the pairwise effects was defined as a latent phenotype $ϕ_{p a i r w i s e}$ ,

ϕ_{p a i r w i s e} (\vec{x}; \vec{θ}) = θ_{0} + \sum_{l = 0}^{L - 1} \sum_{c} θ_{l : c} x_{l : c} + \sum_{l = 0}^{L - 2} \sum_{l^{`} = l + 1}^{L - 1} \sum_{c, c^{`}} θ_{l : c, l^{`} : c^{`}} x_{l : c} x_{l^{`} : c^{`}}

where L is the length of the sequences, C is the total number of amino acids, and $x_{l : c} = {\begin{cases} 1 i f c h a r a c t e r c o c c u r s a t p o s i t i o n l \\ 0 o t h e r w i s e, \end{cases}$ is the one-hot encoding of the sequence at position $l$ when amino acid is $c$ . $\vec{θ}$ represents the weight of the additive and the pairwise effects. Fitness is then modeled as a nonlinear function of the sum of tanh,

\hat{y} = g (ϕ; \vec{α}) = a + \sum_{k = 0}^{K - 1} b_{k} t a n h (c_{k} ϕ + d_{k})

where $K$ specifies the number of ‘hidden nodes’ contributing to the sum, and $\vec{α} = {a, b_{k}, c_{k}, d_{k}}$ are the trainable parameters.

To train the above model, dataset was reformatted to a set of N observations, ${({\vec{x}}_{n}, y_{n})}_{n = 1}^{N}$ , where each observation comprises sequence ${\vec{x}}_{n}$ and its fitness $y_{n}$ . Then, the dataset was randomly divided into a training set, a validation set, and a test set with a ratio of 0.64:0.16:0.2. Model was evaluated using repeated k-fold cross-validation, and hyperparameters were chosen by maximizing the R² of model prediction and Pearson correlation coefficient of model parameters (Otwinowski and Nemenman, 2013). Notably, we find that the R² of our model prediction is insensitive to regularization and largely depends on the quality of the variant fitness data (Figure 3—figure supplement 1).

Estimating epistasis between charged states

Request a detailed protocol

Amino acids are classified according to charges. Positive (+) represents positively charged amino acids (K/R), negative (-) represents negatively charged amino acids (D/E), and neutral (n) represents the remaining amino acids. All variants of the NA antigen region are then converted into charge states, named positive, negative, and neutral for each residue. For example, –344 represents K344 and R344. The epistasis value of a given charge state is the average over the epistasis values of all amino acid variants in the specified charge state.

Analysis of natural sequences

Request a detailed protocol

A total of 66,562 full-length NA protein sequences from human H3N2 were downloaded from the GISAID (http://gisaid.org) (Shu and McCauley, 2017; Supplementary file 2, Figure 1—source data 1). Amino acid sequences of NA residues 328, 329, 344, 367, 368, 369, and 370 in individual strains were extracted. Individual sequences were grouped by the year of isolation, and their mean fitness in different genetic backgrounds is plotted in Figure 2c. The human H3N2 NA protein sequences used in this study are listed in Figure 1—source data 1. The same analysis was performed on the NA homologous region of influenza A H1N1 and influenza B.

Inference of natural coevolution score

Request a detailed protocol

The change in the natural frequency of the charge states (+/-/n) at residue $i$ from year $n - 1$ to year $n$ was computed as

Δ f (s_{i}, n) = \sum_{a a_{i} \in s_{i}} f_{a a_{i}, n} - \sum_{a a_{i} \in s_{i}} f_{a a_{i}, n - 1}

where $s_{i}$ is the charge states (+/-/n) at residue $i$ , and $f_{{a a}_{i}, n}$ is the natural frequency of the amino acid variant ( $a a$ ) at residue $i$ in year $n$ . Charge state ‘+’ included amino acids K and R. Charge state ‘-’ included amino acids D and E. Charge state ‘n’ included the remaining amino acids.

All local maxima were then selected and defined as peak of frequency change using find_peaks module in SciPy. Subsequently, to evaluate the proximity of peaks among two charge states, a coevolution score was calculated as the sum of exponential weights of all pairwise peak distances from a given mutation pair $(i, j)$ :

E (i, j) = \sum_{k}^{N} C e^{- d_{k}}

where $d_{k}$ is the time separation (in years) between two peaks and C is initial value. C = 1 if both peaks have the same sign, otherwise C = –1 (i.e., one peak is positive and the other is negative). Pairs of peaks with a time separation $d_{k}$ of longer than 5 years were discarded. See Figure 5—figure supplement 3 for a schematic overview of calculating the coevolution score.

Construction of evolutionary trajectories based on the fitness data

Request a detailed protocol

To construct an ensemble of evolutionary trajectories using the fitness data, we assume that the same haplotype would not reappear along an evolutionary trajectory. In addition, only those seven-residue variants in our mutant libraries that had a net charge of –1, 0, or 1 were included. Besides, we only included those amino acid variants that were in our mutant libraries. An ensemble of evolutionary trajectories was constructed for all six strains of interest using their corresponding WT sequences as starting points. The number of steps in the constructed trajectory is equivalent to the number of mutations that naturally emerged following the isolation of the focal strain and prior to the point when the subsequent strains of interest were sampled. For example, there were three mutations that naturally merged between 1989 and 1999. As a result, three steps were included for the constructed evolutionary trajectory of Bei89. The fitness value at each step in the evolutionary trajectory was equivalent to that of the corresponding mutant, with the fitness of the WT set as 1 (see section ‘Sequencing data analysis’).

Code availability

Request a detailed protocol

Custom Python scripts for analyzing the fitness landscape data have been deposited to https://github.com/Wangyiquan95/NA_EPI (Wu, 2021; copy archived at swh:1:rev:2126f527add1c02cd9490a520d647b0700a18a2c).

Data availability

Raw sequencing data have been submitted to the NIH Short Read Archive under accession number: BioProject PRJNA742436. Custom python scripts for analyzing the deep mutational scanning data have been deposited to https://github.com/Wangyiquan95/NA_EPI, (copy archived at https://archive.softwareheritage.org/swh:1:rev:2126f527add1c02cd9490a520d647b0700a18a2c).

The following data sets were generated

1. Wang Y
2. Nc Wu
(2021) NIH Short Read Archive BioProject
ID PRJNA742436. Deep mutational scanning of human H3N2 influenza virus NA residues 328, 329, 344, 367, 368, 369, and 370.

https://www.ncbi.nlm.nih.gov/bioproject/742436

References

1. Abed Y
2. Pizzorno A
3. Bouhy X
4. Boivin G
(2011) Role of permissive neuraminidase mutations in influenza A/Brisbane/59/2007-like (H1N1) viruses
PLOS Pathogens 7:e1002431.
https://doi.org/10.1371/journal.ppat.1002431
- PubMed
- Google Scholar
1. Air GM
2. Els MC
3. Brown LE
4. Laver WG
5. Webster RG
(1985) Location of antigenic sites on the three-dimensional structure of the influenza N2 virus neuraminidase
Virology 145:237–248.
https://doi.org/10.1016/0042-6822(85)90157-6
- PubMed
- Google Scholar
1. Air GM
(2012) Influenza neuraminidase
Influenza and Other Respiratory Viruses 6:245–256.
https://doi.org/10.1111/j.1750-2659.2011.00304.x
- PubMed
- Google Scholar
(2006) Protein stability promotes evolvability
PNAS 103:5869–5874.
https://doi.org/10.1073/pnas.0510098103
- PubMed
- Google Scholar
(2010) Permissive secondary mutations enable the evolution of influenza oseltamivir resistance
Science 328:1272–1275.
https://doi.org/10.1126/science.1187816
- PubMed
- Google Scholar
(2005) Thermodynamic prediction of protein neutrality
PNAS 102:606–611.
https://doi.org/10.1073/pnas.0406744102
- PubMed
- Google Scholar
(2012) Epistasis as the primary factor in molecular evolution
Nature 490:535–538.
https://doi.org/10.1038/nature11510
- PubMed
- Google Scholar
(2011) Correlated evolution of nearby residues in Drosophilid proteins
PLOS Genetics 7:e1001315.
https://doi.org/10.1371/journal.pgen.1001315
- PubMed
- Google Scholar
1. Cock PJA
2. Antao T
3. Chang JT
4. Chapman BA
5. Cox CJ
6. Dalke A
7. Friedberg I
8. Hamelryck T
9. Kauff F
10. Wilczynski B
11. de Hoon MJL
(2009) Biopython: freely available Python tools for computational molecular biology and bioinformatics
Bioinformatics 25:1422–1423.
https://doi.org/10.1093/bioinformatics/btp163
- PubMed
- Google Scholar
(1983) Structure of the catalytic and antigenic sites in influenza virus neuraminidase
Nature 303:41–44.
https://doi.org/10.1038/303041a0
- PubMed
- Google Scholar
1. Couch RB
2. Atmar RL
3. Franco LM
4. Quarles JM
5. Wells J
6. Arden N
7. Niño D
8. Belmont JW
(2013) Antibody correlates and predictors of immunity to naturally occurring influenza in humans and the importance of antibody to the neuraminidase
The Journal of Infectious Diseases 207:974–981.
https://doi.org/10.1093/infdis/jis935
- PubMed
- Google Scholar
(2016) Highly Charged Proteins: The Achilles’ Heel of Aging Proteomes
Structure 24:329–336.
https://doi.org/10.1016/j.str.2015.11.006
- PubMed
- Google Scholar
(2015) Site-Specific Amino Acid Preferences Are Mostly Conserved in Two Closely Related Protein Homologs
Molecular Biology and Evolution 32:2944–2960.
https://doi.org/10.1093/molbev/msv167
- PubMed
- Google Scholar
1. Duan S
2. Govorkova EA
3. Bahl J
4. Zaraket H
5. Baranovich T
6. Seiler P
7. Prevost K
8. Webster RG
9. Webby RJ
(2014) Epistatic interactions between neuraminidase mutations facilitated the emergence of the oseltamivir-resistant H1N1 influenza viruses
Nature Communications 5:5029.
https://doi.org/10.1038/ncomms6029
- PubMed
- Google Scholar
1. Echave J
2. Wilke CO
(2017) Biophysical Models of Protein Evolution: Understanding the Patterns of Evolutionary Sequence Divergence
Annual Review of Biophysics 46:85–103.
https://doi.org/10.1146/annurev-biophys-070816-033819
- PubMed
- Google Scholar
1. Fowler DM
2. Fields S
(2014) Deep mutational scanning: a new style of protein science
Nature Methods 11:801–807.
https://doi.org/10.1038/nmeth.3027
- PubMed
- Google Scholar
(2013) Stability-mediated epistasis constrains the evolution of an influenza protein
eLife 2:e00631.
https://doi.org/10.7554/eLife.00631
- PubMed
- Google Scholar
1. Gulati U
2. Hwang C-C
3. Venkatramani L
4. Gulati S
5. Stray SJ
6. Lee JT
7. Laver WG
8. Bochkarev A
9. Zlotnick A
10. Air GM
(2002) Antibody epitopes on the neuraminidase of a recent H3N2 influenza virus (A/Memphis/31/98)
Journal of Virology 76:12274–12280.
https://doi.org/10.1128/jvi.76.23.12274-12280.2002
- PubMed
- Google Scholar
1. Hockenberry AJ
2. Wilke CO
(2019) Evolutionary couplings detect side-chain interactions
PeerJ 7:e7280.
https://doi.org/10.7717/peerj.7280
- PubMed
- Google Scholar
1. Hom N
2. Gentles L
3. Bloom JD
4. Lee KK
(2019) Deep Mutational Scan of the Highly Conserved Influenza A Virus M1 Matrix Protein Reveals Substantial Intrinsic Mutational Tolerance
Journal of Virology 93:e00161-19.
https://doi.org/10.1128/JVI.00161-19
- PubMed
- Google Scholar
(1990) Independent and disparate evolution in nature of influenza A virus hemagglutinin and neuraminidase glycoproteins
PNAS 87:786–790.
https://doi.org/10.1073/pnas.87.2.786
- PubMed
- Google Scholar
(2019) Epistatic interactions can moderate the antigenic effect of substitutions in haemagglutinin of influenza H3N2 virus
The Journal of General Virology 100:773–777.
https://doi.org/10.1099/jgv.0.001263
- PubMed
- Google Scholar
1. Krammer F
2. Fouchier RAM
3. Eichelberger MC
4. Webby RJ
5. Shaw-Saliba K
6. Wan H
7. Wilson PC
8. Compans RW
9. Skountzou I
10. Monto AS
(2018) NAction! How Can Neuraminidase-Based Immunity Contribute to Better Influenza Virus Vaccines?
MBio 9:e02332-17.
https://doi.org/10.1128/mBio.02332-17
- PubMed
- Google Scholar
(2011) Prevalence of epistasis in the evolution of influenza A surface proteins
PLOS Genetics 7:e1001301.
https://doi.org/10.1371/journal.pgen.1001301
- PubMed
- Google Scholar
1. Kumar S
2. Nussinov R
(2002) Close-range electrostatic interactions in proteins
Chembiochem 3:604–617.
https://doi.org/10.1002/1439-7633(20020703)3:7<604::AID-CBIC604>3.0.CO;2-X
- PubMed
- Google Scholar
(2017) Predicting evolution
Nature Ecology & Evolution 1:77.
https://doi.org/10.1038/s41559-017-0077
- PubMed
- Google Scholar
1. Lee JM
2. Huddleston J
3. Doud MB
4. Hooper KA
5. Wu NC
6. Bedford T
7. Bloom JD
(2018) Deep mutational scanning of hemagglutinin helps predict evolutionary fates of human H3N2 influenza variants
PNAS 115:E8276–E8285.
https://doi.org/10.1073/pnas.1806133115
- PubMed
- Google Scholar
1. Li C
2. Hatta M
3. Burke DF
4. Ping J
5. Zhang Y
6. Ozawa M
7. Taft AS
8. Das SC
9. Hanson AP
10. Song J
11. Imai M
12. Wilker PR
13. Watanabe T
14. Watanabe S
15. Ito M
16. Iwatsuki-Horimoto K
17. Russell CA
18. James SL
19. Skepner E
20. Maher EA
21. Neumann G
22. Klimov AI
23. Kelso A
24. McCauley J
25. Wang D
26. Shu Y
27. Odagiri T
28. Tashiro M
29. Xu X
30. Wentworth DE
31. Katz JM
32. Cox NJ
33. Smith DJ
34. Kawaoka Y
(2016) Selection of antigenically advanced variants of seasonal influenza viruses
Nature Microbiology 1:16058.
https://doi.org/10.1038/nmicrobiol.2016.58
- PubMed
- Google Scholar
1. Luksza M
2. Lässig M
(2014) A predictive fitness model for influenza
Nature 507:57–61.
https://doi.org/10.1038/nature13087
- PubMed
- Google Scholar
1. Lyons DM
2. Lauring AS
(2018) Mutation and Epistasis in Influenza Virus Evolution
Viruses 10:E407.
https://doi.org/10.3390/v10080407
- PubMed
- Google Scholar
1. Magliery TJ
2. Regan L
(2004) Beyond consensus: statistical free energies reveal hidden interactions in the design of a TPR motif
Journal of Molecular Biology 343:731–745.
https://doi.org/10.1016/j.jmb.2004.08.026
- PubMed
- Google Scholar
(1994) The structure of a complex between the NC10 antibody and influenza virus neuraminidase and comparison with the overlapping binding site of the NC41 antibody
Structure 2:733–746.
https://doi.org/10.1016/s0969-2126(00)00074-5
- PubMed
- Google Scholar
1. Memoli MJ
2. Shaw PA
3. Han A
4. Czajkowski L
5. Reed S
6. Athota R
7. Bristol T
8. Fargis S
9. Risos K
10. Powers JH
11. Davey RT
12. Taubenberger JK
(2016) Evaluation of Antihemagglutinin and Antineuraminidase Antibodies as Correlates of Protection in an Influenza A/H1N1 Virus Healthy Human Challenge Model
MBio 7:e00417.
https://doi.org/10.1128/mBio.00417-16
- PubMed
- Google Scholar
1. Miton CM
2. Tokuriki N
(2016) How mutational epistasis impairs predictability in protein evolution and design
Protein Science 25:1260–1272.
https://doi.org/10.1002/pro.2876
- PubMed
- Google Scholar
1. Monto AS
2. Petrie JG
3. Cross RT
4. Johnson E
5. Liu M
6. Zhong W
7. Levine M
8. Katz JM
9. Ohmit SE
(2015) Antibody to Influenza Virus Neuraminidase: An Independent Correlate of Protection
The Journal of Infectious Diseases 212:1191–1199.
https://doi.org/10.1093/infdis/jiv195
- PubMed
- Google Scholar
(2008) Energy-dependent fitness: a quantitative model for the evolution of yeast transcription factor binding sites
PNAS 105:12376–12381.
https://doi.org/10.1073/pnas.0805909105
- PubMed
- Google Scholar
1. Narayanan KK
2. Procko E
(2021) Deep Mutational Scanning of Viral Glycoproteins and Their Host Receptors
Frontiers in Molecular Biosciences 8:636660.
https://doi.org/10.3389/fmolb.2021.636660
- PubMed
- Google Scholar
1. Neher E
(1994) How frequent are correlated changes in families of protein sequences?
PNAS 91:98–102.
https://doi.org/10.1073/pnas.91.1.98
- PubMed
- Google Scholar
(2016) Prediction, dynamics, and visualization of antigenic phenotypes of seasonal influenza viruses
PNAS 113:E1701–E1709.
https://doi.org/10.1073/pnas.1525578113
- PubMed
- Google Scholar
1. Neumann G
2. Watanabe T
3. Ito H
4. Watanabe S
5. Goto H
6. Gao P
7. Hughes M
8. Perez DR
9. Donis R
10. Hoffmann E
11. Hobom G
12. Kawaoka Y
(1999) Generation of influenza A viruses entirely from cloned cDNAs
PNAS 96:9345–9350.
https://doi.org/10.1073/pnas.96.16.9345
- PubMed
- Google Scholar
(2013) Universality and predictability in molecular quantitative genetics
Current Opinion in Genetics & Development 23:684–693.
https://doi.org/10.1016/j.gde.2013.11.001
- PubMed
- Google Scholar
1. Otwinowski J
2. Nemenman I
(2013) Genotype to phenotype mapping and the fitness landscape of the E. coli lac promoter
PLOS ONE 8:e61570.
https://doi.org/10.1371/journal.pone.0061570
- PubMed
- Google Scholar
(2018) Inferring the shape of global epistasis
PNAS 115:E7550–E7558.
https://doi.org/10.1073/pnas.1804015115
- PubMed
- Google Scholar
(2013) Modulation of protein stability and aggregation properties by surface charge engineering
Molecular BioSystems 9:2379–2389.
https://doi.org/10.1039/c3mb70068b
- PubMed
- Google Scholar
1. Sandbulte MR
2. Westgeest KB
3. Gao J
4. Xu X
5. Klimov AI
6. Russell CA
7. Burke DF
8. Smith DJ
9. Fouchier RAM
10. Eichelberger MC
(2011) Discordant antigenic drift of neuraminidase and hemagglutinin in H1N1 and H3N2 influenza viruses
PNAS 108:20748–20753.
https://doi.org/10.1073/pnas.1113801108
- PubMed
- Google Scholar
(2019) Modeling Electrostatic Force in Protein-Protein Recognition
Frontiers in Molecular Biosciences 6:94.
https://doi.org/10.3389/fmolb.2019.00094
- PubMed
- Google Scholar
1. Shu Y
2. McCauley J
(2017) GISAID: Global initiative on sharing all influenza data - from vision to reality
Euro Surveillance 22:30494.
https://doi.org/10.2807/1560-7917.ES.2017.22.13.30494
- PubMed
- Google Scholar
1. Soh YS
2. Moncla LH
3. Eguia R
4. Bedford T
5. Bloom JD
(2019) Comprehensive mapping of adaptation of the avian influenza polymerase protein PB2 to humans
eLife 8:e45079.
https://doi.org/10.7554/eLife.45079
- PubMed
- Google Scholar
1. Starr TN
2. Thornton JW
(2016) Epistasis in protein evolution
Protein Science 25:1204–1218.
https://doi.org/10.1002/pro.2897
- PubMed
- Google Scholar
Preprint
(2020) MAVE-NN: Learning Genotype-Phenotype Maps from Multiplex Assays of Variant Effect
bioRxiv.
https://doi.org/10.1101/2020.07.14.201475
- Google Scholar
1. Thyagarajan B
2. Bloom JD
(2014) The inherent mutational tolerance and antigenic evolvability of influenza hemagglutinin
eLife 3:e03300.
https://doi.org/10.7554/eLife.03300
- PubMed
- Google Scholar
(1992) Refined crystal structure of the influenza virus N9 neuraminidase-NC41 Fab complex
Journal of Molecular Biology 227:122–148.
https://doi.org/10.1016/0022-2836(92)90687-f
- PubMed
- Google Scholar
(1983) Structure of the influenza virus glycoprotein antigen neuraminidase at 2.9 A resolution
Nature 303:35–40.
https://doi.org/10.1038/303035a0
- PubMed
- Google Scholar
(1988) Structure of an escape mutant of glycoprotein N2 neuraminidase of influenza virus A/Tokyo/3/67 at 3 A
Journal of Molecular Biology 200:201–203.
https://doi.org/10.1016/0022-2836(88)90344-0
- PubMed
- Google Scholar
(2006) An epidemiologically significant epitope of a 1998 human influenza virus neuraminidase forms a highly hydrated interface in the NA-antibody complex
Journal of Molecular Biology 356:651–663.
https://doi.org/10.1016/j.jmb.2005.11.061
- PubMed
- Google Scholar
1. Weiss CD
2. Wang W
3. Lu Y
4. Billings M
5. Eick-Cost A
6. Couzens L
7. Sanchez JL
8. Hawksworth AW
9. Seguin P
10. Myers CA
11. Forshee R
12. Eichelberger MC
13. Cooper MJ
(2020) Neutralizing and Neuraminidase Antibodies Correlate With Protection Against Influenza During a Late Season A/H3N2 Outbreak Among Unvaccinated Military Recruits
Clinical Infectious Diseases 71:3096–3102.
https://doi.org/10.1093/cid/ciz1198
- PubMed
- Google Scholar
(2012) Genetic evolution of the neuraminidase of influenza A (H3N2) viruses from 1968 to 2009 and its correspondence to haemagglutinin evolution
The Journal of General Virology 93:1996–2007.
https://doi.org/10.1099/vir.0.043059-0
- PubMed
- Google Scholar
(2015) Optimization of an enzyme-linked lectin assay suitable for rapid antigenic characterization of the neuraminidase of human influenza A(H3N2) viruses
Journal of Virological Methods 217:55–63.
https://doi.org/10.1016/j.jviromet.2015.02.014
- PubMed
- Google Scholar
1. Wu NC
2. Young AP
3. Dandekar S
4. Wijersuriya H
5. Al-Mawsawi LQ
6. Wu TT
7. Sun R
(2013) Systematic identification of H274Y compensatory mutations in influenza A virus neuraminidase by high-throughput screening
Journal of Virology 87:1193–1199.
https://doi.org/10.1128/JVI.01658-12
- PubMed
- Google Scholar
1. Wu NC
2. Du Y
3. Le S
4. Young AP
5. Zhang TH
6. Wang Y
7. Zhou J
8. Yoshizawa JM
9. Dong L
10. Li X
11. Wu TT
12. Sun R
(2016) Coupling high-throughput genetics with phylogenetic information reveals an epistatic interaction on the influenza A virus M segment
BMC Genomics 17:46.
https://doi.org/10.1186/s12864-015-2358-7
- PubMed
- Google Scholar
1. Wu NC
2. Xie J
3. Zheng T
4. Nycholat CM
5. Grande G
6. Paulson JC
7. Lerner RA
8. Wilson IA
(2017) Diversity of Functionally Permissive Sequences in the Receptor-Binding Site of Influenza Hemagglutinin
Cell Host & Microbe 21:742–753.
https://doi.org/10.1016/j.chom.2017.05.011
- PubMed
- Google Scholar
1. Wu NC
2. Thompson AJ
3. Xie J
4. Lin CW
5. Nycholat CM
6. Zhu X
7. Lerner RA
8. Paulson JC
9. Wilson IA
(2018) A complex epistatic network limits the mutational reversibility in the influenza hemagglutinin receptor-binding site
Nature Communications 9:1264.
https://doi.org/10.1038/s41467-018-03663-5
- PubMed
- Google Scholar
(2020) Major antigenic site B of human influenza H3N2 viruses has an evolving local fitness landscape
Nature Communications 11:1233.
https://doi.org/10.1038/s41467-020-15102-5
- PubMed
- Google Scholar
Software
1. Wu N
(2021) ANALYSIS FOR H3N2 NA ANTIGENIC region of interest by DEEP MUTATIONAL SCANNING, version swh:1:rev:2126f527add1c02cd9490a520d647b0700a18a2c
GitHub.

https://archive.softwareheritage.org/swh:1:dir:f47ec722f7985f90371e4f6389ca067140011aee;origin=https://github.com/Wangyiquan95/NA_EPI;visit=swh:1:snp:6f04e4bcdaa02f6bdb3708b867b3cda57657bded;anchor=swh:1:rev:2126f527add1c02cd9490a520d647b0700a18a2c
1. Wylie CS
2. Shakhnovich EI
(2011) A biophysical protein folding model accounts for most mutational fitness effects in viruses
PNAS 108:9916–9921.
https://doi.org/10.1073/pnas.1017572108
- PubMed
- Google Scholar
1. Zhu X
2. McBride R
3. Nycholat CM
4. Yu W
5. Paulson JC
6. Wilson IA
(2012) Influenza virus neuraminidases with reduced enzymatic activity that avidly bind sialic Acid receptors
Journal of Virology 86:13371–13383.
https://doi.org/10.1128/JVI.01426-12
- PubMed
- Google Scholar
1. Zhu X
2. Turner HL
3. Lang S
4. McBride R
5. Bangaru S
6. Gilchuk IM
7. Yu W
8. Paulson JC
9. Crowe JE
10. Ward AB
11. Wilson IA
(2019) Structural Basis of Protection against H7N9 Influenza Virus by Human Anti-N9 Neuraminidase Antibodies
Cell Host & Microbe 26:729–738.
https://doi.org/10.1016/j.chom.2019.10.002
- PubMed
- Google Scholar

Decision letter

Richard A Neher

Reviewing Editor; University of Basel, Switzerland
Betty Diamond

Senior Editor; The Feinstein Institute for Medical Research, United States
Richard A Neher

Reviewer; University of Basel, Switzerland
Christopher JR Illingworth

Reviewer; University of Cambridge, United Kingdom
Claus O Wilke

Reviewer; The University of Texas at Austin, United States

Our editorial process produces two outputs: i) public reviews designed to be posted alongside the preprint for the benefit of readers; ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.

Decision letter after peer review:

Thank you for submitting your article "Antigenic evolution of human influenza H3N2 neuraminidase is constrained by charge balancing" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, including Richard A Neher as Reviewing Editor and Reviewer #1, and the evaluation has been overseen by Betty Diamond as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Christopher J R Illingworth (Reviewer #2); Claus O Wilke (Reviewer #3).

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

Essential revisions:

1) The primary conclusion from your work is that evolution of NA is heavily constraint by charge balancing and that this constraint could potentially be used to predict evolution. It would be great if this could be quantified. How does the realized trajectory compare to an ensemble of possible ones? How many of the past substitutions could have been predicted? Are similar constraints seen in homologous regions of A/H1N1 of B? Some additional exploration of this would strengthen the paper.

2) The structural interpretation of the results could potentially be strengthened by analyzing side-chain positions rather than C-α positions. Comparisons of structures with different local charge might also be informative.

3) The variation among the estimates of additive effects is surprising. Is this variation smaller between similar backgrounds (10 years apart) compared to distant backgrounds? How well does an "additive only" model fit the data and how variable are the inferred effects? Some additional discussion of this would help.

4) The experiments presented here are different from most DMS experiments: All possible combinations of a subset of mutations are analyzed, rather than all possible point mutations. It would be good to highlight this difference and possibly use a different term.

The reviews below contain a number of other suggestions that we encourage you to consider.

Reviewer #1 (Recommendations for the authors):

– I have trouble interpreting the result that additive fitness effects are all over the place when comparing backgrounds, while pairwise effects are conserved. How much does this depend on the model choice? You could also use un-ambiguous Fourier transformation on subsets of binary hypercubes. And fit a model with only addititive coefficients. How much would the coefficients vary in that case? How much less variance is captured? Is there a sense that additive fitness effects are more strongly correlated with comparing backgrounds 10y appart then 20, 30, 40y apart?

– The part on analyzing natural variation in charge and prediction could be expanded and explored a little more. Net charge of the regions tends to be either -1,0, or +1, but it is not apriori clear that the realized trajectory is much more narrow in its charge distribution than random permutations of the mutations of the last 50 years are. From the experimental data, a charge of -2 seems better than +1, but +1 has been dominant in the last few years, while -2 is only sporadically observed. There is of course only one realization with a long correlation time. To strengthen the case for predictive power, would it be feasible to repeat this analysis on homologous residues in H1N1 (and maybe B and the H1N1pdm)?

Reviewer #2 (Recommendations for the authors):

1. The authors note that the region studied is subject to antigenic evolution. To what extent do the authors believe that charge balancing imposes a negative constraint upon evolution, as opposed to epistasis driving compensating changes in NA under positive selection following an immunologically driven change in the virus that coincidentally impacts residue charge? This is not to nullify the finding, but of the substiutions observed in this region over time, how many could have been predicted in advance on the basis of fitness landscapes derived here (accounting for their charge or non-charge effects)? Can the authors predict likely changes in this region that are likely to be observed next?

2. With this in mind, I found Figure 5a hard to understand or interpret. Would it be possible to show the data in a way that highlights the frequency of haplotypes by charge against time? I would find this easier to follow.

3. A small amount of work has been done looking at the consequences of charge for the local structure, though sidechain positions might be a little stochastic between crystal structures. Are there any consistent changes in structure arising from changes in local charge? That is, if multiple structures with different charge distributions are locally aligned (e.g. using VMD to perform local structural alignment on backbone atoms and RMSD to assess differences), are any changes in local conformation evident?

4. I am perhaps not sufficiently familiar with the terminology to understand what exactly was done in the virus rescue experiment. Could more details be provided?

5. Is it possible, potentially via a more superficial analysis, to suggest whether the region of charge conservation in NA might be preserved beyond the region under study? I am happy if the authors don't want to go down this route, though I wonder whether a figure along the lines of a revised Figure 5A would show possible routes for future exploration.

Reviewer #3 (Recommendations for the authors):

I think you need to be absolutely clear about whether you performed six experiments of 864 variants or 864 deep mutational scanning experiments of every possible mutation in NA. I think you did the former, but in many places the manuscript reads as if you did the latter. If you did six experiments of 864 variants, I believe you shouldn't use the term "deep mutational scanning" at all.

Neuraminidase has almost 400 residues, so a deep mutational scan of one NA variation would require almost 8000 mutations, ten times more than one of your experiment did (assuming you did only six experiments).

I also would like to emphasize that one of the main strong points of your study (if I understand correctly what you did) is that you systematically explored all possible combinations of a set of mutations. This is very different from deep mutational scanning, which usually only looks at single mutations and maybe sometimes at mutation pairs. By emphasizing deep mutational scanning, you are drawing attention away from this aspect of your study, even though that's the primary selling point of your study in my opinion.

Additionally, redoing the analysis of Figure 4 with side-chain distances should be straightforward. You could use either the smallest distance among all heavy atoms or the distance between the geometric center of the sidechain atoms, whichever is easier for you to calculate. Both should give you approximately the same results.

[Editors’ note: further revisions were suggested prior to acceptance, as described below.]

Thank you for resubmitting your work entitled "Antigenic evolution of human influenza H3N2 neuraminidase is constrained by charge balancing" for further consideration by eLife. Your revised article has been evaluated by Betty Diamond (Senior Editor) and a Reviewing Editor.

The manuscript has been improved but there are some remaining issues that need to be addressed, as outlined below:

Thank you very much for submitting the revision. I still think this study is elegant and important (as did my co-reviewers). Unfortunately, the additional analyses you have done did not provide strong support for the claim that charge conversation can be used to predict influenza evolution, while the central claim of charge being an important constraint remains convincing. I therefore think some parts of the manuscript need to be toned down. The last two sentences of abstract, for example, seem too strong:

"In addition, we show that residue coevolution in this antigenic region can be predicted by charge states and pairwise epistasis. Overall, this study demonstrates the importance of quantifying epistasis and the underlying biophysical constraints for building a predictive model of influenza evolution."

The main evidence that charge constrains natural evolution is the concordance of co-evolution scores and pairwise epistasis. I think it should be better explained what can and what can not be concluded from this and in what sense this is predictive.

The definition of the score is also somewhat confusing and I think there are some problems around lines 470 and 471. Why is s_i summed over when s_i is an argument to \delta f? The text below doesn't help to resolve the matter.

It would be nice if you made it a bit easier for the reader to piece together what exactly happened in natural H3N2 populations in the past 50y. Figure 1c shows the trajectories, but it is at times tough to link colors with amino acids and the associated changes in charge. Highlighting which events change charge would be helpful. These events underlie Figure 5a and ultimately Figure 5c. It should be possible to show more explicitly how these are connected, maybe by combining only charge changing frequency trajectories into one graph or by increasing panel 5a and annotating the curves with the underlying changes in genotype.

https://doi.org/10.7554/eLife.72516.sa1

Author response

Essential revisions:

1) The primary conclusion from your work is that evolution of NA is heavily constraint by charge balancing and that this constraint could potentially be used to predict evolution. It would be great if this could be quantified. How does the realized trajectory compare to an ensemble of possible ones? How many of the past substitutions could have been predicted? Are similar constraints seen in homologous regions of A/H1N1 of B? Some additional exploration of this would strengthen the paper.

Thank you for the comment. In the revised manuscript, we constructed an assemble of evolutionary trajectories using the fitness data, as described in a newly added subsection “Construction of evolutionary trajectories based on the fitness data” in the Methods section:

“To construct an assemble of evolutionary trajectories using the fitness data, we assume that the same haplotype would not reappear along an evolutionary trajectory. In addition, only those seven-residue variants in our mutant libraries that had a net charge of -1, 0, or 1 were included. Besides, we only included those amino acid variants that were in our mutant libraries. An assemble of evolutionary trajectories was constructed for all six strains of interest, using their corresponding WT sequences as starting points. The number of steps in the constructed trajectory is equivalent to the number of mutations that naturally emerged following the isolation of the focal strain and prior to the point when the subsequent strains of interest were sampled. For example, there were three mutations that naturally merged between 1989 and 1999. As a result, three steps were included for the constructed evolutionary trajectory of Bei89. The fitness value at each step in the evolutionary trajectory was equivalent to that of the corresponding mutant, with the fitness of the WT set as 1 (see “Sequencing data analysis” above).”

The results of this analysis are shown in Figure 5—figure supplement 6 and described in the Results section:

“We also attempted to use the fitness data to predict the identity of amino acid mutations along the evolutionary trajectory. Nonetheless, the realized evolutionary trajectories in nature are seldom the most probable ones among an ensemble of trajectories that were constructed from the fitness data (Figure 5—figure supplement 6). As a result, while our analysis demonstrates that a biophysically grounded epistatic landscape is informative for the prediction of the natural coevolution of residues in this NA antigenic region, additional considerations are needed for the prediction of the exact amino acid mutations along the evolutionary trajectory.”

Part of the results are also described in the legend of Figure 5—figure supplement 6:

“Evolutionary trajectories were constructed based on the fitness data (see Methods). The realized trajectory in nature is highlighted in red. This result demonstrates that natural evolutionary trajectory does not always maximize the replication fitness. This discrepancy is likely due to other selection pressures in nature (e.g. immune selection pressure) that were not being captured in our in vitro fitness measurement, or due to changes in the genetic background through epistatic interactions that were not captured by our experiments.”

We also explored whether the charge balancing constraint seen in H3N2 NA can also be observed in the NA homologous region of H1N1 and influenza B. Figure 5—figure supplement 7 in the revised manuscript analyzes the evolution of the NA homologous region in seasonal H1N1 virus, 2009 pandemic H1N1 Influenza A virus, and influenza B virus. The observations are discussed in the Discussion section:

“Although the homologous regions in the H1N1 NA and influenza B NA do not evolve as extensively as human influenza H3N2 NA (Figure 5—figure supplement 7), residue 329 coevolves with 344 in seasonal influenza H1N1 NA, indicating that similar local net charge balancing may also influence the evolution of H1N1 NA and influenza B NA.”

2) The structural interpretation of the results could potentially be strengthened by analyzing side-chain positions rather than C-α positions. Comparisons of structures with different local charge might also be informative.

Thank you for the suggestion. In the revised manuscript, we replaced the Cα analysis in Figure 4d-f and Figure 4—figure supplement 3 by a side-chain analysis to examine the relationship between pairwise epistasis and the side-chain distance. The result of this analysis is described in the Results section:

“Consistently, pairwise epistasis between two amino acid variants does not correlate with their side-chain – side-chain distances (Figure 4d-f and Figure 4—figure supplement 3), further substantiating that direct side-chain – side-chain interaction is not a determinant for epistasis here.”

As indicated in the legend of Figure 4—figure supplement 3:

“The side-chain – side-chain distance was computed based on the geometric center of the sidechain atoms from each residue of interest⁶⁶.”

3) The variation among the estimates of additive effects is surprising. Is this variation smaller between similar backgrounds (10 years apart) compared to distant backgrounds? How well does an "additive only" model fit the data and how variable are the inferred effects? Some additional discussion of this would help.

In the Results section of the revised manuscript, we have made a comment about the relationship between variation of additive fitness contributions and similarity of genetic backgrounds:

“Of note, the variation of additive fitness contributions across genetic backgrounds does not seem to strictly depend on the similarity between genetic backgrounds, since the correlation between the additive fitness contributions of HK68 and HK19 (Pearson correlation = 0.46) is much higher than that of HK68 and Bk79 (Pearson correlation = 0.12), which have a shorter time separation. This observation points at a possibility for complex epistatic interactions between the antigenic region of interest and other regions on NA.”

In the revised manuscript, we added an analysis by fitting an “additive only” model to the data in Figure 3—figure supplement 4.

“Consistently, an “additive only” model, without accounting for epistasis, shows a poor fit to the fitness landscape data compared to the model above with epistasis (Figure 3—figure supplement 4a), despite the fact that the inferred additive fitness effects correlate well between the two models (Figure 3—figure supplement 4b).”

4) The experiments presented here are different from most DMS experiments: All possible combinations of a subset of mutations are analyzed, rather than all possible point mutations. It would be good to highlight this difference and possibly use a different term.

In the revised manuscript, we have avoided using deep mutational scanning to describe our experimental design. We have also highlighted the difference between conventional deep mutational scanning experiments and our approach in the introduction section:

“Unlike saturation mutagenesis in conventional deep mutational scanning, combinatorial mutagenesis enables us to measure the fitness of high-order mutants.”

We also mentioned the difference in the Results section:

“Unlike conventional deep mutational scanning, which studies all possible single amino acid mutations across a protein or domain of interest, our approach analyzes all possible combinations of a subset of mutations.”

The reviews below contain a number of other suggestions that we encourage you to consider.

Reviewer #1 (Recommendations for the authors):

– I have trouble interpreting the result that additive fitness effects are all over the place when comparing backgrounds, while pairwise effects are conserved. How much does this depend on the model choice? You could also use un-ambiguous Fourier transformation on subsets of binary hypercubes. And fit a model with only addititive coefficients. How much would the coefficients vary in that case? How much less variance is captured? Is there a sense that additive fitness effects are more strongly correlated with comparing backgrounds 10y appart then 20, 30, 40y apart?

Thank you for the suggestion. The difference in the additive effect is mainly reflected in the early strains, such as HK68 and BK79, while the difference is smaller for the strains in recent years. In the Results section of the revised manuscript, we have made a comment about the relationship between variation of additive fitness contributions and similarity of genetic backgrounds (see response to Essential Revisions above). We have also added an analysis by fitting an “additive only” model to the data in Figure 3—figure supplement 4 (see response to Essential Revisions above).

– The part on analyzing natural variation in charge and prediction could be expanded and explored a little more. Net charge of the regions tends to be either -1,0, or +1, but it is not apriori clear that the realized trajectory is much more narrow in its charge distribution than random permutations of the mutations of the last 50 years are. From the experimental data, a charge of -2 seems better than +1, but +1 has been dominant in the last few years, while -2 is only sporadically observed. There is of course only one realization with a long correlation time. To strengthen the case for predictive power, would it be feasible to repeat this analysis on homologous residues in H1N1 (and maybe B and the H1N1pdm)?

As a suggested by the reviewer, we explored whether such a charge constrain is conserve in the NA homologous region of A/H1N1 and influenza B. The result is presented as Figure 5—figure supplement 7 in the revised manuscript and described in the Discussion section (see response to Essential Revisions above).

Reviewer #2 (Recommendations for the authors):

1. The authors note that the region studied is subject to antigenic evolution. To what extent do the authors believe that charge balancing imposes a negative constraint upon evolution, as opposed to epistasis driving compensating changes in NA under positive selection following an immunologically driven change in the virus that coincidentally impacts residue charge? This is not to nullify the finding, but of the substiutions observed in this region over time, how many could have been predicted in advance on the basis of fitness landscapes derived here (accounting for their charge or non-charge effects)? Can the authors predict likely changes in this region that are likely to be observed next?

In the revised manuscript, we found that the realized evolutionary trajectories in nature are seldom the most probable ones among an ensemble of trajectories that were constructed from the fitness data (see response to Essential Revisions above). As a result, while our analysis demonstrates that a biophysically grounded epistatic landscape is informative for the prediction of the natural coevolution of residues in this NA antigenic region, additional considerations are needed for the prediction of the exact amino acid mutations along the evolutionary trajectory.

2. With this in mind, I found Figure 5a hard to understand or interpret. Would it be possible to show the data in a way that highlights the frequency of haplotypes by charge against time? I would find this easier to follow.

We have updated Figure 5a to improve the clarify of data presentation.

3. A small amount of work has been done looking at the consequences of charge for the local structure, though sidechain positions might be a little stochastic between crystal structures. Are there any consistent changes in structure arising from changes in local charge? That is, if multiple structures with different charge distributions are locally aligned (e.g. using VMD to perform local structural alignment on backbone atoms and RMSD to assess differences), are any changes in local conformation evident?

Unfortunately, there are only four human H3N2 NA structures available in the PDB. They are from strains A/Memphis/31/98 (Mem98), A/Perth/16/2009 (Perth09), A/Tanzania/205/2010 (Tan10), and A/Minnesota/11/2010 (Min10). While Tan10 has a net charge of -1 in the antigenic site of interest, the other three strains all have a net charge of 0. At the antigenic site of interest, the RMSD between Tan10 and other three strains with net charge of 0 ranges from 0.18 to 0.27 Å, whereas the RMSD between those three strains with net charge of 0 ranges from 0.15 to 0.19 Å. Although these differences are neglibible, we are not comfortable in drawing conclusions regarding the consequences of charge for the local structure due to the small samples size.

4. I am perhaps not sufficiently familiar with the terminology to understand what exactly was done in the virus rescue experiment. Could more details be provided?

In the revised manuscript, details of the virus rescue experiment are described in the methods section under subsection “Recombinant influenza virus”:

“For virus rescue using the eight-plasmid reverse genetics system, transfection was performed in HEK 293T/MDCK-SIAT1 cells (Σ-Aldrich, catalog number: 050715021VL) that were co-cultured (ratio of 6:1) at 60% confluence using lipofactamine 2000 (Life Technologies) according to the manufacturer’s instructions. At 24 hours post-transfection, cells were washed twice with phosphate-buffered saline (PBS) and cell culture medium was replaced with OPTI-MEM medium supplemented with 0.8 μg ml^-1 TPCK-trypsin. Virus was harvested at 72 hours post-transfection. For measuring virus titer by TCID₅₀ assay, MDCK-SIAT1 cells were washed twice with PBS prior to the addition of virus, and OPTIMEM medium was supplemented with 0.8 μg ml^-1 TPCK-trypsin.”

5. Is it possible, potentially via a more superficial analysis, to suggest whether the region of charge conservation in NA might be preserved beyond the region under study? I am happy if the authors don't want to go down this route, though I wonder whether a figure along the lines of a revised Figure 5A would show possible routes for future exploration.

This is a great suggestion. In fact, we are wrapping up another much more in-depth study that is related to charge conservation in another region of NA. As a result, we prefer not to analyze other regions in the current study.

Reviewer #3 (Recommendations for the authors):

I think you need to be absolutely clear about whether you performed six experiments of 864 variants or 864 deep mutational scanning experiments of every possible mutation in NA. I think you did the former, but in many places the manuscript reads as if you did the latter. If you did six experiments of 864 variants, I believe you shouldn't use the term "deep mutational scanning" at all.

Neuraminidase has almost 400 residues, so a deep mutational scan of one NA variation would require almost 8000 mutations, ten times more than one of your experiment did (assuming you did only six experiments).

In the revised manuscript, we have avoided using deep mutational scanning to describe our experimental design. Instead, we described our approach as “a high-throughput experimental approach that coupled combinatorial mutagenesis and next-generation sequencing”.

I also would like to emphasize that one of the main strong points of your study (if I understand correctly what you did) is that you systematically explored all possible combinations of a set of mutations. This is very different from deep mutational scanning, which usually only looks at single mutations and maybe sometimes at mutation pairs. By emphasizing deep mutational scanning, you are drawing attention away from this aspect of your study, even though that's the primary selling point of your study in my opinion.

Thank you for the comment. In the revised manuscript, we have highlighted the difference between conventional deep mutational scanning experiments and our approach (see response to Essential Revisions above).

Additionally, redoing the analysis of Figure 4 with side-chain distances should be straightforward. You could use either the smallest distance among all heavy atoms or the distance between the geometric center of the sidechain atoms, whichever is easier for you to calculate. Both should give you approximately the same results.

Thanks for the suggestion. Geometric center of the side-chain atoms was used in our revised analysis in Figure 4.

[Editors’ note: what follows is the authors’ response to the second round of review.]

Thank you very much for submitting the revision. I still think this study is elegant and important (as did my co-reviewers). Unfortunately, the additional analyses you have done did not provide strong support for the claim that charge conversation can be used to predict influenza evolution, while the central claim of charge being an important constraint remains convincing. I therefore think some parts of the manuscript need to be toned down. The last two sentences of abstract, for example, seem too strong:

"In addition, we show that residue coevolution in this antigenic region can be predicted by charge states and pairwise epistasis. Overall, this study demonstrates the importance of quantifying epistasis and the underlying biophysical constraints for building a predictive model of influenza evolution."

In the revised manuscript, the last two sentences of abstract are toned down:

“In addition, we show that residue coevolution in this antigenic region is correlated with the pairwise epistasis between charge states. Overall, this study demonstrates the importance of quantifying epistasis and the underlying biophysical constraint for building a model of influenza evolution.”

We have also toned down other related sentences in the manuscript (see tracked changes).

The main evidence that charge constrains natural evolution is the concordance of co-evolution scores and pairwise epistasis. I think it should be better explained what can and what can not be concluded from this and in what sense this is predictive.

Edits are made in the revised manuscript to emphasize that our results show a correlation between epistasis and co-evolution, but without mentioning that our results have predictive power.

The definition of the score is also somewhat confusing and I think there are some problems around lines 470 and 471. Why is s_i summed over when s_i is an argument to \Δ f? The text below doesn't help to resolve the matter.

Thank you for catching this problem, which is now resolved in the revised manuscript with additional clarifications.

It would be nice if you made it a bit easier for the reader to piece together what exactly happened in natural H3N2 populations in the past 50y. Figure 1c shows the trajectories, but it is at times tough to link colors with amino acids and the associated changes in charge. Highlighting which events change charge would be helpful. These events underlie Figure 5a and ultimately Figure 5c. It should be possible to show more explicitly how these are connected, maybe by combining only charge changing frequency trajectories into one graph or by increasing panel 5a and annotating the curves with the underlying changes in genotype.

Thanks for the great suggestion. In the revised manuscript, labels are added to Figure 5a to indicate the key events. This is also described in the figure legend:

“Key mutation events that lead to changes in local net charge are indicated.”

https://doi.org/10.7554/eLife.72516.sa2

Article and author information

Author details

Yiquan Wang

Department of Biochemistry, University of Illinois at Urbana-Champaign, Urbana, United States

Contribution
Conceptualization, Formal analysis, Investigation, Methodology, Writing - original draft, Writing – review and editing

Contributed equally with
Ruipeng Lei

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-1954-9808
Ruipeng Lei

Department of Biochemistry, University of Illinois at Urbana-Champaign, Urbana, United States

Contribution
Conceptualization, Investigation, Methodology, Writing - original draft, Writing – review and editing

Contributed equally with
Yiquan Wang

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-4652-3400
Armita Nourmohammad
1. Department of Physics, University of Washington, Seattle, United States
2. Max Planck Institute for Dynamics and Self-Organization, Göttingen, Germany
3. Fred Hutchinson Cancer Research Center, Seattle, United States
Contribution
Formal analysis, Investigation, Methodology, Writing – review and editing

Competing interests
Reviewing editor, eLife

"This ORCID iD identifies the author of this article:" 0000-0002-6245-3553
Nicholas C Wu
1. Department of Biochemistry, University of Illinois at Urbana-Champaign, Urbana, United States
2. Center for Biophysics and Quantitative Biology, University of Illinois at Urbana-Champaign, Urbana, United States
3. Carl R. Woese Institute for Genomic Biology, University of Illinois at Urbana-Champaign, Urbana, United States
4. Carle Illinois College of Medicine, University of Illinois at Urbana-Champaign, Urbana, United States
Contribution
Conceptualization, Formal analysis, Funding acquisition, Investigation, Methodology, Supervision, Writing - original draft, Writing – review and editing

For correspondence
nicwu@illinois.edu

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-9078-6697

Funding

Deutsche Forschungsgemeinschaft (SFB1310)

Armita Nourmohammad

Max Planck Society (MPRG funding)

Armita Nourmohammad

University of Washington (Royalty Research Fund: A153352)

Armita Nourmohammad

National Institutes of Health (R00 AI139445)

Nicholas C Wu

National Institutes of Health (DP2 AT011966)

Nicholas C Wu

National Institutes of Health (R35 GM142795)

Armita Nourmohammad

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

This work was supported by DFG grant (SFB1310) for Predictability in Evolution (AN), MPRG funding through the Max Planck Society (AN), the Royalty Research Fund from the University of Washington (AN), National Institutes of Health (NIH) R35 GM142795 (AN), R00 AI139445 (NCW) and DP2 AT011966 (NCW). We thank Justin Kinney and Ammar Tareen for helpful discussion, and the Roy J Carver Biotechnology Center at the University of Illinois at Urbana-Champaign for assistance with next-generation sequencing.

Senior Editor

Betty Diamond, The Feinstein Institute for Medical Research, United States

Reviewing Editor

Richard A Neher, University of Basel, Switzerland

Reviewers

Richard A Neher, University of Basel, Switzerland
Christopher JR Illingworth, University of Cambridge, United Kingdom
Claus O Wilke, The University of Texas at Austin, United States

Publication history

Preprint posted: July 12, 2021 (view preprint)
Received: July 27, 2021
Accepted: December 7, 2021
Accepted Manuscript published: December 8, 2021 (version 1)
Version of Record published: December 17, 2021 (version 2)

Copyright

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.