Characteristics and evolution of hemagglutinin and neuraminidase genes of Influenza A(H3N2) viruses in Thailand during 2015 to 2018

Ethics statement

The study was conducted under the Declaration of Helsinki and approved by the Institutional Review Board of the Faculty of Medicine, Chulalongkorn University, Bangkok, Thailand (IRB No. 558/62), and the date of approval is September 26, 2019. The Institutional Biosafety Committee (IBC) of the Faculty of Medicine, Chulalongkorn University, approved this study with MDCU-IBC007/2019.

Clinical specimens

In this study, leftover respiratory samples were obtained from the Virology laboratory, Microbiology Department, King Chulalongkorn Memorial Hospital, Thai Red Cross, Bangkok, Thailand. In accordance with IRB regulations, the hospital director’s permission was given without the patient’s consent, which is a prerequisite for IRB approval. A total of 443 leftover respiratory samples that tested positive for Influenza A by QuickNavi Flu+RSV (Denka Seiken, Japan) from 2015 to 2018 and stored at −80 °C after testing were enrolled in this study.

Determination of influenza A virus subtype by multiplex real-time RT-PCR

The protocol for subtyping the influenza A virus into A(H1N1)pdm09 and A(H3N2) including the primers and probes was developed and validated in our laboratory. Viral RNA was extracted from 140 ml of the patient’s respiratory samples using the QIAamp Viral RNA Mini Kit (Qiagen, Hilden, Germany) following the manufacturer’s instructions. For the synthesis of complementary DNA (cDNA), 5 ml of extracted viral RNA was converted into cDNA using the Transcriptor First-Strand cDNA synthesis kit (Roche, Basel, Switzerland) and Uni12 primer (5′-AGCRAAAGCAGG-3′). Subsequently, multiplex real-time PCR was performed using the Luna Universal Probe qPCR Master Mix (New England Biolabs, Ipswich, MA, USA). The sequences of primers and hydrolysis probes were as follows: A(H1N1)pdm09-F (5′ATAYTRAGAACWCAAGAGTCTGAATG 3′), A(H1N1)pdm09-R (5′CCATGCCARTTRTCYCTG 3′), A(H1N1)pdm09-P (5′FAM-CACTAGAATCAGGRTAACAGGAGCA-BHQ1 3′), A(H3N2)-F (5′ ACAGGATTTGCACCTTTYTC 3′), A(H3N2)-R (5′ GGAACACCYAAYTCATTCATC 3′), A(H3N2)-P (5′ HEX-ACCTTATGTGTCATGCGA-BHQ1 3′). Briefly, 5 µl of viral cDNA was added to 20 µl of total reaction volume containing 1x Luna Universal Probe qPCR Master Mix, 0.2 µM of A(H3N2) primers and probe, 0.4 µM of A(H1N1)pdm09 primers, and 0.1 µM of A(H1N1)pdm09 probe. The PCR reaction was done on the QuantStudio 5 Real-Time PCR System (Applied Biosystems, USA). The cycle steps of amplification were as follows: initial denaturation at 95 °C for 60 s, followed by 45 cycles of denaturation at 95 °C for 15 s, and extension at 60 °C for 30 s. Influenza A(H1N1)pdm09 was detected in the fluorescein (FAM; 530 nm) channel, while A(H3N2) was detected in the hexachlorofluorescein (HEX; 560 nm).

Sequencing of Hemagglutinin and Neuraminidase genes

Viral RNA was extracted and converted into cDNA as described in the previous step. The amplification and sequencing primers were followed in the previous study (Lee et al., 2013). The Platinum Taq DNA Polymerase High Fidelity (Invitrogen, Waltham, MA, USA) was used to amplify the HA and NA genes following the manufacturer’s instructions. The cycle steps of PCR amplification were adapted as follows: 2 cycles of denaturation at 95 °C for 5 min, annealing at 47 °C for 30 s, and extension at 68 °C for 1.5 min; 48 cycles of 95 °C for 15 s, at 58 °C (for NA gene) or 61 °C (for HA gene) for 30 s, and at 68 °C for 1.5 min, followed by the final extension at 68 °C for 10 min. All PCR steps were performed using the ProFlex PCR System (Applied Biosystems, Foster City, CA, USA). The amplified product was checked with 1.0% agarose gel electrophoresis (Affymetrix, Santa Clara, CA, USA) and purified using ExoSAP-IT PCR Product Cleanup Reagent (Applied Biosystems) according to the manufacturer’s recommendations. The sequences of the HA and NA genes of influenza A(H3N2) viruses were obtained by bidirectional Sanger sequencing (Bio Basic, Markham, Ontario, Canada). The BioEdit program version 7.2, was used to analyze the data. The nucleotide sequences of the viral HA and NA genes obtained from sequencing were submitted to the GenBank database.

Phylogenetic analysis

The phylogenetic trees of full-length HA and NA genes of influenza A(H3N2) viruses found in this study and the reference strains obtained from the GenBank and Global Initiative on Sharing Avian Influenza Data (GISAID) databases were constructed using the MEGA X program (Kumar et al., 2018). The Maximum Likelihood method was employed to construct the phylogenetic trees using 1,000 bootstrap replications. The Tamura 3-parameter and Gamma distributed (T92 + G) substitution model was used for the NA gene, while the Hasegawa-Kishino-Yano and Gamma distributed (HKY + G) substitution model was used for the HA gene.

The analysis of similarity/identity between the circulating viruses and the vaccine strains

The Matrix Global Alignment Tool (MatGAT v2.01) application (Campanella, Bitincka & Smalley, 2003) to determine the identity or similarity between the HA and NA genes of circulating A(H3N2) viruses and the vaccine strains in the relevant year. The BLOcks Substitution Matrix (BLOSUM) 62 was selected for the scoring matrix and aligned with the sequences.

Homology modeling of protein structure

The hemagglutinin (HA) structure of influenza A(H3N2) was generated utilizing SWISS-MODEL (https://swissmodel.expasy.org) (Waterhouse et al., 2018) along with the template model A/Victoria/361/2011(H3N2) (PDB ID: 4WE8). In a three-dimensional (3D) structure, the amino acid substitutions at the antigenic sites of HA proteins were illustrated utilizing the PyMOL Molecular Graphics System version 2.5.4 (Schrödinger, 2015). The predicted antigenic sites were obtained from previous studies (Popova et al., 2012; Wiley, Wilson & Skehel, 1981; Wilson, Skehel & Wiley, 1981; Wu et al., 2020).

N-glycosylation sites prediction

Utilizing the NetNGlyc-1.0 server (Gupta & Brunak, 2002), the N-glycosylation sites were predicted. All complete amino acid sequences of HA and NA from circulating A(H3N2) viruses and A/Victoria/361/2011 (H3N2; KM821347, KJ942682) were uploaded to the server. With a potential threshold higher than 0.5, asparagine occurring in the Asn-X-Ser/Thr sequences, where X can represent any amino acid except proline, was predicted to be N-glycosylated.

Evolutionary analysis and estimation of the influenza A(H3N2) virus substitution rate

Bayesian evolutionary analysis of molecular sequences utilizing Markov chain Monte Carlo (MCMC) in the Bayesian Evolutionary Analysis Sampling Trees (BEAST) program v.2.7.4 (Bouckaert et al., 2019) was used to estimate the evolutionary rates of the HA and NA genes. The evolutionary rate was analyzed using the birth-death skyline model (BDSKY) with the Hasegawa–Kishino–Yano (HKY) substitution model, which was determined using the maximum likelihood best-fit substitution model (ML) calculation in the MEGA X program (Kumar et al., 2018). The BEAST2 program was run with a strict clock model for 10 million MCMC chains and sampled every 1,000 frequencies. The substitutions per site per year were employed to measure the evolutionary rates. The prior distribution was calculated by estimating the evolution of the influenza A(H3N2) virus using the substitution rate of 4.84 × 10⁻³ substitutions per site per year (Jenkins et al., 2002; Phyu et al., 2022; Westgeest et al., 2014). Tracer v.1.7.2 (Rambaut et al., 2018) was used to represent the estimated evolution parameters derived from BEAST.

Prediction of vaccine effectiveness against influenza A(H3N2)

The predicted vaccine effectiveness (VE) was calculated using pEpitope calculator (version 2.1) (Bonomo, Kim & Deem, 2019) in the MATLAB program (The MathWorks Inc., 2022). The deduced amino acid sequences of the HA gene were submitted to measure the antigenic distance between the circulating viruses and the vaccine strains in the relevant years. The equation of the VE = (−3.32 × p_epitope + 0.66) × 100%, p_epitope is defined as the ratio of the difference between the dominant epitopes of the vaccine strain and the circulating viruses.

Statistical analysis

GraphPad Prism version 10.0.3 (GraphPad Software, La Jolla, CA, USA) was utilized to compute the significant difference values for the vaccine effectiveness through the use of the unpaired t-test followed by the Mann–Whitney test. The results are displayed as the mean ± standard deviation (SD). The p-value was calculated at a significance level of <0.05.

Determination of influenza A virus subtype in clinical samples

To determine the subtype of influenza A virus, a total of 443 samples that tested positive for influenza A viruses were analyzed using multiplex real-time RT-PCR. 147 samples (33.18%) of influenza A(H1N1)pdm09 viruses were found. Concurrently, influenza A(H3N2) viruses were found in 296 samples (66.82%), as indicated in Table 1. Both subtypes, influenza A(H3N2) and A(H1N1)pdm09, were present at the same time. However, influenza A(H3N2) was more prevalent than A(H1N1)pdm09 in each year from 2015 to 2018 (see Fig. S1). No simultaneous infections of both subtypes were identified. It is noteworthy that a high prevalence occurs from the rainy season (June-September) through the winter season (October-January) on an annual basis. Regrettably, recruitment efforts for the samples conducted in 2017 from September to December were unsuccessful.

Phylogenetic analysis of the influenza A(H3N2) viruses

To further characterize the viral HA and NA genes using bidirectional Sanger sequencing, influenza A(H3N2) positive respiratory samples were used. As a representative, two samples were chosen from each month in 2015, 2016, and 2018, while only one sample was obtained each month in 2017 because of the restricted sample size. The selected criterion for representative samples was mainly based on the Ct cycle (less than or equal to 28) from the subtyping influenza A viruses by multiplex real-time RT-PCR experiment. Therefore, a total of 65 samples were recruited from 296 influenza A(H3N2) virus samples (16 (2015), 21 (2016), 21 (2018), and 7 (2017)).

The complete sequence of the HA and NA genes was obtained from two distinct fragments (A and B), as reported in the study by Lee et al. Our study successfully obtained 35 samples (53.85%) for the entire HA and NA genes received 35 samples (53.85%). Altogether, 26 samples (40%), successfully obtained complete HA and NA gene sequences (Table 2). Together with the influenza A(H3N2) virus vaccine strains recommended by the WHO from 2010 to 2020, the sequences of the HA and NA genes from the influenza A(H3N2) positive samples were analyzed. Phylogenetic trees have been constructed for the HA (Fig. 1) and NA (Fig. 2) genes.

The A/Switzerland/9715293/2013 vaccine strain, which was categorized as clade 3C.3a, was distinct from all other influenza A(H3N2) viruses that were circulating in 2015. Nevertheless, they were all members of clade 3C.2a. The vaccine strain used in 2018, A/Singapore/INFIMH-16-0019/2016, was closely related to the circulating viruses in 2016 and 2017, which were classified under sub-clade 3C.2a1. On the other hand, the vaccine strain used in 2016–2017, A/Hong Kong/4801/2014, belonged to clade 3C.2a. 2018 revealed the discovery of two influenza A(H3N2) sub-clades, 3C.2a1b and 3C.2a2. Sub-clade 3C.2a2 was in the same clade as the vaccination strain used in 2019, A/Switzerland/8060/2017, whereas sub-clade 3C.2a1b was closely similar to the vaccine strain used in 2018 (Figs. 1 and 2).

Moreover, to investigate the relationship between the A(H3N2) virus in Thailand and the global influenza A(H3N2) virus. The global influenza strains from 2015 to 2018 were compiled into a phylogenetic tree together with the HA gene of the circulating viruses, based on the continents of Africa, Asia, Europe, North America, South America, and Oceania. As shown in Fig. S2, the results showed that each year, the viruses in Thailand were closely related to viruses found in Asia, North America, and Oceania.

Genetic variations in the HA and NA genes

Based on nucleotide sequencing, 35 sample sequences of each HA and NA gene were further analyzed. The nucleotide alterations encompassed both synonymous and non-synonymous mutations. Therefore, only non-synonymous mutations at the HA antigenic sites were extensively examined.

The MatGAT program was utilized to compute the amino acid similarities by comparing the deduced amino acid sequences of the HA gene found in circulating viruses and vaccine strains from the corresponding years. The outcomes for the years 2015, 2016, 2017, and 2018 were 97.5%, 98.4%, 98.8%, and 98.5%, respectively. HA antigenic site substitutions were identified at antigenic sites A, B, C, D, and E (Table 3). In 2015, amino acid substitutions were observed at all antigenic sites on the viruses, including A138S, R142G, and S144N (A); T128A, Y159S, and T160K (B); H311Q (C); S248T and K171N (D); and D62E and R80Q (E). However, two additional mutations were identified on antigenic sites A (T133N) and B (R197Q) in one circulating strain (A/Thailand/CU-MV34/2015). For 2016, most of the viruses exhibited the following mutations at the antigenic sites: S160K at site A; N137K and N187K at site D. The dominant amino acid substitutions identified in the 2017 viruses were antigenic site A (T151K) and antigenic site D (T137K and T187K). The viruses that exhibited the most prevalent mutations in 2018 were located at antigenic site A (T147K and G158K), antigenic site D (N137K, N187K, and A228T), and antigenic site E (E78K and R277Q). Due to most of the amino acid substitutions being changed in the amino acid group, the HA protein structure was observed. Homology modeling was used to construct the three-dimensional HA protein structure. As a constructed HA model, the influenza A/Victoria/361/2011(H3N2) virus was utilized for demonstrating the antigenic sites A–E. (Figs. 3A and 3B). The constructed HA model was substituted with mutant HA sites that have been identified in circulating viruses. In our study, the mutations seem improbable to affect the tertiary structure of the HA protein (Figs. 3C–3F).

The amino acid similarities of the NA gene were determined similar to the HA gene. The findings indicated significant resemblances of 98.7% (2015), 99.1% (2016), 98.3% (2017), and 98.9% (2018). The amount of amino acid substitutions in 2015, 2016, 2017, and 2018 were 10, 12, 18, and 7 positions, respectively (Table 4).

The prediction of N-glycosylation site in the HA and NA genes

Glycosylation is an important post-translational modification that significantly impacts viruses. Substitutions of amino acids may affect the glycosylation pattern. Therefore, the N-glycosylation sites located in the HA gene have the potential to impact not only the antigenicity of the HA protein but also the viral infectivity efficacy. The N-glycosylation sites were predicted by comparing whole amino acid sequences from 35 circulating viruses to those of the A/Victoria/361/2011. The findings revealed that all circulating viruses had one additional N-glycosylation site (N174) and lacked N160, in contrast to the thirteen N-glycosylation sites found in the HA gene of A/Victoria/361/2011 (Table 5). However, three viruses were identified that possessed distinctive N-glycosylation locations. A/Thailand/CU-MV34/2015 possessed one more additional N-glycosylation site (N112) and lacked N149. A/Thailand/CU-MV22/2017 lacked the N-glycosylation site N149. Lastly, two N-glycosylation sites (N142 and N149) were removed in A/Thailand/CU-MV27/2018. At antigenic site B of the HA gene, a non-synonymous mutation (N160S) was found, which is in line with the outcomes of the prior experiment. This mutation potentially contributes to the absence of the N-glycosylation site N160.

In the NA of A/Victoria/361/2011, seven predicted N-glycosylation sites were identified. Two additional N-glycosylation sites (N245 and N329) were identified in all influenza A(H3N2) viruses in Thailand (Table 6). Correlating with the previous result, the non-synonymous mutations in the NA gene (S245N and N329S) were identified. These mutations may have the potential to impact the addition of N-glycosylation sites N245 and N329.

The evolution rate of influenza A(H3N2) virus

To ascertain the influenza A(H3N2) virus’s evolutionary rate in the HA and NA genes, the Bayesian MCMC approach was used to estimate the evolutionary rate. In the case of the HA gene, the mean evolutionary rate was 3.47 × 10⁻³ substitutions per site per year with 95% highest posterior density (HPD) between 2.55 × 10⁻³ to 4.43 × 10⁻³ substitutions per site per year. For the NA gene, the mean evolutionary rate was 2.98 × 10⁻³ substitutions per site per year with 95% HPD between 2.17 × 10⁻³ to 3.90 × 10⁻³ substitutions per site per year.

The predicted vaccine effectiveness against influenza A(H3N2) from 2015 to 2018

According to the results of our phylogenetic analysis and genetic variations, vaccine protection may be affected by mutations occurring on the epitope of the HA gene. The pEpitope calculator was utilized to compute the predicted vaccine effectiveness (VE) against the circulating influenza A(H3N2) viruses. As shown in Fig. 4, the mean p_epitope values between the vaccine strain and the circulating viruses were 0.17, 0.05, 0.07, and 0.1 whereas the average and standard error of the predicted VE was 10.56 ± 12.28%, 48.48 ± 8.61%, 43.65 ± 8.94%, and 31.63 ± 9.88% in 2015, 2016, 2017, and 2018, respectively. The findings revealed that the predicted VE in 2015 was statistically significantly lower than those in 2016–2018.