A systematic study towards evolutionary and epidemiological dynamics of currently predominant H5 highly pathogenic avian influenza viruses in Vietnam

This study aimed to elucidate virus, host and environmental dynamics of Vietnamese H5 highly pathogenic avian influenza viruses (HPAIVs) during 2014–2017. Epidemiologically, H5 HPAIVs were frequently detected in apparently healthy domestic and Muscovy ducks and therefore these are preferred species for H5 HPAIV detection in active surveillance. Virologically, clade 2.3.2.1c and 2.3.4.4 H5 HPAIVs were predominant and exhibited distinct phylogeographic evolution. Clade 2.3.2.1c viruses clustered phylogenetically in North, Central and South regions, whilst clade 2.3.4.4 viruses only detected in North and Central regions formed small groups. These viruses underwent diverse reassortment with existence of at least 12 genotypes and retained typical avian-specific motifs. These H5 HPAIVs exhibited large antigenic distance from progenitor viruses and commercial vaccines currently used in poultry. Bayesian phylodynamic analysis inferred that clade 2.3.2.1c viruses detected during 2014–2017 were likely descended from homologous clade viruses imported to Vietnam previously and/or preexisting Chinese viruses during 2012–2013. Vietnamese clade 2.3.4.4 viruses closely shared genetic traits with contemporary foreign spillovers, suggesting that there existed multiple transboundary virus dispersals to Vietnam. This study provides insights into the evolution of Vietnamese H5 HPAIVs and highlights the necessity of strengthening control measures such as, preventive surveillance and poultry vaccination.

In 2001, H5N1 subtype highly pathogenic avian influenza viruses (HPAIVs), of the A/goose/Guangdong/1/1996 (Gs/GD) lineage, were detected for the first time in geese in Vietnam 1 . By December 2003, other Gs/GD-lineage H5N1 viruses had been detected in the North of Vietnam and by end-2004 had caused large poultry outbreaks in 57 out of 64 provinces [1][2][3] . Since then, H5 HPAIVs have remained enzootic and primary causative agents of several thousand outbreaks in poultry 4,5 . In addition, Gs/GD-lineage H5 HPAIVs have posed great public health concerns. To date, 127 H5 HPAIV clinical cases have been reported in Vietnam of which 64 were fatal (approximately 50% fatality rate) 6 2,9 . Due to persistent circulation of various H5 HPAIVs in Vietnam and the broader regions, reassortment events resulted from exchanging gene segments between H5 HPAIVs and other avian influenza viruses (AIVs) increases genomic diversification. Approximately 56 distinct genotypes, ranging from VN1-VN56, of Vietnamese H5 HPAIVs have been identified [7][8][9][10] .
To reduce virus loads in poultry and human health risk, mass poultry vaccination against H5 HPAIVs was implemented in Vietnam in mid-2005 11,12 . Vaccination was followed by a 12 month period in which no new outbreaks in poultry were reported and 18 months before a human case was reported 11,13,14 . Currently, two commercial vaccines are used in poultry vaccination programs. The first contains a clade 1 antigen which is produced from inactivated split-virion A/Vietnam/1194/2004 (H5N1). The second is based on a clade 2.3.2 antigen which produced from inactivated rgA/duck/Guangdong/S1322/2010 (H5N1) 5 . Updating of vaccine antigens has not been performed since the clade 2.3.2 vaccine was introduced 2 . Though mass vaccination has been implemented, it was never anticipated that vaccination would eliminate Gs/GD-lineage H5 HPAIVs from Vietnam. Vaccine failure can occur by several factors including insufficient doses, concurrent immunosuppressive disease and antigenic variants 13 . Experimentally, poultry vaccination accelerates selection of antigenic variants of H5 HPAIVs 15 . In addition, novel strains of Gs/GD-lineage H5 HPAIVs that have a poor antigenic match with existing vaccines have been imported to Vietnam 2 . From a general epidemiological perspective, H5 HPAIV persistence is facilitated by agro-ecological condition in which mobile free-ranging ducks, low biosecurity of backyard/small-scale production sectors and poorly bio-secure live bird markets (LBMs) are all considered as key risk factors [16][17][18] . Thus, comprehensive knowledge of viral, host and environmental factors affecting H5 HPAIV evolution is essential for successful control and prevention of H5 HPAIVs. Therefore, this study aimed to investigate virological and epidemiological features of contemporary Vietnamese H5 HPAIVs for the 2014-2017 period, using the results of our active surveillance, global surveillance data and longitudinal literature review.

Results
Surveillance programs and identifying specific hosts of H5 HPAIVs. Profiles gathered from passive and active surveillance programs in poultry in Vietnam were reviewed. The choropleth map in Fig. 1A illustrated space-time distribution of provinces reporting H5 HPAIV outbreaks as denoted via passive notification and provinces conducting active surveillance for the period 2014-2017. Outbreaks of disease associated with H5 HPAIVs were reported in all regions (North, Central and South) in each year, indicating that H5 HPAIVs persisted or occurred in a wide geographic and temporal range. Some variation in the number of provinces reporting outbreaks was apparent, ranging from 8 in 2016 to 35 in 2014 (Fig. 1A and Supplementary Table S1).
Active surveillance programs were deployed in all three regions of Vietnam. The number of H5 HPAIV strains detected in each province that were available in public influenza virus databases were co-plotted in the choropleth map. Within the 4-year period, a substantial number of Vietnamese H5 HPAIV strains (n = 192) were deposited in the public databases. More viruses were deposited in 2014 compared to 2016 in line with the increase in the number of outbreaks (Fig. 1A). These results demonstrate the importance of surveillance for monitoring H5 HPAIVs in Vietnam and its contribution to global influenza surveillance.
Regarding our work, 183 H5 HPAIVs were isolated from 12,585 swab samples from our surveillance program (Supplementary Table S2). Of these, 77 geographically and temporally representative isolates were selected for sequencing (Supplementary Table S3). The results from virus isolation in our active surveillance of apparently healthy birds showed that domestic ducks (51 isolates) and Muscovy ducks (72 isolates) were the species from which H5 HPAIVs were most frequently isolated (Fig. 1B). Likewise, the number of infected birds with other Vietnamese H5 HPAIV strains deposited in global surveillance databases (unknown whether passive or active sources) showed a similar pattern ( Supplementary Fig. S1). Therefore, domestic and Muscovy ducks are preferred species for detection of H5 HPAIVs in active surveillance.  Supplementary Fig. S2A). Similarly, phylogenetic relationships were observed for the NA gene segment. N1 NA genes were inferred to be descendants of N1 genes of H5N1 HPAIVs that previously existed in Vietnam and mostly grouped by the three geographic regions ( Supplementary Fig. S2B Supplementary Fig. S2C).

Distinct phylogeography and genomic diversity of contemporary
Phylogeny of the other six internal gene segments were assigned following previous systems 7-10 , each gene segment (i.e., PB2, PA, NP, M or NS) formed into 2-3 clusters ( Supplementary Fig. S2D,F-I). PB1 genes were likely www.nature.com/scientificreports www.nature.com/scientificreports/ more diverse and formed four distinct clusters in the phylogenetic tree ( Supplementary Fig. S2E). In addition, we identified novel gene cluster(s) in the PB2, PB1, PA, NP and NS gene segments (Supplementary Fig. S2D-G,I).
On the basis of viral reassortment capacity, a total of 132 H5 HPAIV strains (63 strains from our surveillance), which have full genome sequences, were grouped into 12 genotypes, four of which were identified for the first www.nature.com/scientificreports www.nature.com/scientificreports/ time in this study ( Fig. 2A and Supplementary Table S5). Clade 2.3.2.1c viruses had a greater number of genotypes (n = 7) compared to the number of genotypes in clade 2.3.4.4 (n = 4). The result indicates that Vietnamese H5 HPAIVs have continued its diverse genomic evolution although it is not clear whether these viruses have reassorted in Vietnam or elsewhere.
To study genotypic transition of Vietnamese H5 HPAIVs, genotype identification from previous studies for the period 2010-2013 9,10 was collated in our analysis. Figure 22 . The remaining 64% of the viruses lacked 11 amino acids in their stalk region (N6 numbering). It has been demonstrated that NA stalk truncation has arisen from evolutionary adaptation of H5 HPAIVs from wild aquatic birds to terrestrial poultry. Therefore, it might be concluded that recent Vietnamese H5 HPAIVs have retained typical avian-specific motifs and high adaptability in terrestrial avian species although most of these viruses were isolated from aquatic poultry.
Residues in the putative antigenic regions in the HA were highly variable. Antigenic region B had more mutated sites than the A, C, D and E regions (Table 1). Particularly, positions with highly frequent mutations were found at the region A: 125, 144 and 145; region B: 128, 129, 158, 159, 166, 167, 188 and 193; and region C: 271. These highly variable sites were presumably associated with antigenic variation of H5 HPAIVs.

Antigenic diversity of H5 HPAIVs and antibody titers of field vaccinated poultry antisera.
Antigenicity of Vietnamese H5 HPAIVs representing all clades and subclades detected in this study was characterized using cross hemagglutination inhibition (HI) test ( To examine the antigenic evolution of Vietnamese H5 HPAIVs, representative viruses were cross-reacted with laboratory reference antisera and antisera against the commercial clade 1 and 2.3.2 vaccines produced in laboratory condition using HI test. The cross-HI titers were used to produce antigenic cartography. Many of the inactivated viruses and homologous antisera formed independent groups with distances of at least approximately two antigenic units (4-fold HI difference) following distribution of each H5 HA gene assignment (Fig. 3A). For instance, the commercial inactivated Ck/Scot/1959 virus, the reference antigen currently used for assessing vaccination coverage in Vietnam, displayed a large antigenic distance from 2. www.nature.com/scientificreports www.nature.com/scientificreports/ To evaluate antibody titers of vaccinated birds in the field, 60 field antisera samples from vaccinated domestic poultry were tested with the reference antigen and representative inactivated Vietnamese H5 HPAIVs. Detectable titers of the 44 antisera samples were shown in Fig. 3B and Supplementary Table S8. The majority of field antisera exhibited protective HI titers (>8HI) 23    www.nature.com/scientificreports www.nature.com/scientificreports/ their first detection in Vietnam and homologous clade viruses in Laos, Cambodia and southern neighboring provinces of China were included in the analysis. Phylogenetically, clustering of the two clade HA segments in the BP trees was highly identical to ML trees on the basis of clade assignment as above described ( Supplementary  Fig. S3A Supplementary Fig. S4B). Afterward, at least two more introductory linkages were spatiotemporally traced from the southern bordering region of China to northern Vietnam in 2014, 2015 and 2017 (Fig. 4B)

Discussion
Vietnam is regarded as being endemically infected with Gs/GD-lineage H5 HPAIVs. Various control measures have been taken but virus elimination is not regarded as possible in the medium term. Studies aimed to identify viral, host and environmental dynamics of recently circulating H5 HPAIVs are essential to reinforce provisional control measures. The present study was able to systematically elucidate evolutionary and epidemiological dynamics of recent Vietnamese H5 HPAIVs using interdisciplinary analysis.
Regarding the surveillance of H5 HPAIVs in Vietnam, our spatial descriptive analysis illustrated that active and passive surveillance system in Vietnam allowed detection of viruses and disease outbreaks although it also demonstrated that some gaps exist that could result in non-detection of viruses. This work significantly contributed to global influenza surveillance via publishing substantial genetic data of detected viruses (Fig. 1A and www.nature.com/scientificreports www.nature.com/scientificreports/ Supplementary Table S3). However, it is also critical to discuss some constraints in the system. For instance, in the passive surveillance program, genetic information of H5 HPAIVs that caused outbreaks remained limited. The majority of virus sequences available in public databases were derived from active surveillance programs. A human case infected with a clade 1.1.2 virus was reported in a southern province in 2014, but no outbreak case and genetic information of clade 1.1.2 virus were reported in poultry, posing uncertainty about the actual circulation of viruses in this area (Fig. 1A). These examples emphasize the need to further develop surveillance, particularly virus investigation post surveillance and more focus on high risk areas. In addition, close collaboration between veterinary and human health sectors is needed for effective monitoring and control of AIVs.
In our surveillance, H5 HPAIVs were largely and frequently detected in apparently healthy domestic and Muscovy ducks (Fig. 1B). Role of domestic ducks for H5 HPAIV persistence has been demonstrated 24,25 , but the role of Muscovy ducks remains uncertain. Previous studies have reported that Muscovy ducks were more susceptible to Gs/GD-lineage H5 HPAIVs and vaccination was less effective in preventing infection compared with domestic ducks 26,27 . This result reinforces the epidemiological importance of domestic ducks and justifies further exploration of the role of Muscovy ducks. They remain good targets for active surveillance (Fig. 1B and Supplementary Fig. S1). Moreover, all H5 HPAIVs in our surveillance were isolated in LBMs and LBM-like locations, indicating that LBMs play important roles for virus persistence and dissemination 18,28 . These findings confirmed relatively constant epidemiology of H5 HPAIVs in regards that abundance of waterfowl species and live poultry trading activities facilitate H5 HPAIV persistence and these elements should be primary targets of the surveillance and control system.
In this study, we extended the genetic analysis of H5 HPAIVs to include all Vietnamese viruses that were uploaded in global influenza surveillance databases during the study period. This analysis allowed us to identify novel cluster(s) in the viral segments and elucidate the distinct phylogeography of two predominant clade 2.3.2.1c and 2.3.4.4 viruses ( Fig. 2A and Supplementary Fig. S2). In addition to high variability of the HA gene segment, the entire genomic backgrounds of H5 HPAIVs were divergent. A total of four novel genotypes out of 12 circulating genotypes were identified during the period 2014-2017 ( Fig. 2A). Despite multiple reassortment events from exchanges of internal gene segments, only N1 and N6 NA gene segments were consistently assembled with clade 2.3.2.1c and 2.3.4.4 H5 HA segments, generating two high risk H5N1 and H5N6 subtype viruses, respectively. These combinations of the two surface proteins HA/NA probably affords a fitness advantage for H5 HPAIV persistence in Vietnam. From our longitudinal review analysis, genotype transition of Vietnamese H5 HPAIVs was highly dynamic; several genotypes (VN48, VN52, VN53, VN54 and VN55) frequently emerged and then gained dominance from 2012 to 2017 (Fig. 2B). Once a genotype disappeared it did not tend to reappear, suggesting that genomic transition of Vietnamese Gs/GD-lineage H5 HPAIVs followed continuing process in accordance with existence of the HA assignments.
Molecular characterization via viral protein sequences provides information on likely phenotypic characteristics of viruses. Generally, H5 HPAIVs have remained adapted to domestic terrestrial birds and this was also the case with recent H5 HPAIVs in Vietnam which did not contain any mutation makers associated with mammalian adaptation and/or virulence (Table 1). Nonetheless, the threats to human health are always present; in fact, H5 HPAIVs possessing typical avian phenotypes triggered human infections in 2014 (Fig. 1A) [30][31][32] . This implies complexity of virological and epidemiological interactions and gaps of our understanding H5 HPAIV human infection in nature.
Our antigenic analysis results were highly consistent with previous studies, indicating that H5 HPAIVs recently detected in Vietnam had considerable antigenic variation 9,33,34 . The two predominant clade 2.3.2.1c and 2.3.4.4 H5 HPAIVs from 2014 to 2017 antigenically differ from each other and from clade 1.1 H5 HPAIVs (Table 2). Importantly, these predominant H5 HPAIVs had a large antigenic distance with progenitor H5 HPAIVs, the reference antigen and clade 1 vaccine. Recent clade 2.3.2.1c H5 HPAIVs also exhibited minor antigenic variation with clade 2.3.2 vaccine currently used in poultry in Vietnam (Fig. 3A). Thus, current vaccines should be reassessed through challenge trials. Clade 1 vaccine was unlikely to provide appropriate protection from disease 9,33,35 . Clade 2.3.2 vaccine is expected to be effective against recent 2.3.2.1c H5 HPAIVs but this needs to be monitored closely as some antigenic drift was apparent. Besides, vaccination assessment should be performed using both vaccine antigens and field strains and old reference antigens such as Ck/Scot/1959 should be changed to match the antigenic characteristics of circulating viruses.
Little is known about actual immunity of vaccinated poultry in the field against circulating viruses. Therefore, we assessed the field antisera of vaccinated poultry with the reference antigen and recent H5 HPAIVs (we have assumed that the titers are the result of vaccination and not field exposure but this possibility needs to be recognized in interpreting the data). Our results were in line with other reports that vaccination achieved sufficient coverage and protective antibody titers against the reference antigen and classical H5 HPAIVs (Fig. 3B). However, vaccination likely reduced effectiveness to recently circulating viruses 35,36 . Supposedly, 2.3.4.4 viruses likely become more outspread due to its greater ability to escape from current vaccine-induced immunity. Therefore, it is necessary to reconsider selection of vaccine strains that have close antigenic match to currently circulating H5 HPAIVs to enhance vaccination effectiveness.
Strikingly  37,38 . These results underlined the notion that transboundary epidemiological linkages for virus importation to Vietnam exist. In fact, a live poultry movement network between Vietnam and other neighboring countries, including China and Cambodia likely accelerates divergence of Vietnamese H5 HPAIVs 39 . Our speculation is supported by the notion that emergence/re-emergence of H5 HPAIVs in Vietnam is the primary consequence of transboundary spillovers 2 . Yet, we avoided to deeply interpret outwards dispersal linkages derived from Vietnam since interpretation lacks supportive evidence and carries biases. For instance, we only computed H5 HPAIV strains in several Chinese provinces bordering Vietnam, and inferred dispersal linkages of the studied Chinese viruses may be modulated if other geographic Chinese viruses were to be included. However, it is worthy to restate that virus dispersals from China to Vietnam were highly certain ( Fig. 4B and Supplementary Fig. S3B). Finally, we highlight the critical epidemiological role of live poultry movement, particularly cross-border poultry movement, to the transboundary dissemination of H5 HPAIVs and the necessity for preventive control against H5 HPAIV spillovers, including risks of introducing H7 HPAIVs 40 .
Even though we attempted to provide a robust elucidation of H5 HPAIVs in Vietnam, this study remains some unavoidable limitations. First, for instance, the present study investigated all Vietnamese H5 HPAIVs available via public databases but might not capture circulation of all H5 HPAIVs in the field. Second, we were unable to definitively clarify determinants for annual outbreak occurrence (i.e., low outbreak occurrence and virus detection in 2016); this may be due to interconnected effects of viral-agro-ecological factors. Thus, more systematic studies are needed to identify the determinants for outbreak occurrence of H5 HPAIVs.
In summary, the present study indicated that H5 HPAIVs persist in Vietnam and new virus strains are being introduced. During 2014-2017, clade 2.3.2.1c and 2.3.4.4 viruses were predominant and exhibited a large genetic and antigenic diversity. Considerably, these viruses were antigenically distinct from current vaccines and reference antigen used in poultry vaccination programs. Thus, this study highlighted an urgent necessity for reinforcing the provisional control measures in Vietnam including (1) strengthening surveillance for monitoring evolution of H5 HPAIVs, in particular more focus on surveillance and strict poultry movement control in border areas to prevent invasion of H5 HPAIV spillovers and (2) enhancing effectiveness of vaccination programs through updating vaccines and reference antigen to antigenically match currently circulating H5 HPAIVs.

Materials and Methods
Sampling, virus isolation and sequencing in our surveillance program. We conducted active surveillance of AIVs in the period 2014-2017 in representative provinces across the three geographic regions of Vietnam, North, Central and South. A total of 12,585 swabs were sampled from apparently healthy domestic birds and environment in farms, LBMs and other live poultry marketing places (LBM-like locations). The presence of AIVs was tested using M gene detection qRT-PCR. Virus isolation of M gene positive samples using chicken embryos and subtyping were performed following standard protocols 41 . Representative H5 HPAIVs based on geography and year were selected for full genome or single RNA segment sequencing using next generation or conventional Sanger sequencing, respectively 15 . H5 HPAIV strain details and sequence profiles were registered in GenBank/EMBL/DDBJ to obtain the accession numbers listed in Supplementary Table S3.
Data collection of active and passive surveillance programs. H5 HPAI outbreak occurrence in poultry, following passive surveillance in Vietnam during the study period, was obtained from the EMPRES-i, FAO database 42 . Active surveillance program profile of AIVs in Vietnam was obtained from the annual reports of animal diseases by the national Department of Animal Health (DAH) and collated with notifiable outbreak cases as surveillance entity 5 .
Entire Vietnamese H5 HPAIV strains detected for the period 2014-2017 were obtained from the Global Initiative on Sharing All Influenza Data (GISAID) and Influenza Research Database (IRD). Background details of each strain such as infected species, sampling origins and years were available within these databases or shared by our collaborators listed in the acknowledgements section.
Phylogenetic analysis, genotyping and molecular characterization. RNA segment sequences of the 192 representative H5 HPAIVs from our surveillance program and other H5 HPAIVs collated from GISAID and IRD were subjected to phylogenetic analyses. Nine sequence partitions derived from eight viral RNA segments (N1 and N6 NA gene segments were computed separately) were separately aligned using Muscle v3.8.3 43 . The alignments were manually inspected and poorly aligned sequence regions were removed using BioEdit v7.0.5 44 . The best-fit substitution model for each partition was selected based on the Bayesian Information Criterion using jModelTest v2.1.10 45 . Supplementary Table S4 summarizes sequence partition and modeled phylogenetic parameters with general time-reversible (GTR) substitution model following a gamma distribution. An ML tree for each partition was constructed using RAxML v8.2.4 46 with a resampling process of 1,000 rapid replicates. All ML trees were rooted with a corresponding gene segment of the Gs/GD precursor as expected topology. Phylogenetic trees were annotated and visualized using Figtree v1.4.3 47 or ggtree package 48 implemented in R v3.4.4 49 . Clade assignments and genotyping were determined from common phylogenetic relationships observed among the eight gene segments, following the standard criteria of WHO/OIE/FAO nomenclature and previous classification reports 7-10,21 . (2019) 9:7723 | https://doi.org/10.1038/s41598-019-42638-4 www.nature.com/scientificreports www.nature.com/scientificreports/ HI assay, antigenic cartography and antibody assessment. MDCK plaque-cloned viruses representing each H5 HA clade assignment were obtained for antigenic analysis. Antigenicity of these cloned H5 HPAIVs within and between clade assignments was determined by cross-HI test using a panel of chicken antisera as previously described 29,41 . To obtain an overall picture of antigenic evolution, the contemporary Vietnamese viruses, progenitor viruses including A/Hong Kong/483/1997 (H5N1) clade 0 virus; A/tern/South Africa/1961 (H5N3) reference virus; and the A/chicken/Scotland/1959 (H5N1) reference antigen were included for antigenic cartography. Intentionally, all the viruses were formalin-inactivated as a comparable form to the commercial reference antigen purchased from APHA Scientific, UK. These antigens cross-reacted with several reference laboratory antisera 29,41 and antisera collected from two chickens and two ducks immunized with commercial clade 1 and 2.3.2 vaccines under laboratory condition. Each chicken and duck was intramuscularly immunized with 0.5 mL clade 1 or 2.3.2 vaccines as manufacturer's instruction. The resulting data containing cross-HI titers was loaded to the web-based software (http://www.antigenic-cartography.org/) to obtain x/y coordinates of each antiserum and antigen. Antigenic cartography was then visualized using R v3.4.4.
Field antisera from chickens and ducks vaccinated with commercial vaccines in Vietnam were collected to test antibody titers. Vaccination status of individual birds was confirmed by farmers and local veterinarians, as described in Supplementary Table 7. Field antisera were transferred to Japan and tested using the HI test as described above.
Spatiotemporal phylodynamic analysis using Bayesian method. A 2,9,10 ). HA sequences of H5 HPAIV strains in Cambodia, Laos and southern neighboring Chinese provinces were combined with our HA sequence partition of each clade. The sampling dates and locations for each sequence were provided to BEAST to estimate space-time divergence. The lognormal uncorrelated relaxed clock model, the exponential growth tree prior and other parameters were used and described in Supplementary Table S4. For each clade dataset, the Bayesian Markov Chain Monte Carlo (MCMC) chains (one chain per dataset) were run for 90 million states, sampled every 5,000 iterations in a web-based CIPRES Science Gateway v3.3 51 . Convergence of the chains was assessed using Tracer 1.7 52 . Maximum clade credibility (MCC) trees were summarized using Tree Annotator v1.8.4; with 10% of the first sampled trees removed. Each MCMC analysis was independently performed at least twice and compared to ensure reproducibility. The final MCC trees containing space-time divergence were converted to a keyhole markup language file using SPREAD v1.0.7 53 . Then, spatial-temporal dispersal linkages of the two predominant H5 HPAIV population were visualized using several spatial visualization packages (ggplot2, ggmap, rgeos, sp, and maptools) implemented in R v3.4.4.