In November 2017, two novel H5N6 viruses were identified in South Korea: A/duck/Korea/HD1/2017 (HD1) and A/mallard/Korea/Jeju-H24/2017 (Jeju-H24) [8]. Both HD1 and Jeju-H24 are closely related to European clade 2.3.4.4b (group B) H5N8 viruses in all gene segments except neuraminidase (NA), which is derived from a Eurasian low pathogenic avian influenza (LPAI) lineage. In December 2017, we identified A/duck/Korea/H35/2017 (H35) and nine more H5N6 viruses in Korean domestic duck farms and wild-bird habitats (Table 1, Supplementary Table 1). Whole-genome sequencing showed that HD1 and Jeju-H24 had 99.69% to 100% nucleotide-sequence identity, and H35 had 99.43% to 100% identity to the other viruses isolated in December 2017. Sequence identity was lower (96.96–99.38%) between HD1/Jeju-H24 and the other isolates.

To clarify the origins of novel H5N6 reassortants possessing genetic heterogeneity and to estimate the timing of reassortment events that led to their emergence, we performed genetic analyses based on time-scaled phylogenies, gene constellations, and molecular clocks with the available clade 2.3.4.4b H5Nx sequences and other N6 viruses from the Global Initiative on Sharing All Influenza Data (GISAID) and GenBank databases (Supplementary Table 2).

We estimated the time-scaled phylogenies, nucleotide substitution rates, and the most recent common ancestors (tMRCAs) using the Bayesian Markov chain Monte Carlo (MCMC) method implemented in BEAST v.1.8.4 [3]. Genomic sequences were selected on the basis of sampling locations and collection dates. Genotypes were defined on the basis of the gene segment-specific phylogenetic trees for clade 2.3.4.4b 17 H5N6 and 160 H5N8 HPAIVs. All genes of clade 2.3.4.4b H5 viruses were classified into different clusters according to tree topology and a posterior probability > 0.8. The genotype is the combination of the cluster assignment of eight gene segments.

By performing MCMC analysis, the nucleotide substitution rate estimated for the clade 2.3.4.4b H5 gene was 8.803 × 10⁻³ substitutions/site/year (95% highest posterior density [HPD] interval, 6.626 × 10⁻³ to 1.11 × 10⁻² substitutions/site/year), which was higher than a previous estimate for the HA gene of H5N1 viruses from China in 1996 to 2012 (4.378 × 10⁻³ substitutions/site/year) [15], but comparable to that of group A H5N8 viruses from South Korea in 2014 to 2015 (9.23 × 10⁻³ substitutions/site/year) [5]. We found a mean rate of nucleotide substitution in the N8 gene of 9.802 × 10⁻³ substitutions/site/year (95% HPD interval, 6.931 × 10⁻³ to 1.270 × 10⁻² substitutions/site/year). Previously, it was reported that the surface proteins of the influenza A virus evolve in a coordinated way under selection pressures to escape immune responses [11]. The high evolutionary rate of surface genes of these viruses resulting from positive selection by the immune system may have contributed to the observed genetic diversity.

Genes of clade 2.3.4.4b H5 viruses were classified into clusters as described above. The PB2, PB1, PA, and NP genes of these viruses formed six, three, five, and four distinct clusters, respectively (panels D–G in Supplementary Fig. 1), which identified eight H5N8 genotypes and three novel reassortant H5N6 genotypes (Fig. 1). The B0 genotype (including H5N8 viruses isolated from eastern China in 2013 and South Korea in 2014) is the prototype of the currently circulating group B H5 viruses. The B1 viruses were first identified in Russia–Mongolia in May 2016. In August 2016, B1.1 emerged in Russia, with a B1 backbone and a PA gene from Eurasian LPAI viruses. In September 2016, B3.1 (containing novel PB2, PA, and NP segments) emerged in Russia. Our results suggest that at least five H5N8 genotypes (B2, B2.1, B3.1, B3.2, and B3.3) circulated in Europe in the winter of 2016 to 2017 and contributed genes to novel H5N6 viruses.

HD1 and Jeju-H24 have seven genes derived from the B3.2 genotype, along with N6 genes from Eurasian LPAI virus. H35 and additional viruses have seven genes from the B3.1 genotype, along with Eurasian LPAI N6 genes. The N6 gene of HD1-like and H35-like viruses had the highest nucleotide identity with A/barnacle goose/Netherlands/2/2014 (H3N6) (> 97.3%). Compared to wild-bird isolates, we identified unique substitutions K or T133A in HA (H5 numbering) and A350T in NP protein from H35-like viruses isolated from domestic duck. Molecular dating estimates of the mean tMRCAs of HD1-like and H35-like viruses for each gene are shown in Table 2. These tMRCAs overlapped the breeding seasons of many wild birds inhabiting Eurasia. The tMRCA 95% HPD interval for the NA gene of the HD1-like virus was wider than that for the H35-like virus. Previous results suggested that reassortment events of H5N8 viruses in the Netherlands in winter 2016 to 2017 were completed by August 2016, and might have taken place in wild birds in the Russia– Mongolia region [1]. Our results indicated that reassortment of HD1-like viruses (B3.2.1) was complete by September 2016, whereas reassortment of H35-like viruses (B3.1.1) was complete by June 2017. Along with the 2.69% mismatches in NA sequences between HD1 and H35, these results suggest that the reassortment events in these viruses took place separately in wild birds in Eurasia.

In early 2017, the first H5N6 virus outside of Asia (A/chicken/Greece/39_2017/2017; Greece39) was detected in a backyard flock in Greece. It had 98.51% to 99.4% nucleotide identity with HD1. Despite genotypic similarity (Fig. 1), the MCC trees of six genes (panels A, D–G, and I in Supplementary Fig. 1) indicated that HD1 and Greece39 evolved along independent pathways. Furthermore, the differences in the tMRCA of the PB2 segment of Greece39 versus that of HD1 (Table 2) suggest that this genotype has a selective advantage over other competing genotypes, even though the HPD intervals overlap. By contrast, an entire genome sequence of an H5N6 virus isolated in Japan (A/mute swan/Shimane/ 3211A001/2017) and Taiwan (A/spoonbill/Taiwan/DB645/2017) in November 2017 had 99.39% to 100% similarity with HD1, and the two viruses formed a branch cluster together with HD1 in the MCC trees of all genes. In December 2017, H5N6 HPAIVs A/Great Black-backed Gull/Netherlands/1/2017 (GBBG1) and A/Great Black-backed Gull/Netherlands/2/2017 (GBBG2) were isolated in the Netherlands. GBBG1 had a B3.1.1 genotype, but GBBG2 had six genes derived from the B3.3 genotype, along with PB2 and N6 genes that clustered with LPAI viruses identified in Eurasia. These results demonstrated that H5N6 reassortants of various genotypes (B3.1.1, B3.2.1, and B3.3.1) originating from clade 2.3.4.4b H5N8 viruses appeared in East Asia and Europe in late 2017. Migratory birds moving within flyways in Eurasia overlap in circumpolar arctic and sub-arctic breeding areas [12]. Novel H5N6 reassortant viruses seem to have spread between migratory bird species in these breeding sites, subsequently disseminating in different directions during the 2017 migration period.

Since group B1 H5N8 HPAIVs were first detected, these viruses seem to have undergone reassortment with Eurasian LPAI viruses continuously. As a result, at least six genetically distinct H5N8, two H5N5 [1], and one H5N6 reassortant viruses were detected in Europe, Eastern Asia, and North Africa in the winter of 2016 to 2017. Although clade 2.3.4.4 HPAIVs have a tendency to partner with different NA subtypes, H5-NA matching with N6 was only observed in groups C and D until early 2017 [212], following reassortment from precursor H5Nx and H6N6 viruses [14]. The interaction and functional balance between HA and NA in AIVs have crucial roles in viral replication efficiency and fitness [10]. The simultaneous appearance of genetically distinct H5N6 reassortants suggests that this HA/NA genetic compatibility might enhance the replication capacity or the viral fitness related to the H5-N8 genetic combination.

In conclusion, our results show that at least three distinct reassortants of clade 2.3.4.4b H5N6 HPAIVs appeared in poultry and wild birds in East Asia and Europe in late 2017. Novel H5N6 reassortants potentially arose independently, spread through migratory birds at their breeding sites, and then disseminated in different directions during the 2017 migration. These findings provide evidence for more frequent reassortment of clade 2.3.4.4b H5 HPAIVs, which may increase their capacity to disseminate worldwide. Systematic surveillance in the Palearctic breeding grounds of Eurasia and timely analysis of the available data are essential to determine viral genetic diversity and to provide early warnings of novel reassortants of clade 2.3.4.4 HPAIVs.