The evolution and genetic diversity of avian influenza A(H9N2) viruses in Cambodia, 2015 – 2016

Low pathogenic A(H9N2) subtype avian influenza viruses (AIVs) were originally detected in Cambodian poultry in 2013, and now circulate endemically. We sequenced and characterised 64 A(H9N2) AIVs detected in Cambodian poultry (chickens and ducks) from January 2015 to May 2016. All A(H9) viruses collected in 2015 and 2016 belonged to a new BJ/94-like h9-4.2.5 sub-lineage that emerged in the region during or after 2013, and was distinct to previously detected Cambodian viruses. Overall, there was a reduction of genetic diversity of H9N2 since 2013, however two genotypes were detected in circulation, P and V, with extensive reassortment between the viruses. Phylogenetic analysis showed a close relationship between A(H9N2) AIVs detected in Cambodian and Vietnamese poultry, highlighting cross-border trade/movement of live, domestic poultry between the countries. Wild birds may also play a role in A(H9N2) transmission in the region. Some genes of the Cambodian isolates frequently clustered with zoonotic A(H7N9), A(H9N2) and A(H10N8) viruses, suggesting a common ecology. Molecular analysis showed 100% of viruses contained the hemagglutinin (HA) Q226L substitution, which favours mammalian receptor type binding. All viruses were susceptible to the neuraminidase inhibitor antivirals; however, 41% contained the matrix (M2) S31N substitution associated with resistance to adamantanes. Overall, Cambodian A(H9N2) viruses possessed factors known to increase zoonotic potential, and therefore their evolution should be continually monitored.

Introduction Subtype A(H9) avian influenza viruses (AIVs) circulate globally in wild avian species and are endemic in domestic poultry in many Asian, Middle Eastern and African countries [1,2]. Despite its status as low pathogenic avian influenza (LPAI), subtype A(H9) presents a concern for both the agricultural and health sectors. While infected flocks only experience mild respiratory disease with mortality rates generally below 20% [3,4], A(H9) infections decrease body weight of broilers and the egg production and quality in layers and breeders [5,6]. Additionally, infected poultry are rendered more susceptible to secondary infections [7][8][9], potentially increasing flock mortality levels up to 65% [10,11].
Cambodian AIV surveillance started in response to the detection of highly pathogenic avian influenza (HPAI) A(H5N1) in Cambodian poultry in 2004. These viruses now circulate endemically in domestic poultry species [24,25]. In 2013, surveillance efforts expanded to encompass A(H7) and A(H9) viruses [26][27][28]. It is now evident that A(H9) LPAI viruses circulate endemically in Cambodian poultry [24,28,29] similar to Bangladesh, China and Vietnam [30][31][32] and live bird market (LBM) workers are exposed to these viruses [24]. Therefore, A (H9) AIVs pose a significant risk to Cambodian agricultural and public health sectors. In order to further understand the diversity, origins and molecular risk traits of A(H9) AIVs circulating in Cambodian poultry, we analyzed full genomic sequences of A(H9N2) viruses identified in poultry from LBMs between 2015 and 2016.

Sample collection
Active AIV surveillance in Cambodia was performed by Institut Pasteur du Cambodge (IPC) in collaboration with the National Animal Health and Production Research Institute (NAH-PRI) under the direction of the General Directorate for Animal Health and Production, Cambodian Ministry of Agriculture, Forestry and Fisheries (MAFF). Throughout 2015 and 2016 surveillance for AIVs in poultry was conducted at two prominent LBMs in the Cambodian poultry network: i) Orussey market located in Phnom Penh, the capital city of Cambodia; and, ii) Takeo market, a provincial market located in Takeo province (S1 Fig) [33]. In 2015, samples were collected weekly from February to December [28]. In 2016, collections targeted periods known for having high levels of AIV circulation: Lunar New Year (February), Khmer New Year (April) and Pchum Ben (October). Oropharyngeal and cloacal samples from individual birds were pooled and screened for A(H9) AIVs with qRT-PCR assays sourced from the International Reagent Resource (https://www.internationalreagentresource.org/Home.aspx) as described previously [27,28].

Viral isolation
A subset of A(H9) AIV positive samples, with matrix gene qRT-PCR Ct values <30, were isolated in embryonated chicken eggs (ECEs) [34]. Briefly, original pooled samples were diluted 1:1 with a penicillin-streptomycin PBS solution, filtered (0.22 μM) and then inoculated into the allantoic cavity of 10-12 day old ECEs. Allantoic fluid was collected 48-72 hours post inoculation. The presence of AIV was verified by hemagglutination assay with 0.5% chicken red blood cells followed by qRT-PCR for the AIV MP gene. Negative samples were passaged a minimum of three times in ECEs, after which viral isolation was deemed unsuccessful.

Sequencing of viral isolates
A(H9N2) AIV positive isolates were sent to the WHO Collaborating Centre for Reference and Research on Influenza, Melbourne, Australia for whole genome sequencing (WGS). Prior to sequencing, RNA was extracted from samples with the NucleoMag1 VET kit (Macherey-Nagel) and all genomic segments were amplified in a single reverse transcription PCR (RT-PCR) reaction [35]. The standard manufacturer's protocol was used to obtain whole genome sequences using the Ion Torrent PGM. Briefly, the concentration of RT-PCR products was normalised and fragmented to produce a 200 base-read library using the Ion Xpress™ Plus Fragment Library Kit (Life Technologies). Samples were barcoded with unique Ion Xpress barcode adapters (Life Technologies), pooled and purified using AMPure XP reagent (Agencourt). The final libraries were quantified using the Ion Library Quantitation Kit (Life Technologies) and 20 pM was used to perform an emulsion PCR to enrich Ion Sphere Parti-cles™ (ISPs) on the Ion OneTouch 2 instrument (Life Technologies). The ISPs were then loaded onto a 318™ Chip v2 (Life Technologies) and run on the Ion Torrent PGM.
Quality control of the NGS reads was assessed using CLC Genomic Workbench v10 (QIA-GEN). Low quality reads less than 50 nucleotides in length were removed and a minimum base call quality Phred score of 20 was set. Remaining sequences were aligned to a reference genome (GISAID accession numbers: EPI542434, EPI542395, EPI542405-EPI542410). Genes containing gaps or areas of low sequence coverage were completed using Sanger sequencing. Briefly, segment specific primers were used to amplify appropriate regions using conventional PCR [36]. Sequencing was performed using Big Dye Terminator Reaction Mix (Applied Biosystems) on an ABI 3500xL Genetic Analyzer. The NGS and Sanger sequencing data was assembled to produce consensus sequences using Geneious1 9.1.8 (Biomatters Ltd). Sequence accession numbers for viruses generated during this study are available in S1 Table.

Phylogenetic analysis
Sequences included in the phylogenetic analysis were obtained from GISAID [37], GenBank [38] or the Influenza Research Database [39] and aligned along with Cambodian A(H9N2) viral sequences using MAFFT v7.308 [40]. IQ-Tree was used to produce Maximum Likelihood (ML) phylogenetic trees using the best-fit nucleotide substitution model defined by the Akaike Information Criteria [41]. The General Time Reversible nucleotide substitution model with invariant sites (I) and gamma rate of heterogeneity (GTR + I + Γ) was selected as the best-fit for all datasets, except for MP and NS that used the transversion model (TVM) + I + Γ. Topological support was estimated by 1,000 ultrafast bootstrap replicates [42,43].
Time to the most recent common ancestor (TMRCA) for each node of the HA and NA phylogenies was estimated using a Bayesian Markov Chain Monte Carlo (MCMC) approach with BEAST v1.8.4 run on the CIPRES Science Gateway web portal [44]. For each gene the GTR + Γ and SRD06 nucleotide substitution model was used [45] with a relaxed uncorrelated log-normal molecular clock [46]. Population dynamics was investigated by using a coalescent Gaussian Markov random field Bayesian skyride tree prior with time-aware smoothing [47]. Maximum clade credibility (MCC) tree for each gene was the result of two independent analysis runs for 100 million generations, sampled to produce 10,000 states with 10% removed as burn-in and combined using TreeAnnotator v1.8.4. The ML and MCC trees were formatted using FigTree v1.4.3 [48].
Reassortment was visualised by producing a phylogenetic congruency map, whereby segments of the individual viruses were linked across the ML phylogenies. The ML trees were created using IQ-tree as described above with a subset of viruses that had sequencing data available for all eight genomic segments (S2 Table). Incongruence is demonstrated by deviations in the topology and can indicate that reassortment has occurred.

Viral genotyping
For the purpose of this study, the genotyping system used follows conventions initially described by Li et al., 2005 [3].  [49]).

Molecular analysis
The H3 HA and N2 NA numbering systems were used throughout the text, unless otherwise stated. Molecular markers associated with changes in viral fitness or resistance to antivirals were investigated using the molecular inventory matrix produced by Suttie et al., 2019 [50]. Nlinked glycosylation sites for HA and NA were predicted using the NetNGlyc1.0 server with default settings [51].

Selection pressure analysis
Selection pressures acting on each gene of the Cambodian A(H9N2) AIVs was analysed by calculating the ratio of non-synonymous to synonymous mutations (dN/dS ration, ω) using the HyPhy software package [52] accessed via the datamonkey webserver [53]. Selection pressure is interpreted based on the value of ω: ω < 1 indicates negative selection, ω = 1 neutrality and ω > 1 positive selection. Programs used to infer selection pressure included: fixed-effects likelihood (FEL), fast unconstrained Bayesian approximation (FUBAR), mixed effects model of evolution (MEME) and single-likelihood ancestor (SLAC) [54][55][56]. Statistically significant sites detected using two or more programs with a p-value < 0.1 for FEL, MEME or SLAC, or a posterior probability �0.90 for FUBAR were considered valid.

Neuraminidase inhibition assay
The susceptibility of 40 Cambodian A(H9N2) viruses to four NA inhibitors (NAIs; oseltamivir, zanamivir, laninamivir and peramivir) was evaluated using a standard fluorescence based assay that measures NA enzymatic activity to determine the 50% inhibitory concentration (IC 50 ) of viruses. This protocol has been described in detail by Leang et al., 2017 [57]. The IC 50 values were calculated using the JASPR v1.2 software (CDC, Atlanta, USA). NAI susceptibility is described based on the fold change in IC 50 values of the test viruses compared to control viruses; viruses with normal inhibition had <10-fold increase in their IC 50 , reduced inhibition had between 10 to 100-fold increase in IC 50 or highly reduced inhibition >100-fold increase in IC 50 .

Ethical approval
Animal sampling was conducted by the National Animal Health and Production Research Institute under the direction of the General Directorate for Animal Health and Production, Cambodian Ministry of Agriculture, Forestry and Fisheries as part of routine disease surveillance activities; thus, poultry sampling was not considered as experimental animal research. The analysis of poultry samples for avian influenza testing was conducted by the Virology Unit at Institut Pasteur du Cambodge as part of routine activities as a World Health Organization H5 Reference Laboratory.

Phylogenetic diversity of the Cambodian A(H9N2) viruses
A total of 64 A(H9N2) viruses (49 in 2015; and 15 in 2016) were isolated in ECEs, with 43 from chickens and 21 from ducks (S1 Table). Bayesian estimates of the TMRCA for the Cambodian A(H9N2) H9 HA and N2 genes indicate that the H9 HA genes share a more recent common ancestor than the N2 genes, demonstrating there may have been a greater turnover of H9 HA lineages compared to N2. Additionally, Bayesian skyride analysis of genetic diversity shows greater fluctuations in H9 HA diversity compared to N2 (S3A- S3B Fig). Overall, the Bayesian MCC and ML phylogenies reveal there has been reassortment between persistent NA lineages with newly introduced HA segments into Cambodia.
Similar to HA and NA, internal genes of the Cambodian A(H9N2) viruses also split into either two or three main groups (S4A- S4F Fig). However, group composition is not consistent across the different segments, further indicating extensive reassortment between co-circulating lineages. This is also evident by the phylogenetic incongruence shown between segments of many viruses (Fig 3; nucleotide sequence identities: PB2: 83.6-100%, PB1: 92.9-100%, PA: 92.7-100%, NP: 90.8-100%, MP: 93.5-100% and NS: 92.7-100%). For instance, the virus A/ chicken/Cambodia/a58W22M1/2016 displays incongruence (Fig 3). This is a genotype V virus with congruence links in Fig 3 shown in green, however its HA and PB1 genes clearly deviate to cluster with genotype P viruses shown in purple.
All Cambodian A(H9N2) gene phylogenies, excluding NA, cluster with lineages containing zoonotic A(H7N9), A(H9N2) and A(H10N8) AIVs (Figs 1, 2 and S4A-S4F). Typically, Cambodian A(H9N2) AIVs did not cluster closely with the genotype S zoonotic viruses. However, a single isolate, A/chicken/Cambodia/a77W22M1/2016, had PB2, PA, MP and NS genomic segments similar to AIVs that have infected humans (nucleotide sequence identities >98%; S4A, S4C, S4E and S4F Fig). This virus clustered separately from other Cambodian A(H9N2) isolates, representing a possible separate introduction of A(H9N2) into Cambodia. This cluster contained AIVs from China and Japan, as well as waterfowl hosts in Vietnam, and is most closely related to A(H9N2) AIVs detected in Indonesia, with all genomic segments having a nucleotide sequence identity >99%.

Cambodian A(H9N2) viruses have molecular markers associated with viral adaptation to mammals
Screens for mutations known to affect AIV pathogenicity (S3A-S3H Table), are summarised in Table 1 [50]. All Cambodian viruses had HA cleavage site motifs common for A(H9) LPAI  viruses: PSKSSR/GLF (n = 55) and PSRSSR/GLF (n = 9) [60]. The HA genes also contained a number of amino acid substitutions associated with an increase in AIV binding to humantype α2,6 receptors (Tables 1 and 32d), including Q226L in 100% of Cambodian A(H9N2) isolates.
Cambodian A(H9N2) viruses contained numerous substitutions associated with increased polymerase activity of AIVs in mammalian cell lines (Tables 1, S3A-S3C and S3E). However, well known substitutions in PB2, namely E627K and D701N, were not detected. The majority (n = 61) of viruses had full length, 90 aa PB1-F2. A single virus, A/chicken/Cambodia/ a77W22M1/2016, had a minor PB1-F2 truncation and was 76 aa long. Additionally, two viruses had major truncations with premature stop codons at position 9. All Cambodian viruses had full length PA-X proteins at 252 aa in length.
Only a small number of mutations in MP and NS1 were detected that increase AIV virulence in mammalian or avian models (Tables 1 and S3G-S3H). Cambodian A(H9N2) viruses contained either truncated (217 aa, n = 36) or N-terminal elongated (237 aa, n = 28) NS1 protein. NS1 elongation likely arose in a specific cluster within the Cambodian A(H9N2) population (S4F Fig). No molecular markers of concern were identified in NS2 (Tables 1 and S3H). All MP genes contained I43M and T215A, associated with increase virulence in mice (Tables 1  and S3G).

Similarity of Cambodian A(H9N2) viruses to available vaccine strains
Genetic analysis of known A(H9) HA antigenic sites [87] showed high similarity between Cambodian A(H9N2) viruses. Of the 33 sites investigated, Cambodian isolates had a maximum of 4 amino acid differences (S4 Table) Table).  Table). Substitutions at the majority of these sites are not associated with increased AIV fitness. However, substitution of serine for asparagine (S158N) at HA codon 158 has been shown to increase AIV binding to human-type α2,6 receptors in A(H5) viruses [75,87]. Additionally, PB2 K389R increases polymerase activity and replication of A(H7) AIVs in mammalian cell lines [62]. These phenotypes have not been examined in A(H9) AIVs and require further investigation.

The susceptibility of Cambodian A(H9N2) viruses to antiviral drugs
Molecular analysis of the M2 protein from Cambodian A(H9N2) viruses showed 41% (n = 26) had the S31N marker associated with resistance to adamantanes [84]. Mapping this mutation to the MP phylogeny indicates that the S31N has arisen in the population on multiple independent occasions (S4E Fig). None of the Cambodian viruses contained PA I38T/F/M mutations associated with resistance to the recently licenced cap dependent endonuclease antiviral, baloxavir marboxil [88,89]. None of the Cambodian A(H9N2) AIVs harboured known N2 mutations associated with resistance to NAIs. Susceptibility to NAIs was confirmed in 40 Cambodian A(H9N2) viruses with IC 50 values <10 against oseltamivir, zanamivir, laninamivir and peramivir (S7 Table).

Discussion
Since integration of detection protocols in 2013, it is clear that subtype A(H9) AIVs now circulate endemically in Cambodian domestic poultry and they are currently the most prevalent subtype detected [29]. Between 2013 and 2016, the genomic diversity of Cambodian A(H9) HA genes decreased, with 2015-2016 viruses belonging to a single clade, BJ/94 H9-4.2.5 ( Fig  1). In comparison, NA genomic diversity remained relatively stable. However, this could be an artefact caused by the limited number of Cambodian A(H9N2) NA sequences available prior to 2015. The HA, NA and internal gene phylogenies show that these viruses possibly arose from AIVs detected elsewhere in Asia, and are closely related to AIVs concurrently detected in Vietnam. Similar observations between endemically circulating Cambodian and Vietnamese A(H5N1) AIVs also exist. Transmission of A(H5) and A(H9) AIVs between the two countries is likely to be facilitated by cross-border movement of domestic poultry. Transmission by migration of wild birds may also play a minor role in viral diversity in Cambodia. A single A(H9N2) isolate, designated A/chicken/Cambodia/a77W22M1/2016, clustered separately from other Cambodian viruses, possibly representing a separate A(H9) introduction via migration of wild birds along the Eastern Asia/Australasian flyway. This virus clustered with A(H9N2) AIVs detected in poultry from Vietnam between 2014 and 2018, Japan in 2016 and with Indonesian A(H9N2) AIVs detected from 2016 to 2018. The first detection of A(H9N2) in Indonesia occurred during AIV outbreaks in 2016 and was associated with increased flock mortality and an 18% reduction in the egg production of layers [90]. Further surveillance for A(H9) AIVs in domestic poultry and wild birds is vital to understand the dynamics of A(H9) transmission between these countries.
Since 2013, A(H9) genotype S viruses have been the most prevalent AIV lineage detected in Chinese chickens [18]. The enhanced fitness of genotype S in chickens is partly attributable to the acquisition of G1/97-like MP and PB2 genes. These genes increase AIV polymerase activity, replicative capacity, and virulence of AIVs in chickens [91,92]. Both Cambodian genotype P and V viruses contained G1/97-like matrix proteins, but only genotype V also has a G1/ 97-like PB2. Interestingly, the prevalence of genotype V in Cambodian poultry increased from 2015 to 2016; however, this finding may be an artefact of limited sampling performed in 2016 compared to 2015.
The internal genomic cassette of genotype S viruses have been repeatedly donated to zoonotic AIVs such as A(H5N6), A(H7N9) and A(H10N8) [18][19][20][21]. In Cambodia, genotype S viruses have not been detected, although A(H9N2) viruses have multiple segments that cluster with genotype S. This finding is unsurprising considering that the genotype V viruses, originally described in China 2014, arose from the reassortment of genotype S viruses that acquired G9/97 NA genes [1,49]. Interestingly, the Indonesian-like isolate from 2016 (A/chicken/Cambodia/a77W22M1/2016) contained PB2, PA, MP and NS genes that clustered closely with zoonotic A(H7N9), A(H9N2) and/or A(H10N8) AIVs (nucleotide sequence identity >98%). Consequently, introduction of new clades and continual reassortment combined with the endemicity of A(H9N2) viruses in Cambodian poultry raises concerns about the emergence of novel zoonotic AIVs.
Molecular analyses indicate the Cambodian A(H9N2) AIVs were LPAI viruses with multiple molecular markers associated with adaptation to mammalian species. Cambodian A (H9N2) viruses contained multiple HA markers associated with increased binding to "humanlike" α2,6-linked sialic acid residues (Table 1). One well known HA substitution, Q226L, significantly increases the binding of A(H9) AIVs to α2,6-linked sialic acids and enhances AIV transmission in ferrets [30,79]. Between 2013 to 2016, the prevalence of Q226L in Cambodian A(H9N2) viruses increased from 50% to 100% similar to BJ/94 and G1/97 lineage A(H9) viruses in China, Vietnam and the Middle East [1,30,93]. The rising global prevalence of Q226L is concerning considering the enhancement of zoonotic potential. LBM studies performed in Cambodia in 2013 showed 1.8% of LBM workers had seroconverted to the A(H9) subtype [24]. Further work is needed to see if increased subclinical, zoonotic transmission has occurred with this change in receptor specificity.
Other substitutions identified in internal proteins of Cambodian A(H9N2) viruses have also been associated with increased AIV replication and virulence in mammals (Table 1). While major molecular markers of AIV adaptation to mammals in PB2, namely E627K and D701N, were not identified, a small number of isolates contained A588V (n = 1) and Q591K (n = 3) substitutions that partially compensate for the absence of E627K [61,63]. In particular, A588V increases the polymerase activity and replication of A(H7N9) and A(H9N2) viruses in mammalian and avian cell lines, consequently increasing AIV virulence in mammalian models [61]. Since 2013, the number of A(H9N2) and zoonotic A(H7N9) viruses with A588V substitutions has increased significantly at the global level [61].
Candidate A(H9) vaccine viruses are selected as part of pandemic preparedness plans [94]. The Cambodian A(H9N2) viruses have the highest nucleotide identity to A/Hong Kong/308/ 2014 (94.1% to 95.1%), though they are more similar to G9/97 at HA antigenic sites. The most appropriate CVV to protect against Cambodian A(H9) viruses would need to be determined experimentally. Additionally, vaccination of poultry, in combination with biosecurity programs, can effectively decrease AIV circulation in poultry flocks and therefore decrease human exposure. However, the most commonly used commercially available A(H9) poultry vaccines in Asia, are produced from viruses detected in the 1990s [95]. As a result, A(H9) poultry vaccine effectiveness continually decreases due to ongoing antigenic drift [95]. No official poultry vaccination against AIVs occurs in Cambodia. Cambodian viruses have the highest sequence identity to SS/94 (91.8% to 90.5% nucleotide identity), though the HA genes differ at six to nine HA antigenic sites. Therefore, it is possible that this variation could reduce vaccine effectiveness, but further experimental work is warranted.
Assessing the effectiveness of available antivirals is vital for combating zoonotic AIVs as part of influenza pandemic preparedness strategies [94]. Molecular data from Cambodian A (H9N2) viruses indicates that 41% of isolates contained the S31N mutation in the M2 protein, indicative of resistance to adamantanes [84,96]. Furthermore, reassortment between these A (H9N2) viruses and endemically circulating A(H5N1) in Cambodia has resulted in a A(H5N1) clade 2.3.2.1c virus with the MP gene of A(H9) origin containing the M2 S31N marker [97]. As this A(H5) clade does not typically contain markers of adamantine resistance, this is one example of a reassortment event increasing AIV risk [1,30]. However, adamantanes are not recommended for treating influenza as resistance is widespread in both seasonal and other potentially zoonotic AIVs [98]. All Cambodian viruses tested remain susceptible to the first line NAI antivirals: oseltamivir, zanamivir, laninamivir and peramivir.
In summary, Cambodian A(H9N2) AIVs isolated from domestic poultry in LBMs between 2015 and 2016 belonged to the BJ/94 H9-4.2.5 lineage. They were closely related to viruses identified in surrounding countries, suggesting frequent circulation of these viruses possibly through cross-border movement of live poultry. Wild birds may also play a role in A(H9) transmission. The Cambodian A(H9N2) viruses contained multiple genomic segments clustering with genotype S lineage AIVs, with a small number close to zoonotic A(H7N9), A(H9N2) and A(H10N8) viruses. Additionally, all Cambodian A(H9N2) viruses had molecular markers associated with viral adaptation to mammalian species, raising further concerns about their increased zoonotic potential. However, these viruses remained susceptible to therapeutic and prophylactic prevention strategies utilizing NAIs. Overall, the continued endemicity and evolution of A(H9N2) AIVs in Cambodian domestic poultry underlines the need for continued, vigilant surveillance analysis of AIVs within the country.