A total of 1,270 oropharyngeal and cloacal swabs were collected from chickens, ducks, geese and pigeons in a live-poultry market between 2018 and 2021. Each sample was placed in 1 mL of cold phosphate-buffered saline (PBS) containing penicillin (5,000 U/mL) and streptomycin (5,000 U/mL). After mixing and centrifugation at 10,000 × g/min for 5 min, 0.2 mL of supernatant was used to inoculate 9-day-old specific-pathogen-free chicken embryos via the allantoic cavity, followed by incubation at 37°C for 48–72 h. We then harvested the allantoic fluid, and a total of 29 viruses were isolated (Supplementary Table 3).

RNA was extracted from the harvested virus suspensions using the RNeasy Mini Kit (Qiagen, Germany), and the entire full genome sequence was amplified using two-step RT-PCR and the universal primers reported by Hoffmann et al. (2001). The ex Taq premixed enzyme used for the amplification was purchased from the TaKaRa reagent company.¹ PCR products of the HA and NA segments of these viruses were subjected to agarose gel electrophoresis, and the target fragments were recovered using the QiAamp Gel Extraction Kit (Qiagen, Germany)² and sequenced using an ABI3730 DNA Analyzer (Shenggong Bioengineering Co., Ltd., Shanghai, China). Data were merged and assembled using the Lasergene software³ based on the National Center for Biotechnology Information (NCBI) virus database⁴ (accessed on 12 June 2022).

The MAFFT version 7.058 was used to align each of the HA and NA gene segments and eliminate those sequences with less than 95% of the expected segment length. Duplicate sequences in the gene fragment were removed using PhyloSuite. We performed the phylogenetic analysis three times using the ML method in IQ-TREE under the GTR + F + G4 model with 1,000 bootstrap replications. A high-quality visualization of the phylogenetic data was performed using the Interactive Tree of Life (iTOL).

Based on the phylogenetic topologies obtained and their bootstrap values, we selected several representative reference sequences and formed five smaller datasets. TempEst (version 1.5.1) was used to analyze the R² values and correlation coefficient of the temporal signals and the best-fit model in the selected sequences. To estimate the nucleotide substitution rates of HA and NA segments, we used the Bayesian Markov chain Monte Carlo (MCMC) method provided in the Bayesian Evolutionary Analysis Sampling Trees (BEAST) (v1.10.4c) and a relaxed molecular clock model with uncorrelated log-normally distributed rates and a coalescent Bayesian skyline plot. We set the chain lengths to 500 million iterations and performed sampling at every 5,000 steps to obtain an effective sample size (ESS) ≥200, and convergence was assessed using Tracer (v1.7.1). Time-scaled summary of MCC trees with 10% for the post-burn-in posterior were created using TreeAnnotator (v1.10.4), and visualized with FigTree (v1.4.4). The BEAST package was used for the calculation of the tMRCAs of each branch.

The H9N2 viruses that we isolated in the live-poultry market were distributed in three independent branches; therefore, we performed a spatiotemporal analysis of clades A, B, and C, respectively. In order to reduce the potential sampling biases, we randomly subsampled the database in a stratified manner to create a more equitable spatio-temporal distribution of the HA genome sequences of three branches of viruses. To be precise, sequences in each branch were clustered using the CD-HIT program (Huang et al., 2010), and identical sequences within the same time and region were removed. The discrete sampling locations of the clade A, B, and C viruses in this study include Guangdong, Yunnan, Jiangxi, Shandong, Shanghai, Fujian, Jiangsu, Hunan, Henan, Hebei, Hubei, Xinjiang, Ningxia, Chongqing, Guizhou, Guangxi, Sichuan, Shanxi, Beijing, Tianjin, Heilongjiang, and Anhui in China, there are also viruses from Vietnam in clade C. Detailed information regarding the subsampled HA gene sequences of the clade A, B, and C H9N2 subtype viruses used in this study is provide in the Supplementary Table 1.

Time-measured phylogenies were inferred using the Bayesian discrete phylogeographic approach implemented in the BEAST package (v1.10.4). We first performed a regression of root-to-tip genetic distances on the ML tree against exact sampling dates using the TempEst v1.5.3 (Rambaut et al., 2016), which showed a strong temporal signal. Then, we used an uncorrelated lognormal (UCLN) relaxed molecular clock model. In addition, a Bayesian stochastic search variable selection (BSSVS) model with asymmetric substitution was used. For each independent dataset, multiple runs of the MCMC method were combined using LogCombiner (v1.10.4), utilizing 5,000,000,000 total steps for each set, with sampling every 500,000 steps. Subsequently, we used SpreaD3 v0.9.7.1 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 (Bielejec et al., 2016). The BF values >100 indicated robust statistical support, 30 < BF values ≤100 indicated very strong statistical support, 10 < BF values ≤30 indicated strong statistical support, 3 < BF values ≤10 indicated substantial statistical support and BF values <3 indicated poor statistical support (Lemey et al., 2009). We used QGIS Version 3.28 to create plots showing the results of the BF tests.⁵

The selective pressure for each gene segments in clades A, B, and C was determined on the Datamonkey online version of HyPhy package.⁶ The more appropriate number of sequences (n ≥ 200) we selected were depend on a study of Ji-Ming et al. (2022). All the sequences in every branch were analyzed by ratio estimation of non-synonymous (dN) to synonymous (dS) substitutions (ω = dN/dS) on a codon-by-codon basis, and ω < 1 indicates negative or purifying selective pressure; ω = 1 implies neutral evolution; and ω > 1 shows positive selection. BioEdit and MEGA 6.0 (Tamura et al., 2013) were used for the analysis of the sequence format conversion and key amino acid changes, and MegAlign was used to analyze the sequence homologies.⁷ We used Weblogo to analyze the conservation of amino acids.⁸ The prediction of potential N-glycosylation sites was performed using the NetNGlyc server 1.0.⁹

To explore the prevalence of H9N2 in various regions of China, we counted the percentage of H9N2 avian influenza viruses isolated every year in each region. We found that the number of field isolates of H9N2 in East China accounted for 62.1% of the total isolates in China in 2012–2013. Then this percentage began to decline, with isolates dropping to about 42% in 2018–2019. The isolate percentage of H9N2 AIVs in South China increased to 19.6% in 2014–2015, and then decreased and remained at 13%. The percentage in Southwest China increased from 2% in 2012–2013 to 29.3% in 2018–2019 (Figure 1). The decrease in the percentage of isolates in Guangdong is most likely related to live-poultry markets implementing the “1110 system” of “cleaning and disinfection once a day, a comprehensive sweep once a week, off the market once a month and live-poultry retail markets’ zero stock of live birds on the day of the market.”

We isolated 29 strains of H9N2 avian influenza virus from the same live-poultry market, and the nucleotide identity of HA sequences of all isolates ranged from 92.28 to 100% and the amino acid identity ranged from 93.58 to 100%. The isolated strains were clustered into three groups, I–III (Supplementary Figure 1). Genetic evolutionary analysis revealed three clusters located in branches A, B, and C (the average intergroup distance between branches A and B and C was 0.057 and 0.065, respectively, while the average intergroup distance between B and C was 0.055, which can be considered as forming three independent branches) (Figure 2). Our study found that clade A was first isolated from duck flocks in Zhejiang, and this branch was isolated in live-poultry markets in 2018; however, no further strains of this branch were found in a subsequent surveillance. The clade B branch became endemic in central China in 2016 and then spread to Southern and Southwestern China, and multiple strains of clade B were isolated in live-poultry markets in 2019–2020. Clade C strains were detected in increased numbers in Southern China in 2016 and became epidemic in Guangdong starting in 2018, we continued to isolate strains located in clade C in live poultry markets for 4 years (Figure 2A). We also found that Human-derived H9N2 AIVs does not form a separate branch in the phylogenetic evolutionary tree, and multiple H9N2-induced human infection strains occurring across the country in 2019–2021 are located in clades B and C. According to the chronological order of the prevalence of the strains located in these branches in various regions, the H9N2 AIVs prevalent in Guangdong in recent years may be spread mainly by East China, while the H9N2 viruses prevalent in the Foshan live-poultry market are likely to be spread by strains prevalent in other regions that conduct live-poultry trade.

Genetic evolutionary analysis revealed that the evolutionary tree of the NA gene of the H9N2 viruses, which used to be endemic in China, diverged in 2012–2013 and evolved into three major branches, clade B, clade C and clade D (from A/Chicken/SD/6/96). Among them, clade D strains have rarely been isolated since 2017, while clades B and C are widespread in China, with a wide range of host adaptations including poultry, mammals, and even multiple cases of human-derived infections (Figure 2B). We also isolated many strains located in clade A (belonging to branch A/duck/HK/Y439/97), which has been detected in low numbers in China and is present mainly in wild birds and waterfowl. The nucleotide identity of the NA genes of our isolated strains ranged from 81.49 to 100% and amino acid identity from 85.11 to 100% in the clade A and clade B branches, respectively. The 12 isolates located in clade A were all derived from ducks and geese and did not have neck deletions of the NA protein, consistent with the Y439 branch strain being prevalent mainly in waterfowl. The 16 isolates located in clade B were isolated from chickens, with only one from geese, consistent with the prevalence of Shandong/6/96 in poultry and with the finding that the absence of the NA protein neck enhances the adaptation of the virus from wild birds to poultry.

To analyze the evolutionary relationships of clades A, B, and C with the Chinese H9N2 as a whole, we reconstructed the evolutionary dynamics and estimated the divergence times of the three branches. Root-tip regression analysis of the structure revealed a clock-like structure of the HA gene of the H9N2 subtype of avian influenza virus prevalent in China from 2016 to 2021 (n = 189, correlation coefficient = 0.86; R2 = 0.73) (Figure 3A). Our analysis found that the population diversity of the H9N2 subtype prevalent in China was increasing from 2016 to 2017 and peaked in 2017 (Figure 3B), mainly due to the clade 1 branch beginning to diverge into multiple sub-branches during this period. Genetic diversity decreased significantly starting in 2018, likely due to the dominant branch being dominant, corresponding to the widespread clade C branch prevalence. Previous studies have shown that the evolution rate of the H9N2 virus has increased in recent years (Jin et al., 2020; Wang et al., 2021), which implies an increasing rate of mutation. We estimated the overall average evolution rate of HA genes of H9N2 viruses in China from 2016 to 2021 to be 4.538 × 10^–3 substitutions/site/year (95% HPD, 4.00 × 10^–3 to 5.12 × 10^–3 substitutions/site/year) and the estimated time of tMRCA to be September 2008 (95% HPD, July 2006–2010). Meanwhile, we calculated the evolutionary rates of the three different branches where the HA genes of the isolates were located separately and found that the evolutionary rate of clade A is 5.56 × 10^–3 substitutions/site/year (95% HPD, 4.50 × 10^–3 to 6.64 × 10^–3 substitutions/site/year) and the estimated time of tMRCA to be July 2012 (95% HPD, March 2012–November 2012). The evolutionary rate of clade B is 5.32 × 10^–3 substitutions/site/year (95% HPD, 4.44 × 10^–3 to 6.24 × 10^–3 substitutions/site/year), and the estimated time of tMRCA is March 2012 (95% HPD, April 2009–February 2014). The evolutionary rate of clade C is 5.14 × 10^–3 substitutions/site/year (95% HPD, 4.16 × 10^–3 to 6.10 × 10^–3 substitutions/site/year), and the estimated time of tMRCA is August 2015 (95% HPD, November 2014–April 2016) (Table 1 and Figure 3C). Our study shows that 2012–2013 and 2014–2016 are two critical periods for the divergence of the H9N2 subtype AIV in China, and during this period, the evolutionary rate of HA genes was maintained at a high level. The difference between the evolutionary rate of the HA gene as a whole and the evolutionary rate calculated for each branch individually suggests that the branches of H9N2 are not evolving at the same rate; some branches are spreading widely and have a higher probability of communicating with other strains, while the use of vaccines in widespread areas may also lead to a higher probability of the virus to mutate to escape immunity.

We reconstructed the spatial dispersal networks of the three branches based on the HA genes and found that the clade A viruses spread mainly among Chinese provinces, with Shanghai → Hunan (BF = 12), Shanghai → Jiangxi (BF = 33), Shanghai → Zhejiang (BF = 47), Hunan → Shanxi (BF = 207) and Jiangxi → Guangdong (BF = 22) as the main routes (Table 2 and Figure 4A). We also found that Jiangxi → Guangdong (migration rate = 1.54) had a higher migration rate, and that Shanghai, Jiangxi, Jiangsu, and Hunan are the main prevalent geographic locations of this branch and closely connected with other regions (Supplementary Table 2). The clade B viruses spread mainly in the trajectories of Jiangxi → Jiangsu (BF = 3,386), Jiangxi → Guangdong (BF = 849), Jiangxi → Ningxia (BF = 48), Jiangxi → Shanxi (BF = 52), and Hunan → Henan (BF = 112) and Hunan → Guizhou (BF = 62) (Table 2 and Figure 4B). The strain we detected in Foshan was probably transmitted from within Guangdong Province (BF = 849), and it can be shown that Jiangxi and Hunan are the main endemic areas of the clade B viruses, with strong links to other regions. The clade C viruses spread mainly with the trajectories Guangdong → Fujian (BF = 161), Guangdong → Guizhou (BF = 23), Guangdong → Hunan (BF = 39), Guangdong → Yunnan (BF = 58), Fujian → Jiangxi (BF = 15) and Yunnan → Xinjiang (BF = 65) (Table 2 and Figure 4C), and the results show that Guangdong is the epidemiological center of this branch and is closely linked to several regions. The strains isolated from our sampling site in Foshan are also very likely to be locally transmitted in Guangdong. Our findings suggest that different branches of the strain are prevalent in South and East China, respectively, with the clade A and clade B branches predominant in Shanghai, Jiangxi, and Hunan, and the clade C branch predominant in Guangdong, which is consistent with the need for different H9 subtype vaccines in each region of China and may also be the primary reason H9N2 AIV has been difficult to control in China.

Positive selection leads to an increase in the number of genetic variants, allowing the preservation of mutations that favor viral adaptation and survival, while negative selection purifies unfavorable mutations and tends to conserve genes (Yang and Nielsen, 2002; Kosiol et al., 2006). We used the Branch model aBSREL (Smith et al., 2015) to analyze the HA genes clades A, B, and C for selection pressure and found that they were all under purifying selection pressure (Clade A, 0.102; B, 0.117; C, 0.161, p = 1). Then, when we analyzed with the site model MEME (Murrell et al., 2012), FUBAR (Murrell et al., 2013), SLAC (Kosakovsky and Frost, 2005), we identified a total of 12 positively selected pressure sites for the HA gene of H9N2 AIVs prevalent in China from 2016 to 2021, and clades A, B, and C have 9, 6, and 8 positive selection pressure points, respectively (Table 3 and Figure 5). Mutations in amino acids at these sites facilitate the adaptation of H9N2 avian influenza viruses to unfavorable environments.

Receptor-binding preference is important for influenza virus replication and transmission (Herfst et al., 2012), and mutations in multiple amino acid sites located in the HA protein affect the receptor-binding propensity of H9N2 AIVs. We found that 13 of the receptor-binding sites for the HA proteins of H9N2 viruses isolated in China from 2016 to 2021 were non-conserved (Figure 6). Two of the sites are not conserved in clade A; three sites are not conserved in clade B, and up to eight sites are not conserved in clade C. Compared to the other two branches, clade A had the fewest number of mutations at key sites. At the same time, we found that strains located in clade B and clade C have mutations at 145 from D to G, 190 from A to T/V, 160 from A to E/D/N, and 156 from Q to R. These sites are conserved in clade A strains, and mutations increase the binding properties of the strains to human receptors. These findings are consistent with clade B and clade C being responsible for human infection and having a greater prevalence.

We found amino acid polymorphisms in several sites of HA proteins located in different branches (Supplementary Table 4), including key sites 145, 155, 160, and 190 that affect receptor binding, and the presence of polymorphisms at these sites may contribute to the differential ability of the corresponding strains to infect humans. Combined with the selective pressure analysis, we found that the 190 point of clade A, the 160 point of clade B, and the 156 and 190 key receptor-binding sites of clade C are under positive selection. It is possible that the wide- spread use of vaccines for prevention of H9N2 AIVs has resulted in several antigenic sites under positive selection pressure, at the same time, some sites are also receptor binding key sites and their mutations facilitate viral adaptation to new hosts (Supplementary Table 5 and Figure 6).

The cleavage sites of clade B and clade C include the forms PTRSSR↓GLF, PSRYSR↓GLF, RSRYSR↓GLFPSKSSR↓GLF, PSRYSR↓GLF and PSRSGR↓GLF, in addition to the PSRSSR↓GLF possessed by clade A (Supplementary Table 5). Mutations in receptor-binding sites and the diversity of cleavage sites may lead to variability in the transmission efficiency of different branches of the strains across multiple hosts, which in turn affects the prevalence of this branch. we found that the HA proteins of clades A, B, and C were not only missing 215–218 potential glycosylation sites but also a new 313–315 glycosylation site near the cleavage site compared to previous reports (Li et al., 2005; Dong et al., 2011b; Liu et al., 2016). Increased glycosylation sites may significantly increase the ability of viruses to infect mammalian cells and avian species (Peng et al., 2019; Yang et al., 2021), and also increase the ability of viruses to escape immunity to vaccines (Tate et al., 2014).

It has been shown that deletion of the NA protein neck enhances the adaptation of the viruses from wild birds to poultry and also enhances the replication of the viruses in mammals. Phylogenetic analysis of the NA gene showed multiple cases of human infection in both clade B and clade C with neck deletion, but no human infection has been reported in clade A. Meanwhile, our analysis of key sites revealed no mutations at points 274 and 294 in clade A and clade B, maintaining susceptibility to oseltamivir and neuraminidase inhibitors.

We analyzed the HA genes of viruses isolated from live-poultry markets and found that all isolates had the PSRSSR↓GLF cleavage site, consistent with a molecular signature of low pathogenicity. The Q226L and H183N mutations were present in all isolates; the I155T mutation was present in 28 isolates; one isolate had the I155N mutation; and 27 isolates had the A190T mutation. A new glycosylation site, 313–315, was present near the cleavage site, indicating that all these isolates are potential risks for human infection. We found that all 12 isolates located in clade A had ITE at sites 59–61 of the NA protein, and 17 isolates located in clade B had missing sites 59–61, indicating that the live-poultry market has a mixture of NA protein neck-deficient and non-deficient strains of H9N2 AIVs, thereby giving the viruses a greater chance of co-transmission and circulation in wild waterfowl and poultry.