Genomic evolution, transmission dynamics, and pathogenicity of avian influenza A (H5N8) viruses emerging in China, 2020.

Multiple recent outbreaks of highly pathogenic H5N8 viruses originating in aquatic birds frequently occurred in most European countries, Russia, South Korea, and Japan during the winter of 2020-21, and one zoonotic event of poultry workers infected with novel H5N8 viruses were reported in Russia. Strikingly, these novel H5N8 viruses had emerged and been co-circulating in wild birds and poultry in multiple provinces of China during 2020-21. In China, the population of aquatic birds has risen significantly in the past twenty years, and China is regarded as the largest reservoir for influenza viruses carried in aquatic birds across the globe. Hence, the co-circulation of these novel H5N8 viruses poses an alarming threat to not only poultry industry but also human health. In this study, we sequenced full-length genomes of these H5N8 viruses circulating in China. Phylogenetic analysis demonstrated that poultry-origin H5N8 viruses in China fell within wild birds-origin clade 2.3.4.4b H5N8 viruses from Europe during 2020-21, and notably, were genetically closely related to human-infecting H5N8 viruses in Russia. Moreover, they possessed several molecular markers associated with mammalian adaption. Bayesian coalescent analysis showed that these H5N8 viruses might have introduced into China during June-September 2020, suggesting that these H5N8 viruses might have introduced via wild bird migration or poultry trade. Besides, we also found that the effective population size of clade 2.3.4.4b H5N8 viruses dramatically increased during the winter season of 2020/21, as is consistent with previous increase of genetic diversity during the winter seasons of 2013/14 and 2016/17, which indicated that the wild bird migration accelerates the genetic diversity of these H5N8 viruses during the winter season of 2020/21. Notably, these novel H5N8 viruses were lethal to chickens and mice, highly transmissible to ducks, and were antigenically distinct from 2.3.4.4h H5 viruses circulating in China, posing considerable threats to public health. Our findings offer novel insights into the evolution and risk assessment of H5N8 viruses during the winter season of 2020-21.

1. Introduction

The Gs/GD lineage of H5N1 viruses that has evolved rapidly and disseminated across multiple continents (Xu et al. 1999), currently poses great threats to the poultry industry and human health. The global co-circulation of the H5 subtype viruses has evolved in multiple distinct subclades, among which clade 2.3.4.4 has become dominant (Zhang et al. 2020a). Recently, the nomenclature hemagglutinin (HA) of the 2.3.4.4 H5 was divided into subclades 2.3.4.4a to 2.3.4.4h according to the World Health Organization (Smith et al. 2015; Yamaji et al. 2020).

In early 2014, outbreaks of 2.3.4.4b H5N8 viruses were reported in South Korea (Lee et al. 2014), and, subsequently, novel 2.3.4.4b H5N8 viruses disseminated by ways of migratory birds across the globe to Europe, North America, and Africa beginning in late October 2016 (Global Consortium for H5N8 and Related Influenza Viruses 2016). During the winter of 2020–21, multiple outbreaks of H5N8 viruses derived from wild birds and poultry were frequently reported in most European countries, South Korea, and Japan. At the same time, from November 2020 to February 2021, these novel H5N8 viruses also triggered multiple outbreaks among wild swans in Beijing, Shandong, Shanxi, and Jiangsu provinces, China (Ministry of Agriculture and Rural Affairs of the People’s Republic of China 2021). More seriously, on 20 February 2021, one zoonotic event of seven poultry workers infected with novel H5N8 viruses was reported in Russia, and the possibility of human-to-human transmission cannot be ruled out (World Health Organization 2021).

In China, the population of domestic ducks has risen significantly in the past twenty years. As the site where a whopping seventy-five per cent of the aquatic birds are bred, China is regarded as the largest reservoir for influenza viruses carried in aquatic birds across the globe (Van Boeckel et al. 2011). In addition, live poultry markets (LPMs) in China often allow close contact between migratory birds, poultry, and humans (Venkatesh et al. 2018). Therefore, the co-circulation of these novel H5N8 viruses poses an alarming threat to not only poultry industry but also human health. Against that background, we sequenced and analyzed the full-length genomes of H5N8 viruses isolated from poultry in China during December 2020 and explored the antigenicity and pathogenicity of these H5N8 viruses in poultry and mammals. Our finding initially offers novel insights into the genomic evolution, transmission dynamics, and risk assessment of avian influenza A (H5N8) viruses during the winter season of 2020–21.

2. Materials and methods

2.1 Ethics statement and biosafety

All experiments with all available avian influenza A (H5N8) viruses were conducted in an animal biosafety level 3 laboratory and animal facility according to the protocols of South China Agricultural University (SCAU) (CNAS BL0011) protocols. All animals involved in experiments were reviewed and approved by the Institution Animal Care and Use Committee at SCAU and treated in accordance with its guidelines (2017A002). The animals were monitored twice per day for clinical signs, and we provided sufficient chicken feed and water for experimental animals every day.

2.2 Sample collection and virus isolation

From January to December 2020, we collected cloacal and tracheal swab samples from LPMs in southern, southwestern, and northern China. Each sample was placed in 2 ml of the PBS supplemented with penicillin (5,000 U/ml) and streptomycin (5,000 U/ml) and inoculated in the allantoic cavities of ten-day-old specific-pathogen-free (SPF) embryonated chicken eggs at 37 °C. The allantoic fluid was collected and tested via HA assay with one per cent chicken red blood cells before use.

2.3 RNA extraction, RT-PCR, and DNA sequencing

RNA was extracted from the suspension of viral isolates with the RNeasy Mini Kit (QIAGEN) as directed by the manufacturer. Two-step Reverse transcription-PCR (RT-PCR) was conducted with universal primers as previously described (Qi et al. 2018). Briefly, RNA was reverse transcribed into cDNA using the M-MLV reverse transcription (TakaRa) with Uni12 primer (5′-AGCAAAAGCAGG-3′) for 1 h at 42 °C. Then, eight full-length genome sequences of H5N8 viruses were amplified using PrimeSTAR Max DNA Polymerase (TakaRa) with frame-specific primers. PCR program was set as follows: initial denaturation at 95 °C for 5 min, followed by thirty cycles of 95 °C for 30 s, 55 °C for 30 s, and 72 °C for 2 min, and finally extension at 4 °C until use. PCR products were purified with a Gel Extraction Kit D2500 (OMEGA), and the gene segments were sequenced by TSINGKE Co., Ltd. (Guangdong, China).

2.4 Phylogenic analysis

All eight available genome sequences with the complete coding regions of avian influenza A (H5N8) viruses were downloaded from GISAID (http://www.gisaid.org/), after which datasets for the sequences together with new three H5N8 strains were aligned using MAFFT (version 7.149) program (Katoh et al. 2002). Maximum likelihood (ML) phylogenies for the codon alignment of the eight genome sequences were estimated using the GTRGAMMA nucleotide substitution model in the IQ-TREE 1.68 software (http://www.iqtree.org/release/v1.6.0). Node support was determined by nonparametric bootstrapping with 1,000 replicates, and the phylogenetic tree was visualized in the FigTree (version 1.4.3) program (http://tree.bio.ed.ac.uk/software/figtree/).

2.5 Bayesian maximum clade credibility phylogeny

We estimated rates of evolutionary change (i.e. nucleotide substitution) in the eight gene segments of avian influenza A (H5N8) viruses. To ensure their sufficient temporal structure in terms of alignment for reliable rate estimation, we performed a regression of root-to-tip genetic distances in the ML tree against exact sampling dates using the TempEst (Rambaut et al. 2016). To obtain a more robust estimates of the rates, we used the Bayesian Markov chain Monte Carlo (MCMC) method implemented in the BEAST package (version 1.8.2), employing the GTRGAMMA nucleotide substitution model, an uncorrected lognormal relaxed molecular clock model, and a Bayesian skyride coalescent model (Minin et al. 2008). Multiple runs of maximum clade credibility method were combined using LogCombiner (version 1.8.3), with one trillion total steps for each set and sampling every 10,000 steps. The convergence (i.e. effective sample sizes >200) of relevant parameters was assessed using Tracer (version 1.6) (http://beast.bio.ed.ac.uk/). The posterior distribution of trees obtained from BEAST analysis (with 10% of runs removed as burn-in) was also used to obtain the MCC tree for the eight gene segments. After that, we performed BEAST analysis of HA genes of clade 2.3.4.4b H5N8 viruses from 2014 to 2021 with the same parameters used to estimate relative genetic diversity. All phylogenetic trees were visualized in FigTree (version 1.4.3) program (http://tree.bio.ed.ac.uk/software/figtree/) with nodes in ascending order.

2.6 Phylogeographic inference

Time-measured phylogenies were inferred using Bayesian discrete phylogeographic approach implemented in BEAST package (version 1.8.2) (Drummond et al. 2012). The discrete sampling locations of clade 2.3.4.4b H5N8 viruses during 2020–21 include China (n = 4), Czech Republic (n = 9), England (n = 7), Germany (n = 7), Hungary (n = 14), Italy (n = 6), Japan (n = 4); Kazakhstan (n = 5), Netherlands (n = 4), Poland (n = 18), Russia (n = 9), and Sweden (n = 2). We first performed a regression of root-to-tip genetic distances on the ML tree against exact sampling dates using the TempEst (Rambaut et al. 2016), which showed strong temporal signal. Then, we used an UCLN relaxed molecular clock model (Drummond et al. 2006). Besides, a Bayesian stochastic search variable selection (BSSVS) model with asymmetric substitution was also used. For each independent dataset, multiple runs of MCMC method were combined using LogCombiner (version 1.8.3), utilizing 10,000,000,000 total steps for each set, with sampling every 10,000 steps. Subsequently, we used SpreaD3 v0.9.6 to develop interactive visualizations of the dispersal process through time and to compute a Bayes factors (BFs) test to assess the support for significant individual transitions between distinct geographic locations. SpreaD3 takes a rate matrix file for location states generated under the BEAST analysis using the BSSVS procedure (Bielejec et al. 2016). The BF values >100 indicated robust statistical support. BF values >30 and ≤100 indicated very strong statistical support. BF values >10 and ≤30 indicated strong statistical support. BF values >3 and ≤10 indicated substantial statistical support. BF values <3 indicated poor statistical support (Lemey et al. 2009). We used R (https://www.r-project.org) and ArcGIS Desktop 10.4 software (http://www.esri.com/software/arcgis/arcgis-for-desktop/) to create plots showing results of BF tests.

2.7 Antigenic analysis

Hemagglutinin inhibition (HI) antibody titers of the serum of H5 vaccinated chickens (n = 8) and H52001/H5N6 virus derived from clade 2.3.4.4h together with three 2.3.4.4b A/Chicken/China/HB181/2016(H5N8), A/Goose/Jiangsu/JS1603/2014(H5N8), A/duck/southwestern China/B1904/2020(H5N8) viruses in China during 2014–20. In addition, serum from twenty-one-day-old-SPF chickens challenged by H5 subtype (rFJ56, H51901, H52001, H51904, H51906, H51907, H51908, 17598, 181203, and 20053) viruses was prepared from cross-HI assay. Both the antigen and antiserum were available from Guangzhou South China Biological Medicine. The HI test assay was a standard beta test, whereby 4 HA units of influenza viruses in ninety-six-well plates and the two-fold serially diluted serum prepared previously were added. The highest serum dilution that produced complete inhibition of HA activity. The cross-HI data were analyzed using antigenic cartography (http://www.antigenic-cartography.org), which is a method to visualize and increase the resolution of HI results. Each point on the map represents an HA protein antigen. The distance between two HA protein antigens on the map represents the antigenic distance between the two antigens. Points are colored based on categorical hierarchical clustering.

2.8 Animal experiment

In the first chicken experiment, the intravenous pathogenicity index (IVPI) test was performed to determine the pathogenicity of the H5N8 viruses according to the recommendations of the World Organization for Animal Health. In brief, six-week-old SPF chickens (Guangdong Dahuanong Animal Health Products Co., Ltd., Guangdong Province, China) were intravenously inoculated with 0.1 ml of a 1/10 dilution of fresh infectious allantoic fluid. All chickens were examined daily until died. At each observation, each chicken was scored on the basis of its condition, as follows: 0 for normal, 1 for sick, 2 for severely sick, and 3 for dead. IVPI was the mean score per chicken per observation over experimental period.

In the second chicken experiment, groups of ten six-week-old SPF chickens (Guangdong Dahuanong Animal Health Products Co., Ltd., Guangdong Province, China) were inoculated intranasally with 105 fifty per cent TCID50/200 μl of the B1904/H5N8 viruses. Chickens were observed for clinical symptoms until died. Three chickens were dissected to test for virus replication in the organs, including the heart, liver, spleen, lung, kidney, and brain. Tissue samples were collected for virus titration by an EID50 assay.

Groups of eleven Pekin ducks tested negative for H5, H7, and H9 subtype viruses were inoculated intranasally with the B1904 H5N8 viruses. Pekin ducks were observed for clinical symptoms for fourteen days. Five ducks were introduced into the isolators with infected ducks twenty-four hours after the latter ducks were directly inoculated to assess transmission. Cloacal and throat swab specimens of the ducks were collected at three, five, seven, and ten days post infection (dpi) from the infected and naive ducks. All the swabs samples were inoculated in the allantoic cavities of ten-day-old embryonated chicken egg at 37 °C. The allantoic fluid was collected and tested for HA assay with one per cent chicken red blood cells and then used in this study.

Groups of twelve six-week-old female BALB/c mice, obtained from the Vital River Company in Beijing, were anesthetized with isoflurane and inoculated intranasally with 106 TCID50/50 μl B1904 H5N8 viruses. Body weight and clinical signs were monitored daily for fourteen days after infection. To test for virus replication in the organs, three mice were euthanized at 3 and 5 dpi. Tissue samples including heart, liver, spleen, lung, and brain and turbinate samples from each euthanized mouse were collected for virus titration by an EID50 assay.

3. Results and discussion

3.1 Virus sequencing and molecular characterization

In December 2020–February 2021, avian influenza A (H5N8) viruses were detected in wild birds and poultry in northern and southwestern China (Fig. 1D). Considering that event, we isolated three H5N8 viruses, designated A/duck/Northern China/ZGL/2020(H5N8), A/duck/Southwestern China/B1904/2020 (H5N8), and A/duck/Northern China/LSP/2020 (H5N8), and sequenced full-length genomes. All of three H5N8 viruses possessed highly similar identities, and all the sequences in this study were submitted to the GISAID’s EpiFlu Database under accession nos. EPI1844088-EPI1844094, and EPI1844096–EPI1844113. BLAST analysis showed that full-length genomes of H5N8 viruses were closely related to England, Czech Republic, Russia, and Kazakhstan, with the homology of 99.5%–99.8%. All of the H5N8 viruses in this study possessed multibasic cleavage site motif of REKRRKR/GLF motif in HA protein, suggestive of highly pathogenic in chickens (Chen et al. 1998). In addition, these H5N8 viruses contained 226Q in the HA protein (Hay et al. 1986), indicating that the viruses preferentially bind to avian-like influenza receptors, and none of the viruses bear the E627K, A588V, and D701N mammalian adaptation substitutions at the PB2 protein (Bi et al. 2015; Xiao et al. 2016). However, some other mutations in HA, PB1, PA, and M1 proteins associated with increased virulence in mice and transmission in guinea pigs were found in these H5N8 viruses (Table 1), posing increasing threat to human health.

Figure 1.

Phylogeny of avian influenza A(H5N8) viruses in China, 2020. (A) ML tree of HA gene of H5 subtype influenza viruses. All branch lengths are scaled according to the numbers of substitutions per site (subs/site). The tree is rooted using clade 2.3.2 A/duck/Vietnam/LBM360c1-4-1/2013(H5N6), which was collected on 6 February 2013. Red stars indicate newly sequenced viruses isolated from poultry in China. (B) Time-scaled evolution of HA gene sequences of clade 2.3.4.4b H5N8 viruses during 2016–20. A maximum clade credibility tree of the HA gene sequence of H5N8 viruses is shown, H5N8 viruses from different countries are denoted by different colors. Shaded bars represent the 95% highest probability distribution for the age of each node. (C) HA gene tree revealing a single cluster (cluster A/2020) of avian influenza A(H5N8) viruses during 2020. (D) The distribution of 2.3.4.4b avian influenza A(H5N8) viruses in China during 2020/21 winter season. The background indicated the sampling spaces of avian influenza A(H5N8) viruses during 2020/21 winter season in wild birds (red). The map was designed by using ArcGIS Desktop 10.4 software (ESRI, http://www.esri.com).

Table 1.

Key amino acid substitutions of A/duck/Northern China/ZGL/2020(H5N8), A/duck/Southwestern China/B1904/2020(H5N8), and A/duck/Northern China/LSP/2020(H5N8) viruses in China, December 2020.

Protein	Amino acid/motif	Phenotypic consequences
Hemagglutinin	226Q	Preferentially bind to avian-like influenza receptors
	156A	Increased virus binding to α2,6 and increased transmission in guinea pigs
	159P	Antigenic escape mutation
	cleavage site motif -REKRRKR/GLF	High pathogenicity of avian influenza viruses in chickens
Matrix protein 1	30D	Increased virulence in mice
	215A	Increased virulence in mice
Nonstructural protein 1	42S	Increased virulence in mice
Polymerase acidic protein	515T	Increased polymerase activity in mammalian cells
Polymerase acidic protein 2	588A	Associated with avian adaption
	627E
	701D
	526K

3.2 Phylogenetic analysis

In the ML trees, all of our H5N8 isolates were derived from the H5 clade 2.3.4.4b (Fig. 1A), and we found the full-length genome sequences of H5N8 viruses during 2020 mainly formed two independent clusters (cluster A/2020 and cluster B/2020) (Fig. 1B, Supplementary Figs S1–S8). The H5N8 viruses primarily circulating in Hungary, Poland, and Germany belonged to cluster B/2020, while H5N8 viruses in cluster A/2020 were derived from most European countries, China, Kazakhstan, South Korea, Japan, and Russia. The HA genes of our H5N8 isolates were genetically closely related to the A/goose/Russia_Novosibirsk_region/1-12/2020 (H5N8) virus. However, the PB2, NP, NA, M, and NS genes in this study clustered together in cluster A/2020. The time of the most recent common ancestors (tMRCA) of these genes were A/barnacle goose/Sweden/SVA201117SZ0468/KN003355/2020(H5N8) and A/greylag goose/England/033100/2020(H5N8) strains (Fig. 1C), indicating that these H5N8 viruses in China might had been undergoing reassortment with H5N8 viruses from European and Asian countries. Estimating the times of origin of H5N8 viruses during 2020 in two clusters revealed that the tMRCA in cluster B/2020 was much earlier than that of cluster A/2020 (Table 2); however, no significant differences emerged in the substitution rates of the eight gene segments between clusters A/2020 and cluster B/2020 (Supplementary Table S1). It is noteworthy that these H5N8 viruses in China and human-infecting H5N8 virus in Russia grouped together in cluster A/2020, which suggested that these H5N8 viruses in China have the potential threat to infect humans.

Table 2.

tMRCA estimates of full-length genomes of H5N8 viruses during 2020 in different clusters.

Segment	Cluster A/2020	Cluster B/2020	China-cluster/2020
PB2	October 2019 [May 2019, February 2020]	June 2019 [May 2019, September 2019]	June 2020 [May 2020, August 2020]
PB1	January 2020 [October 2019, March 2020]	June 2019 [March 2019, September 2019]	June 2020 [May 2020, July 2020]
PA	January 2020 [September 2019, February 2020]	July 2019 [May 2019, September 2019]	June 2020 [April 2020, July 2020]
HA	February 2020 [November 2019, March 2020]	July 2019 [March 2019, October 2019]	August 2020 [July 2020, September 2020]
NP	October 2019 [April 2019, March 2020]	July 2019 [May 2019, October 2019]	September 2020 [August 2020, November 2020]
NA	September 2019 [March 2019, February 2020]	April 2019 [January 2019, August 2019]	August 2020 [June 2020, October 2020]
M	February 2020 [November 2019, April 2020]	March 2019 [January 2019, May 2019]	August 2020 [June 2020, September 2020]
NS	October 2019 [May 2019, March 2020]	May 2019 [January 2019, August 2019]	August 2020 [June 2020, October 2020]

The estimation of the tMRCAs of eight gene segments indicated that poultry-origin H5N8 isolates in China might have entered China during June–September 2020 (Table 2), which strongly suggests that those H5N8 viruses circulating in China likely originated from Europe or the neighbouring Asian countries via wild bird migration or poultry trade. The migration of H5N8 virus-infected wild birds acts as a route of the long-distance spread of H5N8 viruses into countries with bird populations, and the frequent contact with wild birds remain the most probable cause of viral introduction into domestic poultry (Bouwstra et al. 2015; Hanna et al. 2015; Zhang et al. 2020a). Countries in the Middle East, Europe, and East Asia countries play a particularly important role in the continuing reassortment and global dissemination of avian influenza A (H5N8) viruses, because they lie along major migratory bird pathways–the central Asian, east Africa–west Asian, and Black Sea–Mediterranean routes (Zhang et al. 2020a). The close contact between aquatic birds carrying H5N8 viruses and some mammals including dogs and pigs in LPMs and poultry farms could further accelerate the evolution of these H5N8 viruses (Su et al. 2014; Su et al. 2015; Su et al. 2017), acquiring some mammalian adaption mutations or reassortment. Therefore, prevention and control efforts should be implemented in LPMs and poultry farms to improve their biosecurity measures, including the so-called 1110 strategy for LPMs first implemented in Guangdong Province (http://www.chinanews.com/fz/2014/12-06/6851778.shtml), as a mean to further reduce the spread of H5N8 influenza viruses.

3.3 Evolutionary and transmission dynamics

To estimate the population dynamics of clade 2.3.4.4b H5N8 viruses during 2014–21, we inferred the HA genes of H5N8 virus demographic history using the Bayesian Skyride plots. Among the results, following a sharp increase in winter 2013–14, genetic diversity stabilized during 2014–15 and then declined during 2015–16 (Fig. 2). A second increase was observed in winter 2016–17, when an outbreak occurred in bird population in Europe, after which the genetic diversity plateaued from mid-2017 to mid-2019 (Fig. 2; Supplementary Fig. S9). After that, genetic diversity declined sharply. However, a third increase of genetic diversity did arise in winter 2020–21 (Fig. 2). During the 2014/15 H5 virus epidemic, only few of the reassortment events were observed; however, frequent reassortment of H5 subtype viruses during 2016/17 and 2020/21 epidemic was found in wild bird populations (Global Consortium for H5N8 and Related Influenza Viruses 2016; Lycett et al. 2020; Lewis et al. 2021; Swieton et al. 2017). Combined with our population dynamics analysis, these evidences suggested that the recent outbreaks of H5N8 viruses originating in wild birds and poultry across the globe during 2020/21 winter season accelerated the genetic diversity of H5N8 viruses.

Figure 2.

Bayesian Skyride plot of HA genes of clade 2.3.4.4b H5N8 viruses from poultry and wild birds from 2014 to 2021. A Bayesian Skyride analysis of HA gene of clade 2.3.4.4b H5N8 viruses to display changes in the effective population size over time. The solid red line indicates the median value, and the shaded red area represents the 95% highest posterior density of genetic diversity estimates.

We then estimated the dissemination pathways of clade 2.3.4.4b H5N8 viruses during 2020–2021 followed a phylogeographic approach. Transmission routes with BF values >3 were selected for analysis. There were fifteen significant transmission routes of clade 2.3.4.4b H5N8 viruses. We found that Russia acted with the epicenter for the spread of H5N8 viruses during the winter season of 2020–21. Russia was linked with four locations, including Kasakhstan (BF = 116), England (BF = 24), Italy (BF = 10), and China (BF = 9) (Fig. 3; Table 3), which suggested that the outbreaks of H5N8 viruses in Europe might originated from Russia. In addition, we found that the higher migration rate (1.82) occurred between Russia and Kasakhstan (Fig. 3; Table 3), indicating that these novel H5N8 viruses had been co-circulating locally. Remarkably, we also found that the co-circulation of H5N8 viruses occurred in most European countries, including from England to Netherlands (BF = 33), from Italy to Poland (BF = 12), from Italy to Sweden (BF = 14), from Poland to Czech Republic (BF = 1933), from Poland to Germany (BF = 1186), from Poland to Hungary (BF = 25), and from Sweden to Poland (BF = 10) (Fig. 3; Table 3), for results showing that higher viral migration rates between closer European countries. In addition, we found that China is the main input location of the spread of H5N8 viruses, which was linked with closer countries, including Kasakhstan and Russia. Besides, a strong statistically support (BF = 175) of viral migration of H5N8 viruses was observed from China to Japan (Fig. 3; Table 3), suggesting that the role of migratory birds accelerated the dissemination of these novel H5N8 viruses. However, there are some limitations that some factors such as sampling and sequencing bias could have potentially affected the conclusion.

Figure 3.

Spatiotemporal dissemination of clade 2.3.4.4b H5N8 viruses from 2020 to 2021, which was determined by Bayesian phylogeography inference of HA gene sequences. Curves show the among-province virus lineage transitions statistically supported with BF >3 for H5N8 viruses. Curve widths represent transition rate values; curve colors represent corresponding statistical support (BF value) for each transition rate.

Table 3.

Statistically supported migration rates of clade 2.3.4.4b H5N8 influenza viruses during 2020–21 estimated from HA gene sequences.

From	To	BF	Migration rate
China	Japan	175	0.94
England	Netherlands	33	1.58
England	Poland	11	0.95
Germany	Japan	14	0.89
Italy	Poland	12	1.44
Italy	Sweden	14	1.10
Kasakhstan	China	8	0.96
Russia	China	9	0.82
Poland	Czech	1933	1.73
Russia	England	24	1.12
Poland	Germany	1186	1.18
Poland	Hungary	25	0.86
Russia	Italy	10	1.01
Russia	Kasakhstan	116	1.81
Sweden	Poland	10	1.32

The spring and autumn waves of H5N8 virus in Europe might derive from different sources, with wild birds coming from different directions. Previous study showed that the spread of H5 viruses along the autumn migration routes from Russia and Kazakhstan to the Black Sea region (Gilbert et al. 2006), and the close contact of poultry and wild birds in Romania accelerated the outbreak of H5 viruses during autumn (Bourouiba et al. 2010). In addition, multiple avian species migrated north in spring and gathered together in the northern breeding sites, and then migrated south to the wintering sites during autumn (Liang et al. 2010; Tian et al. 2015; Xu et al. 2016). We hypothesized that these H5N8 viruses circulating in China might originate from Europe or Russia along the multiple migratory flyways via southward autumn migration.

3.4 Pathogenicity of H5N8 viruses in chickens, ducks, and mice

The IVPI test in SPF chickens obtained for the A/duck/Southwestern China/B1904/2020(H5N8) virus (hereafter B1904) was 3.0, and all of the chickens died within twenty hours, suggestive of high pathogenicity in chickens (Chen et al. 1998). We then tested virulence of chickens challenged with the B1904 virus. Ten chickens were intranasally inoculated with 105 TCID50/200 μl of the B1904 virus. All of infected chickens died within ninety-six hours (Fig. 4B). The organs, including the heart, lungs, kidneys, brain, spleen, and liver of the three dead chickens all tested positive for H5N8 viruses, with an average of 106.5, 105.26, 106, 106.25, 106.17, and 105.92 EID50/200 μl (Fig. 4A), suggestive of the broader tissue tropism in chickens.

Figure 4.

The pathogenicity of avian influenza A(H5N8) viruses in chickens and mice. (A) Viral titers in the heart, liver, spleen, lung, kidney, and brain samples of three infected chickens. (B) Percent survival of infected chickens. (C) Percent survival of infected mice. (D) Viral titers in the heart, liver, spleen, lung, turbinate, and brain samples of three infected mice on 3 dpi. (E) Viral titers in the heart, liver, spleen, lung, turbinate, and brain samples of three infected mice on 5 dpi. (F) Body weight of infected and uninfected mice. (G) HI antibody titers of the serum of H5/H7 chickens (n = 8) vaccinated with clade 2.3.2 and 2.3.4.4 H5 vaccines with H52001/H5N6 virus derived from clade 2.3.4.4h and rHy/H5 virus derived from clade 2.3.2 together with three 2.3.4.4b H5N8 viruses (HB181/H5N8, JS1603/H5N8, and B1904/H5N8) in China during 2014–2020. The name of viruses is the abbreviation of the original name for each virus. H52001/H5N6, A/Chicken/China/H52001/2020(H5N6); HB181/H5N8, A/Chicken/China/HB181/2016(H5N8); JS1603/H5N8, A/Goose/Jiangsu/JS1603/2014(H5N8); B1904/H5N8, A/Duck/Southwestern China/B1904/2020(H5N8). (H) Antigenic map of A/Duck/Southwestern China/B1904/2020(H5N8) viruses together with H5 subtype viruses circulating in China by cartography. The HI data were analyzed by using antigenic cartography (http://www.antigenic-cartography.org), which is a method to visualize and increase the resolution of HI results. Each point on the map represent an HA protein antigen. The distance between two HA protein antigens on the map represents the antigenic distance between the two antigens. Points are colored based on categorical hierarchical clustering. Viral titers in panels A, D, and E were presented as the mean ± SD and analyzed using GrapPad Prism 5.0. The dashed lines indicated the lower limit of detection. EID50, 50% egg infective dose.

Meanwhile, although duck challenged with the B1904 virus showed no clinical symptoms, cloacal and trachea virus shedding was observed in all of the infected ducks at 3 and 5 dpi. Added to that, all of the contacted ducks showed virus shedding on 3, 5, and 7 dpi (Supplementary Table S2), indicating that these H5N8 viruses were highly transmissible in ducks. Although these H5N8 viruses were not lethal to ducks, secondary bacterial or viral infections can lead to the high rates of mortality in domestic and wild aquatic birds (Morris, Cleary, and Clarke 2017), which warrants greater attention.

In addition, 106 TCID50/50 μl B1904 virus was intransally inoculated into six-week-old female BALB/c mice. All of the mice began to lose weight on 2 dpi (Fig. 4F), and all of the mice died on 7 dpi (Fig. 4C), thereby indicating that these H5N8 viruses were lethal to the mice. Virus isolation revealed that the lung and turbinate samples were all positive for H5N8 viruses on 3 and 5 dpi, and heart samples were all positive for H5N8 viruses on 5 dpi. The virus titers in lung samples were far higher than those in the other organs, with averages of 103.08 and 103.5 EID50/200 μl on 3 and 5 dpi, respectively (Fig. 4D–E), which suggests that these H5N8 viruses pose an increasing threat to humans.

3.5 Antigenic properties

We then determined HI titers in the serum of H5/H7 vaccinated chickens (n = 8) and circulating 2.3.4.4h strain H52001 together with three 2.3.4.4b H5N8 strains during 2014–2020. Serum from chickens vaccinated with H52001 showed high titers (7–11 log2) and low titers (0–6 log2) to the 2.3.4.4b H5N8 viruses from 2014 to 2016, and serums from chickens vaccinated with H52001 reacted especially poorly to the B1904 H5N8 virus in 2020 (Fig. 4G; Supplementary Table S3). Beyond that, the B1904 H5N8 virus exhibited antigenicity that is distinct from that of the commercial H5 vaccine and circulating H5 subtype viruses in China (Fig. 4H; Supplementary Table S4), suggesting that antigenic drift of circulating H5 subtype viruses and decreased antibody response against the prevailing H5 subtype viruses. However, the efficiency of protective efficacy of chickens immunized with current commercial vaccines requires further investigation. In China, vaccinating of chickens is the most effective strategy for controlling of avian influenza virus. After an influenza H5/H7 bivalent vaccine for poultry was used in September 2017, the prevalence of the H7N9 viruses in poultry and humans decreased dramatically (Shi et al. 2018; Zeng et al. 2018; Zhang et al. 2020b), and for that reason, new vaccine candidates in China need to be updated to further reduce the dissemination of these H5N8 viruses in poultry. Although vaccinating of chickens can largely control the influenza outbreak; it can also accelerate the generation of novel influenza variants under the vaccine selection pressure. Therefore, the enhancement of biosecurity precautions is also essential to reduce the dissemination of avian influenza viruses.

4. Conclusions

Our results confirmed that H5N8 viruses were co-circulating at low levels in wild birds and domestic poultry in multiple provinces of China from December 2020 to February 2021. Our full-length genome phylogenic analysis revealed that these poultry-origin H5N8 viruses belonged to clade 2.3.4.4b and clustered with H5N8 viruses from some European and Asian countries and were genetically closely related to human-infecting H5N8 virus in Russia. Notably, these H5N8 viruses were not only lethal to chickens and mice but also highly transmissible in ducks. Added to that, they were antigenically distinct from the circulating H5 subtype viruses and commercial vaccine strains in China, which signals alarming threats to the poultry industry and public health. Although the European-origin clade 2.3.4.4b H5N8 viruses were currently observed in a few provinces of China, the risk that they spread to other localities through the transportation of live poultry and wild bird migration and the threat posed to humans clearly merits serious attention in the future.

Click here for additional data file.