Dynamics of invasion and dissemination of H5N6 highly pathogenic avian influenza viruses in 2016-2017 winter in Japan.

Large highly pathogenic avian influenza (HPAI) outbreaks caused by clade 2.3.4.4e H5N6 viruses occurred in Japan during the 2016-2017 winter. To date, several reports regarding these outbreaks have been published, however a comprehensive study including geographical and time course validations has not been performed. Herein, 58 Japanese HPAI virus (HPAIV) isolates from the 2016-2017 season were added for phylogenetic analyses and the antigenic relationships among the causal viruses were elucidated. The locations where HPAIVs were found in the early phase of the outbreaks were clustered into three regions. Genotypes C1, C5, and C6-8 HPAIVs were found in specific areas. Two strains had phylogenetically distinct hemagglutinin (HA) and non-structural (NS) genes from other previously identified strains, respectively. The estimated latest divergence date between the viral genotypes suggests that genetic reassortment occurred in bird populations before their winter migration to Japan. Antigenic differences in 2016-2017 HPAIVs were not observed, suggesting that antibody pressure in the birds did not contribute to the selection of HPAIV genotypes. In the late phase, the majority of HPAI cases in wild birds occurred south of the lake freezing line. At the end of the outbreak, HPAI re-occurred in East coast region, which may be due to the spring migration route of Anas bird species. These trends were similar to those observed in the 2010-2011 outbreaks, suggesting there is a typical pattern of seeding and dissemination of HPAIV in Japan.

After an H5N1 highly pathogenic avian influenza virus (HPAIV) caused a highly pathogenic avian influenza (HPAI) outbreak on a goose farm in China in 1996, its virus progenies (Gs/Gd-like viruses) have experienced antigenic mutations and genetic reassortment [42]. Gs/Gd-like viruses are classified into 10 clades ranging from 0 to 9 according to the phylogenetic lineage of H5 subtype HA [40]. Some of these viruses have spread to Asia, Europe, North America, and Africa, causing extensive damage to the poultry industry [1]. Moreover, the H5 HPAIV returned to the natural host, migratory water birds [4, 22], increasing the possibility of continuous spread of virus worldwide. More recently, clade 2.3.4.4 HPAIVs have become common, subdivided into eight groups, a to h [41].

Japan has experienced outbreaks of HPAI every few years. Of these, the outbreaks in the winters of 2010–2011 and 2016–2017 occurred across an unusually wide area. In the 2010–2011 season, clade 2.3.2.1 H5N1 HPAIVs were isolated from 63 wild birds including migrating and resident birds, and caused outbreaks in 24 chicken farms [29, 38]. The causal viruses were classified into three subgroups according to their hemagglutinin (HA) gene lineage, suggesting that virus dissemination was associated with multiple flyways of migratory waterbirds around Japan [34]. In 2016–2017, 218 wild/captive birds were revealed to be infected with clade 2.3.4.4e HPAIVs, which caused further outbreaks in 12 chicken/duck farms. Genotypes C1–C8 (excepting C3, reported in Korea [18]) were identified from the causal viruses [26, 35]. The first report regarding the genetic information of the causal viruses in the outbreak in Japan was co-released by several institutes conducting HPAI diagnosis on commission: Hokkaido University, National Institute of Animal Health, Kyoto Sangyo University, Tottori University, and Kagoshima University [25]. Subsequently, each institute reported specific cases; Takemae et al. classified the genetic background of the causal HPAIVs in poultry farms and Kyoto Racecourse [35]. Hiono et al. reported on characteristics of HPAIVs isolated from wild/captive birds in Northern Japan [9]. Successive HPAI outbreaks in wild birds in limited areas (Kagoshima and Ibaraki Prefectures) were described by Ozawa et al. and Tsunekuni et al., respectively [26, 36]. Usui et al. investigated the dynamics of virus dissemination in captive birds in zoos [39]. While these reports have collated important information from each specific HPAI situation, a comprehensive study including geographical and time course validations is required to maximize the understanding of the 2016–2017 experience, to better respond to expected future cases. HPAI occurred in Japan in 2017–2018 [23, 37] and 2020–2021 [12, 15, 30] by clade 2.3.4.4b H5N6 and H5N8 viruses, respectively, and is likely to continue. In the present study, 58 HPAIVs isolated at Tottori University add to the phylogenetic analyses together with other reported strains to understand overall HPAI circumstances in 2016–17 in Japan. Additionally, reference strains were selected based on phylogeny, and antigenic relationships were elucidated.

MATERIALS AND METHODS

Isolation and identification of viruses

Fifty-eight viruses that were examined in the present study are shown in Supplementary Tables 1 and 2 (listed as “present study” in the references column). These viruses were obtained following definitive diagnosis and ex-post epidemiological investigation conducted by Tottori University in collaboration with Ministry of Environment and Ministry of Agriculture in the winter of 2016–2017. For virus isolation, tracheal/cloacal swabs of the dead birds or environmental samples on the farms were collected in 2 ml of phosphate buffered saline. The samples were passed through 0.45 µm pored filter (Sartorius, Göttingen, Germany) and inoculated into the allantoic cavities of 10-day-old embryonated chicken eggs (Aoki breeder farm, Tochigi, Japan). After incubation for 48 hr at 35°C, the allantoic fluids were harvested and their HA activities were examined [10]. Subtypes of the HA and neuraminidase (NA) of influenza virus isolates were identified by hemagglutinin inhibition (HI) test [31] and neuraminidase inhibition test [2] using the reference antisera of avian influenza viruses. The viruses were stored as the seed virus stock and used in subsequent experiments.

Sequencing

Viral RNA was extracted from chicken embryos infected with viruses using QIAamp Viral RNA Mini kit (Qiagen, Hilden, Germany) and reverse transcribed with the Uni12 primer [11] and PrimeScript reverse transcriptase (Takara Bio Inc., Kusatsu, Japan). The full length of eight gene segments was amplified by polymerase chain reaction with universal [11] and/or gene-specific primers. Direct sequencing of each gene segment was performed using BigDye Terminator version 3.1 Cycle Sequencing Kit with 3130 Genetic analyzer (Applied Biosystems, Waltham, MA, USA). For sequencing some isolates, next-generation sequencing was applied as follows. MiSeq libraries were prepared using KAPA RNA Hyper Prep Kit (Illumina, Inc., San Diego, CA, USA) and KAPA Dual-Indexed Adapter Kit (Roche, Basel, Switzerland). The prepared MiSeq libraries were sequenced on a MiSeq by using MiSeq Reagent kit v3 (Illumina, Inc.) with 2 × 300 bp paired-end read length. The sequencing data were analyzed using GENETYX Network version 12 (Genetyx Co., Tokyo, Japan), GeneStudio (https://genestudio.com/), or CLC Genomics Workbench 12 (Qiagen). The gene sequences obtained in the present study have been registered at the DNA Data Bank of Japan and were already shared by the National Center for Biotechnology Information and GISAID databank (Supplementary Tables 1 and 2).

Phylogenetic analysis

Combined eight-gene segments, each viral gene segment, and the HA protein sequences of the virus strains, whose full genotypes were determined previously or in the present study, were aligned using the MUSCLE program [6], and the best-fit substitution model was selected using MEGA X: Molecular Evolutionary Genetics Analysis across computing platforms [17]. Genetic data of some HPAIVs isolated by the other institutes could not be obtained (designated as “unknown” in this report). The confirmed models of nucleotide substitution are shown in Supplementary Table 3. Maximum likelihood (ML) trees were constructed using the MEGA X software, with a resampling process of 1,000 replicates. A maximum clade credibility (MCC) tree of the combined-eight gene segments was constructed according to the Bayesian Markov chain Monte Carlo (MCMC) method using the Bayesian Evolutionary Analysis by Sampling Trees (BEAST) 2 software package version 2.6.6 [3]. A strict clock model was applied to estimate the substitution rates. MCMC chains were run for 1 × 109 steps in length, where the number of steps was determined to obtain an explained sum of squares of >200 (as assessed by Tracer version 1.7.2 [27]). The MCC tree was generated from the MCMC samples using the TreeAnnotator software program, version 2.6.6, of the BEAST 2 package. The constructed tree was exported to FigTree v1.4.4. (http://tree.bio.ed.ac.uk/software/figtree/) and edited.

Antigenic analysis

Representative viruses were selected from genotypes C1–C6 except C3 (Table 1). These viruses together with two Niigata strains with specific HA or NS genes (Supplementary Table 4) and three Japanese clade 2.3.4.4 HPAIV isolates in other seasons, A/mandarin duck/Gifu/2112D001/2014 (H5N8), A/mute swan/Shimane/3211A001/2017 (H5N6), and A/water/Tottori/NK1209/2020 (H5N8) were used to prepare a panel of chicken polyclonal antisera. Chickens were immunized with the viruses that had been inactivated with formalin (FUJIFILM Wako Pure Chemical Corp., Tokyo, Japan; final concentration, 0.1%) and purified by sucrose density-gradient ultracentrifugation [16]. Antigenic characterization of the clade 2.3.4.4 H5 AIVs and an additional 6 strains which were assigned to significantly independent groups based on the phylogenic trees of their HA genes and/or proteins (Supplementary Fig. 1D, 1I) was performed by cross-HI test using the polyclonal antisera and formalin-inactivated virus with glycerol (stabilizer, FUJIFILM Wako Pure Chemical Corp.) and sodium azide (antiseptic, Nacalai Tesque, Kyoto, Japan) at final concentrations of 50% and 0.08%, respectively.

Table 1.

Cross-reactivity of Japanese clade 2.3.4.4 H5 highly pathogenic avian influenza virus isolates in hemagglutination inhibition test

Strain	HA clade	Gene conste- llation	Chicken anti-serum against
			2014–2015	2016–2017							2017–2018	2020–2021
			2.3.4.4c	2.3.4.4e							2.3.4.4b	2.3.4.4b
				C1	C2	C4	C5	C6	1511C003	5112007
2014–2015 HPAIV
A/mandarin duck/Gifu/2112D001/2014 (H5N8)	2.3.4.4c		512*	16**	16**	16**	32**	64**	16**	128**	8**	64
2016–2017 HPAIV
A/coot/Shiga/2501T010/2017 (H5N6)	2.3.4.4e	C1	128	512	1,024	512	1,024	1,024	1,024	2,048	64**	128
A/teal/Tottori/1/2016 (H5N6)	2.3.4.4e	C2	16**	128	512	128	512	512	512	1,024	16**	32**
A/Tundra swan/Niigata/5112004/2016 (H5N6)	2.3.4.4e	C2	16**	128	512	1,024	512	512	512	512	32**	64
A/environment/Aomori/4/2016 (H5N6)	2.3.4.4e	C2	8**	128	256	512	256	256	256	256	32**	64
A/snowy owl/Akita/0051D007/2016 (H5N6)	2.3.4.4e	C2	4**	128	512	1,024	512	512	512	512	32**	32**
A/environment/Saga/4/2017 (H5N6)	2.3.4.4e	C4	32**	128	256	128	256	256	256	1,024	8**	16**
A/chicken/Kumamoto/45/2016 (H5N6)	2.3.4.4e	C4	16**	128	512	512	256	256	512	512	32**	64
A/teal/Tottori/2/2016 (H5N6)	2.3.4.4e	C5	128	512	1,024	512	512	512	512	2,048	32**	64
A/duck/Tottori/E9/2016 (H5N6)	2.3.4.4e	C5	16**	128	256	512	512	256	256	512	32**	32**
A/mute swan/Hyogo/2801ITM015/2017 (H5N6)	2.3.4.4e	C5	16**	128	512	512	256	256	256	512	32**	64
A/common porchard/Yamaguchi/3501B002/2017 (H5N6)	2.3.4.4e	C6	32**	128	256	128	512	512	512	1,024	32**	32**
A/Tundra swan/Niigata/1511C003/2016 (H5N6)	2.3.4.4e	Other	64**	256	512	256	1,024	1,024	1,024	1,024	64**	32**
A/Tundra swan/Niigata/5112007/2016 (H5N6)	2.3.4.4e	Other	32**	128	256	128	256	512	512	1,024	32**	32**
2017–2018 HPAIV
A/mute swan/Shimane/3211A001/2017 (H5N6)	2.3.4.4b		128	64**	256	64	512	1,024	256	512	1,024	256
2020–2021 HPAIV
A/water/Tottori/NK1209/2020 (H5N8)	2.3.4.4b		128	16**	128	128	128	256	64**	128**	512	256

* Homologous titer. ** The titer eight times and less lower than homologous titer.

RESULTS

Phylogenetic analysis of the Japanese H5N6 HPAIVs in 2016–2017

Full genome sequences of 58 Japanese H5N6 HPAIVs isolated in 2016–17 were determined for the first time (Supplementary Tables 1 and 2). Together with the previously published 130 strains with full genome, an ML tree was constructed using combined eight gene segments of each strain (Fig. 1). The majority of isolates were classified into groups C2 and C5. Group C7 and C8 only contained the strains from the Kagoshima Prefecture reported by Ozawa et al [26]. Group C2 comprised of strains from the reported clusters in the Ibaraki Prefecture (groups 1 and 2) [36] or zoos (Omoriyama and Higashiyama) [39]. Two strains from group C1 which had minor NS gene segments (NS-II; Supplementary Fig. 1H and Supplementary Table 4) were included in the C2 cluster (Fig. 1). Similarly, two strains isolated in the Niigata Prefecture, A/Tundra swan/Niigata/1511C003/2016 (H5N6) and A/Tundra swan/Niigata/5112007/2016 (H5N6) had phylogenetically distinct HA and NS genes from the other C2 strains (Fig. 1, Supplementary Fig. 1D, 1H, and Supplementary Table 4; designated as “other” in the present study). HA genes of the 1511C003 strain showed the highest homology (97%, 56 nucleotide differences) with the #5112007 strain under web BLAST search ( https://blast.ncbi.nlm.nih.gov/Blast.cgi). NS genes of the #5112007 strain was most similar to Indonesian H5N1 HPAIVs (data not shown). A/common pochard/Yamaguchi/3501B002/2017 (H5N6) was additionally classified into the group C6 together with the reported isolates from the Kagoshima [26] and Miyazaki Prefectures [35] (Fig. 1). Group C5 included many isolates from the Hyogo Prefecture together with mute swan isolates from Kyoto Racecourse as reported previously [35]; 15 strains isolated from mute swans in Koya Pond, Itami city, Hyogo (ITM strains) were phylogenetically similar with the other strains in the same prefecture.

Fig. 1.

Phylogenetic tree of the combined sequences of eight gene segments of H5N6 highly pathogenic avian influenza viruses isolated in the winter of 2016–2017 in Japan. Horizontal distances are proportional to the minimum number of nucleotide differences required to join nodes and sequences. Numbers at the nodes indicate confidence levels in a bootstrap analysis with 1,000 replications. Bootstrap values of 70% or more are shown at each branch. Each strain name is colored by isolated region as shown in the left side map. Some previously proposed clusters (Groups 1 and 2 in Ibaraki Prefecture [36], zoos [39], Kyoto Racecourse [35], and C7 [26]) were compressed.

Bayesian analysis provided a similar classification of the Japanese H5N6 HPAIV 2016–2017 winter strains as the ML method (Fig. 1 and Supplementary Fig. 2). According to the MCC tree, groups C2 and C4–C8 were estimated to diverge genetically before 2016 (Supplementary Fig. 2A). Furthermore, groups C1 and C other were estimated to have diverged from the C2 cluster on April 7, 2016 (95% highest posterior density [HPD]: February 18 to June 10, 2016) and July 20, 2016 (95% HPD: June 12 to August 24, 2016), respectively (Supplementary Fig. 2B).

Epidemiology

The locations and times of HPAI outbreaks in wild birds and poultry are shown in Fig. 2. Group C2 HPAIVs were maintained all over Japan from the beginning to the end of the outbreak; the C2 strain outbreaks were observed in the North-Eastern region in the early stage (Supplementary Fig. 3A), and then moved to the Central region (Supplementary Fig. 3B). Group C5 strains were isolated until January 2017 predominantly in specific prefectures (Aichi, Kyoto, Hyogo, and Tottori) in the Central region (Fig. 2). Similarly, group C1, C4, and C6 HPAIVs were found in limited areas; the C1 strains were isolated in adjacent prefectures, Gifu and Shiga, and isolation sites of the C4 and C6 strains were concentrated in the South-Western region. Multiple group HPAIVs were isolated in singular locations at the Miyagi, Niigata, and Tottori Prefectures as well as the Kagoshima Prefecture reported by Ozawa et al [26]. Two HPAIVs with phylogenetically distinct HA or NS genes, as mentioned above, were found together with C2 strains in migratory mute swan flocks in Hyoko Lake, Niigata Prefecture (Fig. 2). Group C2 and C5 strains were isolated from teal feces on the same day at Nikko pond, Tottori Prefecture (Fig. 2 and Supplementary Table 1). C4 and C5 strains were found sporadically in the Miyagi Prefecture (Fig. 2).

Fig. 2.

Map of H5N6 virus infection in wild/captive birds (circles) and poultry farms (triangles) in the winter of 2016–2017 in Japan. Each case is represented by different colors according to the times of outbreak and genetic background of the causal viruses.

Antigenic relationships among Japanese clade 2.3.4.4 HPAIVs

Based on the genetic classification of Japanese HPAIVs in 2016–2017, reference strains for each group were selected from stocked viruses at Tottori University (Table 1). Chicken antisera against these strains were prepared for cross-HI tests to evaluate antigenic relationships among Japanese clade 2.3.4.4 HPAIVs (Table 1). Clade 2.3.4.4e HPAIVs in 2016–2017 were antigenically different from clade 2.3.4.4c and 2.3.4.4b HPAIVs in other seasons; typically, antisera against clade 2.3.4.4c and 2.3.4.4b HPAIVs showed more than eight-fold reduced reactivity to clade 2.3.4.4e HPAIVs. Conversely, no considerable antigenic differences were observed among clade 2.3.4.4e HPAIVs.

DISCUSSION

The present study clarified the genetic background, epidemiology, and antigenicity of Japanese HPAIVs isolated in the winter of 2016–2017. Locations where HPAIVs were found in the early phase of the outbreak (Supplementary Fig. 2A) were clustered into three regions: North-Eastern, Central, and South-Western, as observed in the 2010–2011 outbreak [34]. These trends are likely related with multiple winter migration routes of birds from the Asian continent to Japan, as reviewed previously [43]. Similar trends were observed in subsequent HPAI outbreaks; in the outbreaks of 2017–2018 [13] and 2020–2021 [12], the first HPAIV isolates were found in the Shimane (Central region) and Hokkaido (North-Eastern), respectively. Multiple genotype HPAIVs were isolated in the Miyagi, Niigata, Tottori, and Kagoshima Prefectures (Fig. 2), which have major wintering sites for Anseriformes and Gruiformes [8]. No antigenic differences were observed among these HPAIV isolates (Table 1). Bayesian analysis suggested that genetic reassortment among the viruses circulating at that time had already occurred before the 2016–2017 winter migration season (Supplementary Fig. 3). These data suggest that the HPAIVs diverged independently of antibody pressure in the birds, and that the strains were individually transmitted in Japan. However, the possibility that further viral genetic reassortment occurred during the 2016–2017 winter season in Japan cannot be ruled out; we actually found a dead mallard (Anas platyrhynchos) that had contracted multiple HPAIVs during the 2020–2021 outbreak (unpublished data).

In the late phase of the 2016–2017 outbreak, most of the HPAI cases in wild birds occurred south of the lake freezing line [7] (Supplementary Fig. 3B). Significant HPAI outbreaks at Koya Pond, Hyogo Prefecture in mute swans (Supplementary Table 1) and jungle crows [33] occurred in the 2016–2017 and 2017–2018 winters, respectively, indicating that specific sites have a continuous risk of HPAIV invasion and perpetuation. Genotype C1 and C6 HPAIVs were isolated in limited areas (Fig. 2). Of these, C1 HPAIVs were found in both wild bird and poultry farm, supporting the notion that wild birds play a key role in direct/indirect transmissions of HPAIV to poultry farms. Interestingly, genotype C4 HPAIV was isolated in the Miyagi Prefecture on 23rd March, far from the South-Western region where the C4 virus caused HPAI previously (Fig. 2). On the same day, HPAI also occurred in the Chiba Prefecture, and these were the last cases of the 2016–2017 winter (Supplementary Table 2). The final HPAI case of the 2010–2011 winter was also reported in Chiba Prefecture [29]. The reason why HPAI have usually occurred in the East coast region in final phase of winter season is likely due to spring migration routes of common Anas; many stopover sites of mallards (Anas platyrhynchos) during spring migration are distributed in the North-Eastern area [44]; some Eurasian wigeons (Mareca penelope) moved inside of Japan in a northerly direction from the South-Western region in spring [5]. Additionally, north of the freezing line, a limited number of HPAIVs were isolated mainly from whooper swans (Supplementary Fig. 3B and Supplementary Table 1). The North-Eastern region has been a major wintering site for swan species [32]; therefore, HPAIVs might have been maintained in these areas until beginning of spring migration, around late February.

Fortunately, immediate definitive diagnoses of HPAI in 2016–2017 were performed by each institute involved, however it caused massive economic losses in the poultry industry [24]. Clade 2.3.4.4b HPAIVs in subsequent seasons were antigenically distinct from clade 2.3.4.4e HPAIVs in 2016–2017 (Table 1). Clade 2.3.4.4 HPAIVs have genetically and antigenically evolved [41], so that invasion of a novel phenotype of HPAIV remains a concern. Fast detection of HPAIV in early winter and sharing and disclosure of viral information is essential for rapid and correct diagnoses. Based on such a policy, Tottori University and Hokkaido University published a press release or article regarding the first isolated HPAIVs in 2017–2018 and 2020–2021 winters, respectively [12, 13]. Further, HPAI have usually been reported in neighboring countries such as Korea prior to occurrence of outbreaks in Japan [14, 18,19,20]. Consequently, an ongoing HPAI emergency action level announcing system by the Ministry of Environment in Japan is important. This should play a key role to assist in rapid and proper prevention countermeasures hereafter, and criteria for setting the action level should be continuously discussed and updated according to HPAI conditions in the future.

The present study followed up the dynamics of HPAIV invasion and dissemination by migratory birds in Japan. According to recent climate change estimates, migratory flyway and wintering period of migratory birds are likely to change [21, 28]. In order to eradicate future HPAI in Japan, comprehensive studies of HPAIV dynamics should be continued.

POTENTIAL CONFLICTS OF INTEREST

The authors declare no conflicts of interest associated with this manuscript.