Genetic and antigenic divergence in the influenza A(H3N2) virus circulating between 2016 and 2017 in Thailand

Influenza virus evolves rapidly due to the accumulated genetic variations on the viral sequence. Unlike in North America and Europe, influenza season in the tropical Southeast Asia spans both the rainy and cool seasons. Thus, influenza epidemiology and viral evolution sometimes differ from other regions, which affect the ever-changing efficacy of the vaccine. To monitor the current circulating influenza viruses in this region, we determined the predominant influenza virus strains circulating in Thailand between January 2016 and June 2017 by screening 7,228 samples from patients with influenza-like illness. During this time, influenza A(H3N2) virus was the predominant influenza virus detected. We then phylogenetically compared the hemagglutinin (HA) gene from a subset of these A(H3N2) strains (n = 62) to the reference sequences and evaluated amino acid changes in the dominant antigenic epitopes on the HA protein structure. The divergence of the circulating A(H3N2) from the A/Hong Kong/4801/2014 vaccine strain formed five genetic groups (designated I to V) within the 3C.2a clade. Our results suggest a marked drift of the current circulating A(H3N2) strains in Thailand, which collectively contributed to the declining predicted vaccine effectiveness (VE) from 74% in 2016 down to 48% in 2017.


Introduction
Influenza A virus is an important respiratory pathogen responsible for the annual influenza outbreak and considerable socio-economic burden on the public healthcare system [1]. The multivalent influenza virus vaccine administered annually can reduce the risk of morbidity and mortality, but it is dependent on how well the chosen strains included in the vaccine match the strains in circulation [2]. The commonly circulating seasonal influenza A subtypes are A(H1N1)pdm09 and A(H3N2), of which the latter is reportedly associated with a high rate of hospitalization and mortality in the United States in the 2016-2017 flu season [3].
The hemagglutinin (HA) surface glycoprotein of A(H3N2) possesses defined antigenic and receptor-binding sites [4,5]. The HA diversity resulting from accumulated mutations facilitates a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 viral escape from the host immune response [6,7]. The HA protein is proteolytically processed into two subunits (HA1 and HA2) [8]. The globular HA1 domain contains five antigenic sites (A through E) [9][10][11], while the HA2 stem domain mediates fusion and viral uncoating [12,13]. Sequence drifts on the HA from accumulated mutations are observed more frequently in the A(H3N2) than A(H1N1), which often lead to the gradual reduction of the vaccine effectiveness (VE) over time [14][15][16]. As a result, influenza virus strains most suitable for vaccine production are carefully evaluated each year [17,18].
Timely analysis of the genetic variations on the HA1 sequence of the circulating influenza virus strains is crucial for the prediction of VE. We therefore determined whether A (H3N2) was regionally predominant in the current influenza season and compared the genetic composition of the circulating A(H3N2) to the current vaccine strain A/Hong Kong/4801/2014.

Ethical approval
Respiratory samples from patients with influenza-like illness were analyzed in the Center of Excellence in Clinical Virology at King Chulalongkorn Memorial Hospital as part of the routine influenza surveillance program. This study was approved by the Institutional Review Board (IRB) of the Faculty of Medicine at Chulalongkorn University (IRB No. 377/ 57). The IRB waived the need for consent because the samples were de-identified and anonymous.

Samples and HA gene amplification
A total of 7,228 samples obtained between January 2016 and June 2017 in Bangkok and Khon Kaen province were collected from patients with fever >38˚C and respiratory symptoms such as sore throat, nasal congestion, cough and runny nose. These individuals sought medical care at King Chulalongkorn Memorial Hospital, Bangpakok 9 International Hospital, and Chum Phae Hospital. Samples were subjected to viral RNA extraction (GeneAll Biotechnology, Seoul, Korea) according to the manufacturer's instructions. We used a previously described real-time reverse-transcription polymerase chain reaction (RT-PCR) to identify influenza virus A(H1N1pdm09), A(H3N2), and influenza B virus [19]. Influenza B virus-positive samples were subjected to cDNA synthesis using ImProm-II reverse transcription system (Promega, Madison, WI) and primer FluB (5'-AGCAGAAGCA-3') [20], followed by lineage determination using multiplex PCR and melting curve analysis [21,22]. Among A(H3N2)-positive samples, 62 strains (approximately 4 strains per month) were randomly selected for cDNA synthesis using primer Uni12 (5'-AGCAAAAGCAGG-3') [23] and the entire HA gene amplified using published primer sets [15]. Briefly, the reaction mixture consisted of 5 μl PRIME MasterMix (5Prime, Hamburg, Germany), 0.25 mM of MgCl 2 , 0.5 μM each of forward and reverse primers, 2 μl of cDNA template, and nuclease-free water to a final volume of 25 μl. Amplification in a thermal cycler was performed under the following conditions: initial denaturation for 3 minutes at 94˚C, followed by 40 cycles of 30 seconds at 94˚C for denaturation, primer annealing for 30 seconds at 55˚C, 90 seconds at 72˚C for extension, and 7 minutes of final extension at 72˚C. Amplicons were agarose gel-purified using Expin Combo GP kit (GeneAll Biotechnology, Seoul, Korea) and the HA gene sequenced. A(H3N2) nucleotide sequences were assembled using the SeqMan Pro software (DNASTAR, Madison, WI) and deposited in the GenBank database under the accession numbers (MF673231-MF673292) (S1 Table).

Phylogenetic analysis
A total of 91 A(H3N2) HA nucleotide sequences, 62 obtained from this study and an additional 29 sequences identified in Thailand publicly available from the NCBI (http://www.ncbi. nlm.nih.gov) and GISAID (http://platform.gisaid.org) databases, were aligned using ClustalW and translated into amino acid residues using BioEdit Software version 7.0.9.1 (S2 Table). Phylogenetic analysis was performed using the maximum likelihood method and HKY+G model implemented in MEGA 6 [24]. Bootstrapping was done in 1,000 replicates and values >70% were shown. Potential N-linked glycosylation sites on the HA was determined using the NetN-Glyc 1.0 server with the threshold value of >0.5 for the mean potential score [25]. The selective pressure or the proportion between non-synonymous to synonymous substitutions (dN/dS, defined as ω) observed on all HA sequences were analyzed using the single likelihood ancestor counting (SLAC), fixed effects likelihood (FEL), and mixed effects model of evolution (MEME) algorithms implemented in the HYPHY software [26]. Positively selected codon was considered significant at P-value of 0.

Estimates of vaccine effectiveness
The predicted vaccine effectiveness (VE) was estimated using the P epitope model, which characterizes the antigenic distance between the A(H3N2) vaccine and circulating strains. Antigenic distance defined by the P epitope was calculated from the fraction of substituted amino acid residues in the dominant HA epitope [28]. For A(H3N2), the association between the VE and the P epitope is given by VE = -2.47 × P epitope + 0.47 in which VE is 47% when P epitope = 0.

Statistical analysis
Statistical analyses were performed using the Statistical Package for Social Sciences version 22.0 (SPSS Inc., Chicago, USA). The one-way ANOVA test was used to analyze VE divergence, and p < 0.05 was considered statistically significant.

Results
In all, 15 Residue variations on the HA affecting antigenic epitopes, receptorbinding site, and potential glycosylation strains were defined by residue pattern S47T+T131K+R142K+S219F+R261Q. Meanwhile, group IV strains differed from the vaccine strain by N121K+S144K. Most interestingly, the novel group V formed by strains identified in this study possessed a unique constellation of residue variations N31S+D53N+R142G+S144R+N171K+I192T+Q197H (Fig 3). A(H3N2) strains circulating in Thailand collectively showed a total of 39 residue substitutions in the antigenic epitopes. Some of the strains possessed as many as seven amino acid changes (S1  [15], substitutions at N122D (1/91) and S124N (2/91) found in some strains in this study eliminated the potential glycosylation site at residue 122 (S1 Fig). Additionally, there was a loss of potential glycosylation site at residue 126 due to N126S substitution (1/91), at residue 133 due to T135K (9/91), and at residue 158 due to N158K (3/91) and T160K (2/91) substitutions.
The receptor-binding sites of A(H3N2) are highly conserved at positions 98, 136,153,183,190,194,195 and 228 on HA1 [13]. We observed T135K residue change adjacent to the receptor-binding region for all 3C.2a group II strains. Evolutionary selective pressure on the entire HA amino acid sequence analyzed using the dN/dS ratio showed an average of 0.305 (SLAC algorithm) among the A(H3N2) strains in this study, suggesting no positively selected site. Meanwhile, the FEL method revealed that residues 142, 171, 261 and 406 (HA numbering) were under positive selection pressure. Moreover, the MEME method indicated five positively selected sites (131, 144, 171, 261 and 416). Therefore, these data provided strong evidence of positively selected sites within epitope A (131, 142 and 144), D (171) and E (261) (S3 Table).

Predicted vaccine effectiveness
The average P epitope value was 0.049 (n = 53) for all of the A(H3N2) strains identified in 2016, suggesting a predicted VE against the virus of 74.17%. Meanwhile, the average P epitope was 0.099 (n = 38) for strains identified in 2017, suggesting a predicted VE against the virus of 47.87% (S4 Table). Although the predicted VE between the strains identified in this study and the A/Hong Kong/4801/2014 vaccine strain ranged between 17.02 and 87.18%, there was an overall decline in predicted VE each quarter (Fig 4). These results suggest that the A(H3N2) strains circulating in Thailand drifted from the vaccine strain and effectively reduced the VE.

Discussion
In this study, we aimed to characterize the circulating A(H3N2) in Thailand beginning in 2016 by determining clade designation and identifying mutations in the antigenic sites impacting the predicted VE. We categorized the strains identified in this study and those reported elsewhere between 2016 and 2017 using the HA sequences. Due to evolving genetic variations of A Dominant epitope (A-E) determined for each sequence are marked with different colored dots (A = yellow, B = pink, C = purple, D = green, and E = blue). The signature amino acid substitutions (colored) occurring on the antigenic epitopes are also shown. Scale bar represents approximately 0.5% nucleotide difference between close relatives. Residue numbers are specific for HA1 (color-coded by epitope) and HA2 (dark blue). The vaccine strain A/Hong Kong/4801/2014 belonged to 3C.2a clade. Shaded area (Group V) highlights strains of interest. https://doi.org/10.1371/journal.pone.0189511.g002 Influenza A(H3N2) virus circulation in Thailand between 2016-2017 (H3N2) and the lagging nomenclature standard, we attempted to reconcile the most recent 3C.2a strains identified to date. Group I and II strains in this study collectively represented clade 3C.2a1 designated by the WHO [32]. Group III strains in this study comprised some newer Israel strains and the proposed 3C.2a1 cluster III from the U.K., while group IV described additional Israel strains and the proposed 3C.2a2 clade [30,31]. Additionally, newly emergent A(H3N2) strains identified in this study necessitated further subdivision into a separate group V according to a number of significant variance from the A/Hong Kong/4801/2014 vaccine strain and other previously described variants. Members of group V were characterized by seven amino acid substitutions at N31S, D53N, R142G, S144R, N171K, I192T and Q197H. These residue variants have a number of implications including alterations of the antigenic epitopes and immune escape. For example, residue 144 in antigenic site A is adjacent to residue S145 implicated in receptor-binding [33].
New antigenic variants emerge when at least one substitution occurred in the antigenic sites [11,34]. Several strains circulating in Thailand (5/91) belonging to group II and IV demonstrated F193S change located in epitope B, which is one of seven mutation sites representing a major antigenic transition cluster [5]. Both T135K and N126S in some A(H3N2) strains resulted in the loss of N-linked glycosylation, an important observation to consider since the gain or loss of N-linked glycosylation can affect influenza virus virulence and recognition by neutralizing antibody [35]. T135K located adjacent to the receptor-binding site in epitope A has also been identified in other studies and is implicated in decreasing vaccine effectiveness [36,37]. Taken together, these variations underscore the rapid evolution of A(H3N2) influenza virus in circulation.
Southeast Asian countries use the influenza vaccine formulated for the southern hemisphere, which for some years have different inclusion strains in the vaccine than those formulated for the northern hemisphere. A/Hong Kong/4801/2014 belonging to clade 3C.2a is a component in both northern and southern hemispheres since 2016. This vaccine strain was well-matched for the circulating A(H3N2) in Thailand that first season (predicted VE of approximately 80%) [38]. Since then, A(H3N2) strains circulating in Europe and Canada have genetically drifted away from the vaccine strain consistent with the observed antigenic drift and decreasing predicted VE we found in this study for each quarterly period beginning in 2016 [39][40][41]. As a result, A/Hong Kong/4801/2014 may not be as effective in eliciting immunity against future circulating A(H3N2) in the next influenza season despite the decision to include it in the 2017-2018 vaccine for northern hemisphere [42].
This study had several limitations. We were unable to ascertain the vaccination status of the individuals in which the samples were derived, which would have improved the evaluation of the predicted VE. The scope of this study did not allow us to determine additional evolutionary relationships among the strains identified in different countries, which would have required additional nucleotide sequences from around the world not yet deposited in the databases and a longer study time frame. The antigenic drift and predicted VE were estimated from the accumulated mutations on the antigenic epitopes and would benefit from additional antigenic characterization such as hemagglutination inhibition assay. Finally, immunity against the neuraminidase contributing to the antigenic drift was not assessed in this study. Nevertheless, any genetic surveillance of influenza viruses will continue to be an important component in influenza prevention and vaccine improvement.