Genetic clades determination of the influenza virus in Japan has been performed routinely as the part of the work of the National Epidemiological Surveillance of Infectious Diseases in Japan¹ and the Global Surveillance of Influenza in the WHO Reference Laboratories². Information on the specimen collection and weekly report of HA type is available on the National Epidemiological Surveillance of Infectious Diseases³. Briefly, RNAs were extracted from viruses by using QIAamp viral RNA kit (QIAGEN, Dusseldorf, Germany). The HA genes were amplified from extracted RNAs by RT-PCRs using gene-specific primers (the primer sequences are available upon request) and SuperScript III One-step RT-PCR system with Platinum Taq (Invitrogen, Carlsbad, CA, USA). Sequencing reactions were performed with BigDye terminator kit (Applied Biosystems) and sequences were determined using 3730xl DNA analyzer (Applied Biosystems, Foster City, CA, USA). The genetic clades were determined by the generation of phylogenetic trees of the HA genes. The phylogenetic trees were constructed using MEGA 6 software (Tamura et al., 2013) with the neighbor-joining method. The numbers of HA sequences obtained in the individual periods are 39, 83, 125, 108, 70, and 77 for the September 2013 – Jan. 2014, February 2014 – August 2014, September 2014 – January 2015, February 2015 – August 2015, September 2015 – January 2016, and February 2016 – August 2016, respectively. The nucleotide sequences used in this study are registered at GISAID, a publicly accessible influenza virus database⁴. Accession number at GISAID for the HA sequence used for the molecular modeling in this study is EPI543763 (A/Switzerland/9715293/2013 strain).

Three-dimensional (3-D) models for glycosylated extracellular domains of HA trimers of influenza A (H3N2) in the ligand-free state were constructed by homology modeling with Molecular Operating Environment (MOE) (Chemical Computing Group Inc., Montreal, QC, Canada). The crystal structure of the HA trimer of the influenza A/Victoria/361/2011 (H3N2) virus (PDB code: 4O5N; resolution: 1.75 Å; amino acid residues 4–325 and 330–502 for HA1 and HA2 peptides, respectively) was used as the modeling template. Obtained models were optimized by energy minimization using MOE and an Amber10: Extended Huckel Theory (EHT) force field implemented in MOE, which combines Amber10 and EHT bonded parameters for the large-scale energy minimization (Gerber and Muller, 1995; Case et al., 2005). The high-mannose oligosaccharide Man₅GlcNAc₂ was added to potential N-glycosylation sites in HA using Online Glycoprotein Builder⁵.

Glycosylated HA trimer models in a ligand-free state were subjected to MD simulation essentially as described for simulations of HIV-1 gp120 (Yokoyama et al., 2016). MD simulations were performed by the PMEMD (Particle Mesh Ewald Molecular Dynamics) module in the AMBER 14 program package (Case et al., 2014), employing the Amber ff99SB-ILDN force field, a protein force field with improved side-chain torsion potentials (Lindorff-Larsen et al., 2010), the GLYCAM06 force field, a biomolecular force field for glycans (Kirschner et al., 2008), and the TIP3P water model for simulations of aqueous solutions (Jorgensen et al., 1983). Bond lengths involving hydrogen were constrained with SHAKE, a constraint algorithm to satisfy a Newtonian motion (Ryckaert et al., 1977), and the time step for all MD simulations was set to 2 fs. A non-bonded cutoff of 10 Å was used. After heating calculations for 20 ps until 310 K using the NVT ensemble for the constant volume, temperature, and numbers of particles in the system, simulations were executed using the NPT ensemble for the constant pressure, temperature, and numbers of particles in the system at 1 atm, at 310 K, and in 150 mM NaCl for 100 ns. Root mean square deviations (RMSDs) between the heavy atoms of the two superposed proteins were used to measure the overall structural differences between the two proteins (Case et al., 2005). The RMSD was calculated using the cpptraj module in AmberTools 14, a trajectory analysis tool (Case et al., 2014). We used “Computer System for the Prediction of Mutations of Pathogens” at Research Center for Zoonosis Control, Hokkaido University for the MD simulations.

We calculated RMSFs of individual components of the high mannose oligosaccharides Man₅GlcNAc₂ around the receptor binding site between 50 to 100 ns of MD simulations to quantify structural dynamics of glycans during the MD simulations. The average structures during the last 50 ns of MD simulations were used as reference structures for RMSF calculation. RMSFs were calculated as previously described (Naganawa et al., 2008; Yokoyama et al., 2012, 2016; Kuwata et al., 2013) by using the ptraj module in Amber, a trajectory analysis tool (Case et al., 2005).

Over the 2014/15 influenza season, influenza surveillance reports from reference laboratories warned of a global epidemic of a new A(H3N2) variant population, termed 3C.2a. In Japan, the 3C.2a was first detected as a relatively minor population over the season from February to August in 2014, during which four A(H3N2) clades co-existed (Figure 1, upper panel). However, the 3C.2a rapidly became dominant at the beginning of the winter season in 2015, displaced pre-existing clades, and has continually predominated at the collection sites in Japan ever since, representing 76 of 77 (98.7%) clades in the season from February to August in 2016 (Figure 1, lower panel). In parallel, another A (H3N2) clade, 3C.3a, which co-existed in the same period as the 3C. 2a over the February to August 2014 season, became minor after February 2015 (Figure 1). These results are consistent with the reports on the global epidemic of the 3C.2a from the WHO Reference Laboratories⁶, and suggest the selective advantage of the 3C.2a for human-to-human transmission during the study period. Notably, hemagglutination activity of the 3C.2a was significantly attenuated when measured with a conventional hemagglutination assay (Skowronski et al., 2016). Together, these data suggest that certain structural changes occurred in the HA protein to confer selective advantages on the 3C.2a.

Molecular epidemiological data from above study and the global surveillance of influenza in the WHO Reference Laboratories suggest that the A (H3N2) clades, 3C.2a and 3C.3a, were diversified from 3C.2 and 3C.3, respectively and thereafter the 3C.2a had dominated over the 3C.3a. Therefore, we examined differences in HA proteins between the clades 3C.2a and 3C.3a. The HA protein of the 3C.2a population has multiple amino acid substitutions as compared with the clade 3C.3a (amino acid numbers 3, 128, 138, 142, 144, 159, 160, 311, 326, and 489). In this study, we focused on the four substitutions around the receptor-binding surface on the tip of the HA protein, i.e., Ala128Thr, Asn144Ser, Ser159Tyr, and Lys160Thr (Figure 2). The Ala128Thr and Lys160Thr create new potential N-glycosylation sites, Asn-X-Ser/Thr, whereas Asn144Ser results in the loss of a single N-glycosylation site. Ser159Tyr is located adjacent to the Lys160Thr. The Ala128Thr substitution initially detected in the 3C.2 and had been preserved in its descendent, 3C.2a, whereas the other three substitutions newly emerged in the 3C.2a (Figure 2A). All the mutations are placed at or near the major antigenic elements for H3 HAs, such as antigen sites A (amino acids 140-146) and B (amino acids 155-160 and 188-198) (Webster and Laver, 1980; Gerhard et al., 1981; Wiley et al., 1981; Knossow et al., 1984; Wiley and Skehel, 1987) (Figure 2B).

To understand the structural impacts of the four mutations on the tip of the 3C.2a HA protein, we constructed molecular models of the glycosylated HA trimers with and without the mutations. Two structural models of the glycosylated HA trimers in a ligand-free state were constructed by homology modeling: a model for the A(H3N2) 2015/16 vaccine strain (A/Switzerland/9715293/2013) that belongs to the perishing clade 3C.3a (Figures 1, 2) and a model for A/Switzerland/9715293/2013 possessing the four 3C.2a-specific mutations on the HA globular heads (Ala128Thr, Asn144Ser, Ser159Tyr, and Lys160Thr). The obtained models were optimized by energy minimization using the MOE:EHT force field (Gerber and Muller, 1995; Case et al., 2005) and were subjected to the MD simulations using the Amber ff99SB-ILDN force field (Lindorff-Larsen et al., 2010) and the GLYCAM06 force field (Kirschner et al., 2008) as described previously (Yokoyama et al., 2016).

The structural dynamics of the glycosylated HA proteins in solution were monitored by RMSD between the initial model structure and the structures at given time points of the MD simulation (Figure 3). For the protein portions of the glycosylated HA molecules, the RMSDs sharply increased in the beginning and reached a near plateau after 20 ns of the MD simulations in both the 3C.3a and its mutant with 3C.2a mutations (Figures 3A,B). This RMSD profile was similar to that obtained with the glycosylated HIV-1 gp120 protein (Yokoyama et al., 2016). In addition, we monitored the RMSDs of the glycan portions of the glycosylated HA molecules. We obtained an RMSD profile that was similar to that of the protein portions: there was a sharp increase soon after the MD simulation onset followed by a near plateau after 20 ns (Figures 3C,D). These results suggest that the structural distortions of the amino acid residues and glycans of the initial models were relieved shortly after the start of MD simulation under thermodynamic driving forces in solution. The data also predict that the HA structure can reach a state of thermodynamic equilibrium in solution.

Finally, we compared the 3-D structure of the HA trimer in a state of thermodynamic equilibrium during the MD simulations (100 ns) between the 3C.3a and its mutant (Figure 4A). A marked structural difference between the 3C.3a and its mutants was detected on the globular heads of the HA protein at 100 ns of two MD simulations: the glycan moieties were placed more compactly in the mutant than in the 3C.3a (Figures 4A,B). Consequently, the space for the access of ligands around the receptor-binding site was more confined in the mutant (Figure 4C). The apical spaces for ligand binding could be influenced not only by the arrangements of glycans but also by fluctuation of the glycans. Therefore, we examined structural fluctuations of the individual glycans around the receptor binding site by calculating RMSF using snapshots of structures from 50 to 100 ns of each MD simulation. In the present HA models, the high mannose oligosaccharides Man₅GlcNAc₂ were attached on the asparagine residues of potential N-glycosylation sites around the receptor binding site (Figure 5A). We calculated RMSFs of individual components of the oligosaccharides, such as GlcNAc and Man, at positions 1 to 7 (Figure 5B). Notably, RMSFs of the glycan components at given positions, and the profiles of RMSFs at the seven positions were similar among the four N-glycosylation sites. The data indicate that all glycans fluctuated in a similar fashion around the receptor binding site during the MD simulations. These results suggest that the apical spaces for ligand binding are mainly influenced by arrangements of glycans and that the spaces are constantly different between the two models during the MD simulations.