Mutational fitness landscape of human influenza H3N2 neuraminidase

Abstract

Influenza neuraminidase (NA) has received increasing attention as an effective vaccine target. However, its mutational tolerance is not well characterized. Here, the fitness effects of >6,000 mutations in human H3N2 NA are probed using deep mutational scanning. Our result shows that while its antigenic regions have high mutational tolerance, there are solvent-exposed regions with low mutational tolerance. We also find that protein stability is a major determinant of NA mutational fitness. The deep mutational scanning result correlates well with mutational fitness inferred from natural sequences using a protein language model, substantiating the relevance of our findings to the natural evolution of circulating strains. Additional analysis further suggests that human H3N2 NA is far from running out of mutations despite already evolving for >50 years. Overall, this study advances our understanding of the evolutionary potential of NA and the underlying biophysical constraints, which in turn provide insights into NA-based vaccine design.

Keywords: CP: Molecular biology; deep mutational scanning; evolution; influenza; neuraminidase; protein language model; protein stability; protein structure.

Figure 1.. Deep mutational scanning of Mos99 NA head domain

The fitness of individual mutations in the Mos99 NA head domain was measured by deep mutational scanning and is shown as a heatmap. Wild-type (WT) amino acids are indicated by a black circle. Mutations in gray were excluded in our data analysis due to low input count or high occurrence in the WT sample (see STAR Methods). See also Figures S1 and S2 and Table S1.

Figure 2.. Biophysical determinants of mutational tolerance

(A) Mutational tolerance was computed by averaging the fitness of mutations at each residue. The distribution of mutational tolerance for each residue category is shown as a violin plot and a boxplot. For the boxplot, the middle horizontal line represents the median. The lower and upper hinges represent the first and third quartiles, respectively. The upper whisker extends to the highest data point within a 1.5× inter-quartile range (IQR) of the third quartile, whereas the lower whisker extends to the lowest data point within a 1.5× IQR of the first quartile. Each data point represents the mutational tolerance of one residue. (B) The differences in mutational tolerance between residue categories were computed by two-tailed Student’s t tests. The p value for each pairwise comparison is shown. (C) Mutational tolerance of each residue is shown on Mos99 NA (PDB: 7U4F). One protomer is in surface representation with individual residues colored according to their mutational tolerance. The other three protomers are in gray cartoon representation. (D) The relationship between relative solvent accessibility (RSA) of each residue on the Mos99 NA tetramer (PDB: 7U4F) and its mutational tolerance is shown. Each data point represents one residue and colored according to the residue category. The Spearman’s rank correlation coefficient (ρ) is indicated. See also Figure S3 and Table S2.

Figure 3.. A cluster of residues with low mutational tolerance at the membrane proximal side of the NA head domain

(A) Mutational tolerance of each residue is shown on Mos99 NA (PDB: 7U4F). One protomer is in surface representation with individual residues colored according to their mutational tolerance. The other three protomers are in gray surface representation. The C-terminal of the stalk domain, which connects the NA head domain to the membrane, is indicated. A cluster of residues with low mutational tolerance is circled. (B) Residues within the cluster in (A) are shown. Electrostatic interactions are represented by dashed lines. (C) The natural variants of NA residues 283, 288, 304, 355, and 383 (N2 numbering) are represented by sequence logos. Top: 58,937 human H3N2 strains from 1968 to 2020. Middle: 6,900 avian N2 strains. Bottom: representative strains from N1-N9 subtypes. See also Table S2.

Figure 4.. Deep mutational scanning result correlates with mutational stability and natural evolution

(A) The relationship between predicted stability effect (ΔΔG) by FoldX^, and fitness measurement by deep mutational scanning is shown. Mutations with predicted ΔΔG >30 kcal mol⁻¹ are shown as 30 kcal mol⁻¹. A smooth curve was fitted by loess. (B) The relationship between fitness inferred from natural sequences using MSA transformer and fitness measured by deep mutational scanning is shown. The Spearman’s rank correlation coefficient (ρ) is indicated. (C and D) Same as (B) except only those mutations that are (C) observed in nature or (D) not observed in nature are shown. (E) Distributions of fitness effects of mutations that have been observed in nature (observed) and those that have not been observed in nature (unobserved) are shown. See also Figure S4 and Tables S2 and S3.