Evolution, Transmission, and Pathogenicity of High Pathogenicity Avian Influenza Virus A (H5N8) Clade 2.3.4.4, South Korea, 2014–2016

During 2014–2016, clade 2.3.4.4 H5N8 high pathogenicity avian influenza virus (HPAIV) caused the largest known avian influenza epidemic in South Korea. Based on data from earlier H5N8 outbreaks, primitive H5N8 virus in South Korea was classified into five subgroups: C1, C2, C3, C4, and C5. The present study investigated the pathogenic and molecular epidemiologic characteristics of H5N8 viruses obtained from 388 cases of poultry farms and 85 cases of wild bird infections in South Korea during 2014–2016. Representative viruses of subgroups C1, C2, and C4 showed significant pathobiological differences in specific pathogen-free (SPF) chickens, with the H1731 (C1) virus showing substantially lower infectivity, transmissibility, and pathogenicity than the H2102 (C2) and H1924 (C4) viruses. Full genome sequence analysis showed the number of mutations that significantly increased in domestic duck-origin H5N8 HPAIVs compared to the viruses from gallinaceous poultry. These differences may have been due to the long-term circulation of viruses in domestic duck farms. The same mutations, at positions 219 and 757 of PB1, independently evolving in the C0, C1, and C2 subgroups may have been positively selected, resulting in convergent evolution at the amino acid level. Bayesian discrete trait phylodynamic analysis (DTA) indicated multiple introductions of H5N8 HPAIV from wild birds into domestic poultry in various regions in South Korea. Following initial viral introduction into domestic duck farms in the western part of Korea, domestic ducks played a major role in viral transmission and maintenance. These findings highlight the need for continued genomic surveillance and pathobiological characterization of HPAIV in birds. Enhanced biosecurity in poultry farms should be implemented to prevent the introduction, maintenance, and spread of HPAIV.


INTRODUCTION
The H5 subtype of high pathogenicity avian influenza viruses (HPAIVs) belonging to the Goose/Guangdong (Gs/Gd) lineage have been circulating since 1996 and evolving into genetically diverse clades and subclades (1)(2)(3). The H5N8 subtype of Gs/Gd HPAIV belonging to the subclade 2.3.4.4 was initially identified in domestic ducks from eastern China in 2010 (4) and novel reassortant H5N8 HPAIVs were detected in live poultry markets in eastern China in late 2013. Subsequently, reassortant H5N8 HPAIVs were introduced into South Korea and Japan in early 2014 (5)(6)(7). In late 2014, several countries in Europe and North America experienced an invasion of HPAI H5Nx viruses (8)(9)(10)(11). The transcontinental spread of these viruses was associated with dissemination in multiple directions by migratory wild birds (12).
In  (13). The fifth HPAI outbreak started in January 2014 and lasted for 28 months (14). During this outbreak, South Korea experienced four waves of H5N8 outbreaks with 393 cases of poultry farms and 58 cases of wild bird infections. Moreover, two genetically distinct groups of clade 2.3.4.4 H5N8 viruses were identified: Buan2-like and Gochang1-like. The Buan2-like H5N8 viruses predominated during the epidemic and further diverged into five distinct subgroups: C1, C2, C3, C4, and C5 (8,14). After the first wave, from January to July 2014, subgroups C1 and C5 were detected mainly in poultry farms in South Korea, whereas the other three subgroups were detected in wild birds and spread over long distances via wild birds including subgroups C2 in South Korea and Japan, C3 in North America and Japan, and C4 in Europe, Japan, and South Korea (8)(9)(10)(11). The C2 and C4 subgroups were maintained in wild birds and re-introduced into South Korea during the fall migration season of wild birds in 2014.
Little is currently understood about the evolution and spread of H5N8 HPAIV during the 2014-2016 outbreak in South Korea. Previous studies investigating the origin and transmission of H5N8 HPAIV in South Korea included limited datasets, such as partial genome sequences of early isolates or selected viruses (8,14,15). The present study analyzed the complete genome sequences of H5N8 HPAIVs identified in 388 (of 393 total cases) poultry farms and 53 (of 58 total cases) wild bird cases during the 2014-2016 outbreak in South Korea. The pathogenic and molecular epidemiologic characteristics of the viruses were determined to identify the underlying causes of the prolonged outbreak and to reconstruct the evolutionary and transmission dynamics of H5N8 HPAIVs in South Korea.

Samples
Fecal samples, bird carcasses, and oropharyngeal and cloacal swabs were collected from poultry farms and wild bird habitats in South Korea through active and passive surveillance from January 2014 to June 2016. Viruses were isolated by inoculating samples into 9-11-day-old specific-pathogen-free (SPF) embryonated chicken eggs after incubation for 4 days at 37 • C.

Genome Sequencing
Viral RNA was extracted from the infected allantoic fluids of SPF chicken eggs using Maxwell R RSC simplyRNA Tissue Kits (Promega, Madison, WI, USA). Complementary DNA was synthesized by reverse transcription reaction using Superscript III First-Strand Synthesis System (Invitrogen, Carlsbad, CA, USA) and eight segments of each virus were amplified by PCR using AccuPrime Pfx DNA Polymerase (Invitrogen, Carlsbad, CA, USA), as described (16). The complete genome sequences of 441 H5N8 HPAIVs isolated from 388 domestic poultry farms and 53 wild bird habitats were determined by nextgeneration sequencing with the Illumina Miseq platform using the Nextera DNA Flex Library Prep Kit (Illumina, San Diego, CA, USA) according to the manufacturer's instructions. Data were analyzed by CLC Genomics Workbench (Qiagen, Valencia, CA, USA). Nucleotide sequences were deposited under accession nos. EPI_ISL_157609-EPI_ISL_410213.

Genetic Subgrouping by Maximum-Likelihood Phylogenetic Analysis
The complete genomes of the 441 H5N8 HPAIVs sequenced in this study were subjected to phylogenetic analyses, along with 32 previously determined whole genome sequences of clade 2.3.4.4A H5N8 HPAIV identified in wild waterfowl in South Korea that had been deposited in the GISAID EpiFlu database (http://www. gisaid.org). Furthermore, our data set was added 20 sequences of the H5N8 viruses isolated from poultry and wild waterfowl in Japan. Maximum-likelihood (ML) phylogenies of the 493 concatenated complete genome sequences were generated using the Randomized Axelerated Maximum Likelihood (RAxML 8.2.12) (17) method, involving the general time reversible model of nucleotide substitution and the Gamma model of amongsite rate heterogeneity model with 1,000 bootstrap iterations. A genetic subgroup was defined as a monophyletic cluster of sequences with high bootstrap support (>70%).

Animal Experiments
The objective of this study was to evaluate the infectivity, transmissibility, and pathogenicity of the first detected viruses (index virus) of three subgroups (C1, C2, and C4 except for C5 with only four cases) of H5N8 HPAIVs in chickens. Viruses evaluated in chickens included A/broiler duck/Korea/H1731/2014 (subgroup C1, abbreviated H1731), A/mallard/Korea/H2102/2015 (subgroup C2, abbreviated H2102), and A/mallard/Korea/H1924/2014 (subgroup C4, abbreviated H1924). To determine the mean bird lethal dose (BLD 50 ), five 4-week-old SPF white leghorn chickens each were inoculated intranasally with serial 10-fold dilutions (10 3 to 10 6 mean egg infectious dose [EID 50 ]) of each virus. Viral transmissibility was evaluated by housing three virus-naïve chickens with the chickens challenged with 10 6 EID 50 /0.1 ml of each virus 8 h after infection. The chickens were monitored daily for clinical signs for 14 days. Oropharyngeal and cloacal swabs were collected at 1-7, 10, and 14 days post-infection (dpi). Serum was collected from surviving chickens at 14 dpi to determine seroconversion by hemagglutination inhibition test (HI) (18). The tissue tropism and dissemination of the viruses were evaluated by intranasally inoculating three chickens with 10 6 EID 50 /0.1 mL of each virus. These chickens were necropsied at 3 dpi and 12 organs (trachea, thymus, liver, pectoral major muscle, spleen, proventriculus, pancreas, cecal tonsil, lung, kidney, heart, and brain) were collected for viral titration. Virus titers were measured in DF-1 cells incubated with 10% tissue homogenates. Experiments in animals were reviewed and approved by the Institutional Animal Care and Use Committee of the Animal and Plant Quarantine Agency (APQA) (approval no: 2018-398 and 2018-412). All procedures were carried out in a biosafety level three facility at the APQA.

Shannon Entropy and Selective Pressure Analysis
Genomic variability between H5N8 HPAIVs derived from ducks and other hosts was evaluated by genome-wide Shannon entropy analysis using the entropy tool in the HIV sequence database (https://www.hiv.lanl.gov/content/sequence/ENTROPY/ entropy.html). Shannon entropy differences (Hdiff) were calculated for each gene segment, with differences in variability considered significant when p < 0.05.
Gene-and site-specific selection pressures for all segments of the H5N8 viruses isolated in South Korea during 2014-2016 were measured as the ratio of non-synonymous substitutions (dN) to synonymous substitutions (dS) per site. These analyses were performed using a combination of four methods: single-likelihood ancestor counting (SLAC), fixed-effects likelihood (FEL), fast unconstrained Bayesian approximation (FUBAR), and a mixed effects model of evolution (MEME) available at the Datamonkey online version of the Hy-Phy package (19)(20)(21). The selective pressures of each gene of the C0, C1, C2, and C4 subtypes were compared and statistically significant positive selection of amino acid sites was estimated (SLAC, FEL, MEME, p-value < 0.1; FUBAR, posterior probability >0.9).
DTA was also performed to determine the transmission dynamics between host species, including "Wild Bird" (ntax = 115 sequences), "Duck" (ntax = 294 sequences), "Chicken" (ntax = 71 sequences), and "minor poultry" (ntax =27 sequences). The Wild Bird, Duck, and Chicken datasets were subsampled using CD-HIT, based on nucleotide identities of 99.95, 99.87, and 100%, respectively (22). This subsampled DTA dataset consisted of 204 sequences, including 65 wild bird, 61 duck, 51 chicken, and 27 minor poultry sequences. For both analyses, Bayesian relaxed clock phylogenies were reconstructed using BEAST version 1.10.4 (23). The Hasegawa, Kishino, and Yano nucleotide substitution model with an uncorrelated log-normal distribution relaxed-clock method was used, along with a Gaussian Markov Random Field (GMRF) Bayesian skyride coalescent prior (24). The best supported viral transitions between discrete categories was identified by calculating the Bayes factor (BF) using SPREAD version 1.0.7 (25). A transition was defined as significant when BF >3 and posterior probability >0.5. The rate and number of viral transitions between discrete categories (Markov jump) and the time spent in a state between transitions (Markov reward) were estimated using stochastic mapping (26). The Markov Chain Monte Carlo (MCMC) was run in parallel for three chains, each with 100 million steps and samples across chains combined after 10% burn-in. The parameters, each of which had effective sample sizes >200, were analyzed with TRACER v1.5 (http://tree.bio.ed. ac.uk/software/tracer/). A maximum clade credibility (MCC) tree was generated using TreeAnnotator and visualized using FigTree 1.4.4 (http://tree.bio.ed.ac.uk/software/figtree/). The number of transition events and the discrete state proportion over time in monthly intervals were calculated based on analyses of posterior trees using the program PACT (http://www.trevorbedford.com/ pact).
The impact of viral source/sink sample sizes on Bayesian discrete inference datasets was determined using generalized linear model (GLM) analysis as an extension of DTA (23). GLM calculates the virus migration rates as a linear combination of coefficients and coefficient indicators, and predictors. The sample size of each discrete state was used to predict the viral transition rates of two separate DTAs. The coefficient indicator and coefficient describe if and to what degree each predictor contributes to the migration rate.

ML Phylogenetic Analysis
A ML phylogenetic tree of concatenated whole genome sequences of 493 H5N8 HPAIVs isolated in South Korea and Japan during 2014-2016 was generated to determine the genetic relationship among isolates (Figure 1). The viruses isolated from wild bird and poultry in South Korea were classified into multiple genetic clusters, including primitive (subgroup C0) and further classified (subgroups C1, C2, C4, and C5) viruses.
Viruses of subgroup C3, which were detected in North America, Taiwan, and Japan during 2014-2015, were not detected in South Korea.

Distribution in Geographic Region and Host
Of the 273 viruses belonging to subgroup C0 detected in South Korea from January to June 2014, 58 were detected in wild birds and 215 in domestic poultry, including in 162 ducks,

Pathogenicity and Transmissibility in Chickens
To investigate the differences in pathobiology among viral subgroups, chickens were intranasally inoculated with H1731, H2102, or H1924 virus, defined as the index viruses of the C1, C2, and C4 subgroups, respectively ( Table 1). All chickens inoculated with 10 6 EID 50 /0.1 ml of the H1731 virus died within 6 days with a mean death time (MDT) of 5 days. By contrast, all chickens inoculated with 10 6 EID 50 /0.1 ml of the H2102 and H1924 viruses died in less than 4 days with MDTs of 3.8 and 3.7 days, respectively. The BLD 50 of the H1731, H2102, and H1924 viruses were 10 4.4 , 10 4.5 , and 10 3.5 EID 50 /0.1 ml, respectively. Viral shedding of H1731 into the oropharynx and cloaca was detected from 2 to 6 dpi, whereas viral shedding of H2102 and H1924 were detected from 1 to 4 dpi (Figures 2A,B). The peak virus titers in oropharyngeal swabs were lower in H1731inoculated chickens (3.4 log 10 median tissue culture infectious dose [TCID 50 /0.1 ml]) than in chickens inoculated with H2102 (5.25 log 10 TCID 50 /0.1 ml) and H1924 (4.3 log 10 TCID 50 /0.1 ml). None of the surviving SPF chickens shed virus or showed an HI response.
Analysis of viral transmission showed that one of the three chickens cohoused with the H1731-inoculated chickens died and shed virus (1.3-3.3 log 10 TCID 50 /0.1 ml) into the oropharynx (Figures 2C,D). By contrast, all of the chickens cohoused with the H2102 and H1924-inoculated chickens died and shed virus into both the oropharynx (0.8-3.75 log 10 TCID 50 /0.1 ml) and cloaca (1.0-2.5 log 10 TCID 50 /0.1 ml). Chickens cohoused with the H1731-and H2102-inoculated chickens survived longer, with MDTs of 8 and 7.9 days, respectively, than chickens cohoused with H1924-inoculated chickens, with an MDT of 6.1 days.
All viral strains replicated in all collected tissues, including brain, heart, kidney, lung, cecal tonsil, pancreas, proventriculus, spleen, muscle, liver, thymus, and trachea (Figure 3). Viral titers in tissue samples were substantially lower in H1731-inoculated than in H2102-and H1924-inoculated chickens. In particular, the mean viral titers in the hearts of H1731-inoculated chickens were significantly lower than the viral titers in the hearts of H2102-and H1924-inoculated chickens (p < 0.001 each). In addition, viral titers in the kidney, lung, and cecal tonsil were significantly lower in H1731-than in H1924-inoculated chickens (p < 0.05 each), and viral titers in the thymus of H1731-inoculated chickens were significantly lower than titers in the thymus of chickens inoculated with H1924 (p < 0.01) and H2102 (p < 0.05) viruses.

Shannon Entropy Analysis and Selective Pressure
Calculations of Shannon entropy for all positions of the coding regions identified 53 dominant mutations in duck-derived H5N8 HPAIVs when compared with other host-derived H5N8 viruses ( Table 2). H1731 (C1) virus had multiple mutations with significance, whereas Buan2 (C0), H2102 (C2), and H1924 (C4) viruses had no related mutation. Shannon entropy analysis of H1731 (subgroup C1) virus identified 18 amino acid substitutions (T184A, A442S, and D678G in PB2; V219I in PB1; V323I, V327G, and T618K in PA; S141P, S163N, D376N, Q455K, and A528V in HA; L23I, S69N, and T329I in NA; N80S and D171E in NS1; and M14K in NS2). Entropy analysis of the Buan2 (subgroup C0), H2102 (subgroup C2), and H1924 (subgroup C4) viruses, however, did not identify any amino acid substitutions. All of the mutations were founded in subgroup C1 and took         (Figure 7A). Collectively, these results indicate that domestic ducks and minor poultry played major roles in the maintenance and transmission of H5N8 HPAIV in South Korea after its multiple introductions by wild birds. GLM analysis revealed that sample size had a negligible effect on potential bias toward the origin or destination state in this analysis (source indicator: 0.156, source coefficient: 0.05453, sink indicator: 0.43, sink coefficient: 0.185).

Transmission History Between Regions
The DTA among provinces in South Korea involved seven discrete nominal categories: CB, CN, GG, GS, JB, JN, and WB. GLM analysis showed that sample size had a negligible effect on potential bias toward the origin or destination state in this analysis (source indicator: 0.13, source coefficient: 0.0601, sink indicator: 0.388, sink coefficient: 0.0996). The DTA indicated multiple introductions of H5N8 HPAIV from wild birds into domestic poultry in various regions in South Korea (Figures 8-10; Table 4). DTA also showed that virus was transmitted from WB to JB (TR: 1

DISCUSSION
AIVs evolve through the gradual accumulation of mutations and genome reassortment. H5N8 HPAIVs likely had opportunities to reassort with various AIV gene pools during the 28 months of  (28). The representative challenge viruses of the C1, C2, and C4 subgroups differed significantly in infectivity, transmissibility, and pathogenicity in SPF chickens. H1731 (C1) had a longer MDT, significantly lower viral titers in tissues, and lower transmissibility in chickens compared with H2102 (C2) and H1924 (C4). The virulence of H1731 in chickens was also lower than that of the buan2 (the index virus of C0) (29). In particular, the amounts of H1731 virus shed into the oropharynx and cloaca of inoculated SPF chickens ranged from 10 1.5 to 10 3.4 (TCID 50 /0.1 ml) and from 10 1.5 to 10 2.8 (TCID 50 /0.1 ml), respectively, whereas the amounts of Buan2 virus shed into these tissues ranged from 10 4.2 to 10 6.3 (TCID 50 /0.1 ml) and from 10 3.1 to 10 7.5 (TCID 50 /0.1 ml), respectively. Viral titers in tissue samples were also lower in chickens inoculated with H1731 than with buan2 virus (30). These findings suggest that viruses of the C1 subgroup have lower adaptation and transmissibility in chickens. Because viruses of the C1 subgroup had circulated mainly in domestic duck farms, they likely adapted to domestic ducks. AIVs acquire diverse mutations during adaptation to a particular host (31,32). Shannon entropy analysis identified several dominant amino acid substitutions in H5N8 HPAIVs from domestic ducks compared with viruses from other species, with H1731 being the only challenge virus to have amino acid substitutions, further indicating that viruses of the C1 subgroup had become adapted to domestic ducks.
Long-term circulation of AIVs in a certain host species allows viruses to adapt to that species, affecting their virulence and host specificity via genetic evolution. Two amino acid residues (positions 219 and 757 in PB1) in C0, C1, and C2 were found to be under positive selection pressure, with Shannon entropy analysis showing significant mutations at these positions in H5N8 HPAIVs from ducks. Because viruses of the C0, C1, and C2 subgroups mainly circulated in domestic duck farms during the 2014-2015 outbreak, these results suggest that these viruses may have undergone positive selection in domestic ducks. Amino acid substitutions that occurred during the adaptation of H5 HPAIVs to a particular host were found to be related to viral pathogenicity in that host (33). H5N8 HPAIVs have been reported to be transmitted efficiently, with subclinical infection in domestic ducks contributing to the silent spread of virus to other hosts (34,35). The linkage between long-term circulation of H5N8 HPAIV and asymptomatic infection in domestic ducks during the 2014-2016 outbreak has not been fully determined. Analyses of the putative role of selected amino acid substitutions in domestic ducks would require a determination of the molecular mechanism underlying viral adaptation to host and changes in pathobiology.
The estimated mean tMRCA of all H5N8 HPAIVs identified in Korea was July 2013, indicating that the H5N8 HPAIVs of wild bird origin detected in South Korea most likely emerged during the breeding season of wild waterfowl in 2013. These findings were consistent with other studies on the estimated tMRCA of H5N8 HPAIV (36,37). Phylogeographic analysis suggests multiple introductions of the H5N8 HPAIVs from wild birds into JB province, followed by subsequent spread to other provinces in South Korea during the winters of 2013-2014 and 2014-2015. These findings highlight the major role of JN, JB, CN, and CB provinces in dissemination and persistence of the virus after its introductions from wild birds. The regional proportion of ancestral discrete states through time and Markov reward analyses also indicate that the H5N8 viruses were maintained predominantly in western South Korea, specifically in JN, JB, GG, CN, and CB provinces (Supplementary Figure 3A).
The western region of South Korea contains abundant habitats for wild waterfowl and a high density of domestic duck farms (8), creating an environment vulnerable to the transmission of HPAIVs (38,39). In addition, surveillance studies have suggested that domestic ducks are a major contributor to the spread and maintenance of HPAIVs in South Korea (40). Consistent with previous findings, the present study showed that the transmission of H5N8 HPAIV from domestic ducks to other host species occurred frequently, further indicating that domestic ducks play major roles in the maintenance, amplification, and spread of HPAIVs (38)(39)(40). Although a previous study reported unidirectional transmission of H5N8 viruses from wild waterfowl to domestic ducks (40), the present study found exchange of the viruses between domestic ducks and wild birds (Supplementary Figure 3B). Geographical transmission dynamics showed dissemination of the virus into wild birds from JN, an area with large numbers of wild bird habitats adjacent to domestic duck farms. These findings suggest the need for enhanced levels of surveillance and biosecurity measures at domestic duck farms in this region to effectively monitor and prevent the introduction and spread of HPAIV. Viral phylodynamic analysis also showed exchanges of HPAIVs between minor poultry populations and domestic ducks in South Korea, indicating that minor poultry may serve as reservoirs to maintain and disseminate HPAIV. Minor poultry are generally raised in small outdoor operations and backyards, and are usually marketed in live bird markets (LBMs), highlighting their importance as intermediary hosts in virus transmission under poor biosecurity conditions (41). LBMs include a wide variety of live poultry species, providing an ideal environment for the introduction, maintenance, and adaptation of viruses, as well as potential conditions for transmissions between duck farms and chicken farms (40).
Collectively, the present study presents unique molecular epidemiology and pathobiology of the H5N8 HPAIV outbreak in 2014-2016 in South Korea. Our data showed multiple introductions of the H5N8 HPAIV isolated in South Korea during 2014-2016 from wild waterfowl to poultry farms in multiple provinces. The virus, initially introduced into the western part of South Korea, which contains large populations of domestic ducks, was subsequently disseminated into other regions throughout the country. Furthermore, domestic ducks likely played a pivotal role in the persistent circulation of the virus, with minor poultry also serving  as a source population. In addition, sequence analysis and in vivo experiments support that the possible adaptation of H5N8 HPAIV in domestic ducks likely reduced its virulence in chickens. Enhanced genomic surveillance and pathobiological characterization of the viruses are essential for better understanding of HPAI epidemiology and the design of prevention strategies.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found at: GISAID EpiFlu database (http://www.gisaid.org) and accession no. EPI_ISL_157609-EPI_ISL_410213.

ETHICS STATEMENT
Experiments in animals were reviewed and approved by the Institutional Animal Care and Use Committee of the Animal and Plant Quarantine Agency (APQA) (approval no: 2018-398 and 2018-412).