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) n class="Species">avian influenza viruses (n>n class="Gene">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-9], potentially increasing flock mortality levels up to 65% [10,11].

A(H9) n class="Gene">AIVs in poultry pose a n>n class="Disease">zoonotic risk to humans [12]. Since the first official report of a human A(H9N2) case in 1998, 59 additional cases have been reported from Bangladesh, China, Egypt, Pakistan and Oman [13]. Human A(H9) infections are generally asymptomatic or manifest as mild respiratory disease. Only one A(H9N2)-associated death has been reported, likely due to comorbidities [14]. All A(H9N2) viruses detected in humans are genetically similar to poultry viruses from the same period [13,15-17]. A(H9) AIVs also commonly donate internal genomic segments to non-A(H9) viruses through reassortment, creating unique AIV genotypes that may present a threat for agricultural sectors and human health [18-21]. This phenotype has been observed for A(H5N1), A(H5N2), A(H5N6), A(H7N7), A(H7N9) and A(H10N8) subtype AIVs [18,22,23]. Of particular interest, A(H5N6), A(H7N9) and A(H10N8) AIVs all acquired internal segments from genotype 57 (G57) A(H9) viruses before emerging as zoonotic viruses [18-22]. G57 is generally equivalent to genotype S [4], the most prevalent A(H9) genotype detected in China since 2013 [1,18].

Cambodian n class="Gene">AIV surveillance started in response to the detection of highly pathogenic avian n>n class="Species">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-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-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.

Materials and methods

Sample collection

Active n class="Gene">AIV surveillance in Cambodia was performed by Institut Pasteur du Cambodge (IPC) in collaboration with the National Animal Health and Production Research Institute (NAHPRI) 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 n>n class="Gene">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) n class="Gene">AIV positive samples, with matrix gene qRT-PCR Ct values <30, were isolated in embryonated n>n class="Species">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) pan class="Gene">AIV positive isolates were sent to the WHO Collaborating Centre for Reference and Research on pan class="Species">Influenza, Melbourne, Australia for whole genome sequencing (WGS). Prior to sequencing, RNA was extracted from samples with the NucleoMag® 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 Particles™ (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 (QIAGEN). 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 Geneious® 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 pan class="Species">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]. Using this system A(H9) n class="Gene">AIVs are sepan>rated into genotypes designated A through to W based on the phylogenetic relationship of genes to eight reference viruses: A/turkey/Wisconsin/1/1966 (WI/1/66), A/duck/Hong Kong/d73/1976 (HK/d73/76), A/n>n class="Species">chicken/Beijing/1/1994 (BJ/1/94), A/chicken/Korea/38349-p96323/1996 (KR/96323/96), A/chicken/Hong Kong/G9/1997 (G9/97), A/duck/Hong Kong/Y439/1997 (Y439/97), A/quail/Hong Kong/G1/1997 (G1/97) and A/chicken/Shanghai/F/1998 (F/98) (for review see reference [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]. N-linked 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) n class="Gene">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 (n>n class="Disease">FEL), fast unconstrained Bayesian approximation (FUBAR), mixed effects model of evolution (MEME) and single-likelihood ancestor (SLAC) [54-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 (n class="Chemical">NAIs; n>n class="Chemical">oseltamivir, zanamivir, laninamivir and peramivir) was evaluated using a standard fluorescence based assay that measures NA enzymatic activity to determine the 50% inhibitory concentration (IC50) of viruses. This protocol has been described in detail by Leang et al., 2017 [57]. The IC50 values were calculated using the JASPR v1.2 software (CDC, Atlanta, USA). NAI susceptibility is described based on the fold change in IC50 values of the test viruses compared to control viruses; viruses with normal inhibition had <10-fold increase in their IC50, reduced inhibition had between 10 to 100-fold increase in IC50 or highly reduced inhibition >100-fold increase in IC50.

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 pan class="Species">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.

Results

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 n class="Species">chickens and 21 from n>n class="Species">ducks (S1 Table). Further analyses were performed on full genomes of 63 isolates and the partial genome of one isolate (A/duck/Cambodia/Z495W31M3/2015). Phylogenetic analysis of the Cambodian A(H9N2) viruses showed that the diversity of A(H9) AIVs decreased from 2013, when all three main H9-HA lineages (Y439/97, G1/97 and BJ/94) were identified (Fig 1). Only the BJ/94 lineage was detected in 2015 and 2016. Despite arising from the same H9-4.2.5 clade [58], the 2015/16 HA genes (nucleotide identities 95.5–100%) did not cluster with 2013 viruses indicating that HA replacements occurred between the surveillance periods.

Fig 1

Maximum likelihood phylogenetic tree of the A(H9) HA gene.

The tree was produced using IQ-Tree with the GTR+ I + Γ model. Phylogenetic support was estimated using 1,000 ultrafast bootstrap replicates. Bootstrap values greater than 70% are displayed on branches. Cambodian A(H9) viruses are coloured according to the year they were detected: 2013 is shown in orange, 2015 in light blue and 2016 in dark blue. AIVs detected in humans are indicated in red and G57 viruses, as determined by Pu et al., 2015 [18], are shown in pink. Candidate vaccine viruses are indicated by an asterisks (*). Reference viruses are indicated in bold and viral lineages are indicated on the right of the phylogeny. Scale bars indicate the number of nucleotide substitutions per site.

The tree was produced using IQ-Tree with the GTR+ I + Γ model. Phylogenetic support was estimated using 1,000 ultrafast bootstrap replicates. Bootstrap values greater than 70% are displayed on branches. Cambodian A(H9) viruses are coloured according to the year they were detected: 2013 is shown in orange, 2015 in light blue and 2016 in dark blue. pan class="Gene">AIVs detected in n>n class="Species">humans are indicated in red and G57 viruses, as determined by Pu et al., 2015 [18], are shown in pink. Candidate vaccine viruses are indicated by an asterisks (*). Reference viruses are indicated in bold and viral lineages are indicated on the right of the phylogeny. Scale bars indicate the number of nucleotide substitutions per site.

A single virus, designated A/n class="Species">chicken/Cambodia/a77W22M1/2016 (Group A), clustered independently with a lineage predominantly containing A(H9N2) viruses collected from: Vietnamese poultry from 2014 to 2018, a single virus detected in a Japanese n>n class="Species">chicken from 2016 and Indonesian poultry from 2016 to 2018. The remaining viruses (n = 63, Group B) formed a monophyletic lineage with Vietnamese viruses identified from 2014 to 2018 with multiple, independent, region specific clusters indicating frequent cross-border movement of A(H9N2) (Fig 1). All HA genes diverged from a common ancestor around September 2011 with a 95% Highest Posterior Density (HPD) interval from January 2011 to June 2012 (S2A Fig) indicating dispersion from an original Asian source. However, predominant Group B viruses diverged more recently during May 2013 (95% HPD, September 2012—December 2013).

The NA genes of Cambodian A(H9N2) viruses showed a greater genetic diversity (mean TMRCA April 2005, 95% pan class="Disease">HPD August 2004–December 2005) (S2B Fig) with three main clusters (Fig 2), all in the G9/97 NA reference lineage. Group A contained all Cambodian A(H9N2) clade h9-4.2.5 viruses from 2013 (n = 5), and a subset of 2015 (n = 25) and 2016 (n = 1) viruses. Persistence of the N2 lineage containing viruses from 2013, but not associated HA lineages, suggests reassortment between the previously circulating and newly introduced h9-4.2.5 viruses. N2 Group A stemmed from A(H9N2) viruses identified in poultry from Hong Kong n>n class="Species">SAR (China) and mainland China in 2008/09, as well as a single A(H9N2) virus identified in a human from Hong Kong SAR (China) in 2008. Within this cluster, Cambodian A(H9N2) AIVs detected between 2013 to 2016 formed a monophyletic lineage with viruses identified in Vietnamese poultry from 2012 to 2018. N2 Group B contained a single Cambodian isolate, A/chicken/Cambodia/a77W22M1/2016, which clustered similarly to HA Group A with viruses from Indonesia, Japan and Vietnam. Finally, N2 Group C contained some Cambodian A(H9N2) viruses from 2015 (n = 24) and the majority of 2016 isolates (n = 13). Similarly to N2 Group A, this clade stemmed from A(H9N2) viruses identified in China and clustered closely with viruses identified in Vietnamese poultry from 2014 to 2018 (nucleotide identity 89.3–100%).

Fig 2

Maximum likelihood phylogenetic tree of the N2 NA gene.

The tree was produced using IQ-Tree with the GTR+ I + Γ model. Phylogenetic support was estimated using 1,000 ultrafast bootstrap replicates. Bootstrap values greater than 70% are displayed on branches. Cambodian A(H9) viruses are coloured according to the year they were detected: 2013 is shown in orange, 2015 in light blue and 2016 in dark blue. AIVs detected in humans are indicated in red and G57-like viruses, as determined by Pu et al., 2015 [18], are shown in pink. Candidate vaccine viruses are indicated by an asterisks (*) to the right of the taxa name. Reference viruses are indicated in bold and viral lineages are indicated on the right of the phylogeny. Scale bars indicate the numbers or nucleotide substitutions per site.

The tree was produced using IQ-Tree with the GTR+ I + Γ model. Phylogenetic support was estimated using 1,000 ultrafast bootstrap replicates. Bootstrap values greater than 70% are displayed on branches. Cambodian A(H9) viruses are coloured according to the year they were detected: 2013 is shown in orange, 2015 in light blue and 2016 in dark blue. pan class="Gene">AIVs detected in n>n class="Species">humans are indicated in red and G57-like viruses, as determined by Pu et al., 2015 [18], are shown in pink. Candidate vaccine viruses are indicated by an asterisks (*) to the right of the taxa name. Reference viruses are indicated in bold and viral lineages are indicated on the right of the phylogeny. Scale bars indicate the numbers or nucleotide substitutions per site.

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%, pan class="Gene">PB1: 92.9–100%, n>n class="Chemical">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.

Fig 3

Phylogenetic congruency of Cambodian A(H9N2) viruses.

ML trees were produced for all genomic segments using all Cambodian A(H9N2) viruses detected from 2015 to 2016 and key reference viruses. The phylogenies of the Cambodian viruses have been linked across all eight segments and colour coded according to viral genotype, genotype P viruses are shown in purple and genotype V in green. Viral designations are listed on the left hand side of the HA phylogeny. Incongruence is demonstrated by deviations in topology across the segments and can indicate viral reassortment has occurred. Scale bars are shown below each phylogeny and indicate the number of nucleotide substitutions per site.

All Cambodian A(H9N2) gene phylogenies, excluding NA, cluster with lineages containing n class="Disease">zoonotic A(H7N9), A(H9N2) and A(n>n class="Species">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%.

Genotypes of Cambodian A(H9N2) viruses

The 64 Cambodian A(H9N2) viruses isolated during this study belonged to two genotypes, P (n = 28) and V (n = 36). Initially described in Chinese poultry, each genotype stems from multiple reassortment events of A(H9) viruses from the n class="Gene">G1/97, G9/97, BJ/94 and F/98 viral lineages (Fig 4) [1,59]. Genotype P viruses, represented by A/n>n class="Species">chicken/Zibo/B1/2008, were first identified in 2008. They had an F/98 backbone, a G9/97 N2 segment and a G1/97 MP segment [4,59]. Genotype V viruses, represented by A/chicken/Guangxi/SIC19/2014, were first identified in 2014 and are the result of reassortment of genotype S viruses with the G9/97 lineage [1,4]. These viruses have an F/98 backbone with PB2 and MP genes from the G1/97 lineage and a G9/97 NA segment. In Cambodia, genotype P and V were both detected in 2015 and 2016 (Fig 4), and have previously been detected in Vietnam [30]. In Cambodia, the detection rates of each genotype were similar in 2015, P (n = 27) and V (n = 22). However, in 2016, genotype V predominated (n = 14) versus P (n = 1).

Fig 4

Genotypes of A(H9N2) Cambodian viruses detected from 2015 to 2016.

The eight genomic segments are colour coded based on the putative ancestral reference AIVs. The genomic constellations of relevant reference lineage viruses (BJ/94, G1/97, G9/97 and F/98) are shown in the bottom panel. Segments from BJ/94 viruses are shown in black, G1/97 in purple, G9/97 in orange and F/98 in green.

The eight genomic segments are colour coded based on the putative ancestral reference n class="Gene">AIVs. The genomic constellations of relevant reference lineage viruses (BJ/94, n>n class="Gene">G1/97, G9/97 and F/98) are shown in the bottom panel. Segments from BJ/94 viruses are shown in black, G1/97 in purple, G9/97 in orange and F/98 in green.

Cambodian A(H9N2) viruses have molecular markers associated with viral adaptation to mammals

Screens for mutations known to affect n class="Gene">AIV pathogenicity (S3A–S3H Table), are summarised in Table 1 [50]. All Cambodian viruses had HA cleavage site motifs common for A(H9) Ln>n class="Chemical">PAI 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 human-type α2,6 receptors (Table 1 and 32d), including Q226L in 100% of Cambodian A(H9N2) isolates.

Table 1

Summary of amino acid substitutions in Cambodian A(H9N2) viruses associated with an increase in viral fitness.

Protein	Phenotype	Mutation/Motif	Residues for Cambodian Isolates (%)	References
PB2	Increased polymerase activity and replication in mammalian cell lines, increase virulence in mice	A588V	V (2)	[61]
		V598T	T (44)	[62]
		Q591K	K (5)	[63]
		E627K	E (100)	[64,65]
		D701N	D (100)	[65]
PB1	Increased polymerase activity and virulence in mice	D622G	G (100)	[66]
PB1-F2	Truncations to the 90 aa protein increase AIV pathogenicity in chickens	90aa	90aa (95)	[67]
			76aa (2)
			8aa (3)
PA	Increased polymerase activity and replication in mammalian cell lines, increased virulence in mice	V63I	I (100)	[68]
		K356R	R (2)	[69]
PA-X	Full length PA-X increases A(H9) virulence in mice	252aa	252aa (100)	[70]
HA	Multibasic cleavage site can increase viral pathogenicity	Mono-basic	PSKSSR/GLF (86)	[71,72]
			PSRSSR/GLF (14)
	Increased specificity for α2,6 human-type receptors	D101N	N (98)	[73]
		I155T*	T (100)	[74]
		S158N	N (36)	[75]
		E/T190V	V (9)	[76]
		N248D	D (3)
		R496K*	K (100)	[77]
	Increased specificity for α2,6 human-type receptors, enhanced replication in mammalian cells and ferrets, enhanced contact transmission in ferrets	Q226L	L (100)	[78,79]
NP	Increased virulence in chickens	M105V	V (91)	[80]
M1	Increased virulence in mice, chickens and ducks	I43M	M (100)	[81]
	Increased virulence in mice	T139A	A (100)	[82,83]
M2	Increased resistance to amantadine and rimantadine	S31N	N (41)	[84]
NS1	Decreased antiviral response and increased virulence in mice	P42S	S (100)	[85]
	Decreased interferon response and increased virulence in chickens	V149A	A (100)	[86]

* Amino acid (aa) mutations where the phenotype is produced only when in combination with other mutations

Cambodian A(H9N2) viruses contained numerous substitutions associated with increased polymerase activity of n class="Gene">AIVs in n>n class="Species">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 n class="Gene">NS1 were detected that increase n>n class="Gene">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). Compared to available poultry vaccine viruses (A/n class="Species">chicken/Shandong/6/96 (6/96), A/n>n class="Species">chicken/Guangdong/SS/94 (SS/94), F/98), Cambodian A(H9N2) AIV HA was most similar to SS/94 (nucleotide sequence similarity: 90.5% to 91.8%, 6 to 9 antigenic aa differences; S4 Table). Compared to the human candidate vaccine viruses (CVVs) (G9/97, A/Hong Kong/1073/99, A/Hong Kong/33982/2009, A/Hong Kong/308/2014, A/Anhui-Lujiang/39/2018), Cambodian viruses were most similar to A/Hong Kong/308/2014 at a nucleotide level (sequence similarity: 94.1% to 95.1%). Although, they were more similar to G9/97 at H9 HA antigenic sites with 4 to 7 aa differences, compared to A/Hong Kong/308/2014 that had 8 to 11 aa differences.

Post translational modifications

In the Cambodian A(H9N2) viruses, seven conserved HA glycosylation sites were predicted, including: 11NSTE, 280NTTL, 287NVSK/R, 295NCSK and 474NGTY in the HA stalk domain, as well as 123NVSY and 200NRTF in the HA globular domain (S5 Table). More variation was observed in glycosylation sites predicted for N2. Three glycosylation sites (positions: 44, 61 and 146) were conserved in the N2 gene of all Cambodian viruses. Variable sites were identified at positions: 48, 69, 70, 86, 234 and 306. The majority of isolates contained 6 predicted glycosylation sites in NA (n = 40), though isolates with 7 (n = 15), 5 (n = 8) and 4 (n = 1) sites were also identified.

Selection pressure

Selection pressure analysis showed that the majority of sites in each gene of Cambodian A(H9N2) viruses were under neutral or negative selection. A small number of codons in HA (78, 158), NA (9, 12, 416), PB2 (389, 758), M2 (21, 97) and NS2 (119) were under positive selection (S6 Table). Substitutions at the majority of these sites are not associated with increased pan class="Disease">AIV fitness. However, substitution of n>n class="Chemical">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 n class="Mutation">S31N marker associated with resistance to n>n class="Chemical">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 IC50 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) n class="Gene">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, n>n class="CellLine">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/pan class="Species">chicken/Cambodia/a77W22M1/2016, clustered sepan>rately from other Cambodian viruses, possibly representing a sepan>rate A(H9) introduction via migration of wild birds along the Eastern Asia/Australasian flyway. This virus clustered with A(H9N2) n>n class="Gene">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 n class="Gene">AIV lineage detected in Chinese n>n class="Species">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 n class="Disease">zoonotic AIVs such as A(n>n class="Species">H5N6), A(H7N9) and A(H10N8) [18-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) n class="Gene">AIVs were Ln>n class="Chemical">PAI viruses with multiple molecular markers associated with adaptation to mammalian species. Cambodian A(H9N2) viruses contained multiple HA markers associated with increased binding to “human-like” α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 n class="Gene">AIV replication and virulence in mammals (Table 1). While major molecular markers of n>n class="Gene">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 n class="Gene">AIV circulation in poultry flocks and therefore decrease n>n class="Species">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 n class="Disease">zoonotic AIVs as part of n>n class="Species">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) n class="Gene">AIVs isolated from domestic poultry in LBMs between 2015 and 2016 belonged to the n>n class="CellLine">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.

Map of sample collection sites in Cambodia from 2015 to 2016.

Samples were collected from two LBMs during this period, Orussey market (Phnom Penh) and Takeo market (Takeo), indicated by red dots. The map was produced with QGIS v2.18.4 using public domain data obtained from Natural Earth (http://www.naturalearthdata.com/).

(TIF)

Click here for additional data file.

Bayesian maximum clade credibility (MCC) phylogenetic trees of Cambodian A(H9N2) isolates identified from 2013 to 2016 produced from the a) HA and b) NA genomic segments.

Cambodian pan class="Gene">AIVs are coloured based on the year of detection: 2013 is shown in orange, 2015 in green and 2016 in blue. Trees were produced using BEAST v1.84 with GTR + Γ and SRD06 nucleotide substitution model. The nodes are labelled with the proposed TMRCA and node bars show corresponding 95% pan class="Disease">HPD ranges. Branch lengths are time-proportional and the time scale is shown on the x-axis.

(PDF)

Click here for additional data file.

Bayesian skyride analysis of Cambodian A(H9N2) genetic diversity.

Genetic diversity of the HA and NA A(H9N2) genes was estimated using the Gaussian Markov Random Field (GMRF) model. The x-axis measures time in years and the y-axis is an estimate of genetic diversity calculated from Neτ (effective population size and the generation length in years) shown in log scale. The median estimate of genetic diversity over time is shown as a solid black line and the purple shading represents the 95% pan class="Disease">HPD intervals.

(TIF)

Click here for additional data file.

Maximum likelihood phylogenetic trees for all genomic segments of Cambodian A(H9N2) viruses isolated from 2015 and 2016.

a) PB2 b) n class="Gene">PB1 c) n>n class="Chemical">PA d) NP e) MP and f) NS. Trees were generated with IQ-Tree using the GTR+ I + Γ model with 1,000 ultrafast boostrap replicates. Cambodian viruses are coloured based on the year of detection: 2015 is light blue and 2016 dark blue. AIVs identified in humans are coloured red, and G57 lineage viruses (as defined by Pu et al., 2015) are pink. Candidate vaccine viruses are indicated by an asterisks (*) next to the taxa name. Viruses from A(H9) reference lineages are shown in bold and the lineage Cambodian A(H9N2) viruses fall under is indicated on the right hand side of the tree. Bootstrap values of 70 or greater are displayed on branches. The scale bar indicates number of nucleotide substitutions per site.

(PDF)

Click here for additional data file.

Summary of Cambodian A(H9N2) viruses identified in 2015 and 2016 that were analysed as part of this study with sample collection details, genotyping information and sequence accession numbers listed.

(XLSX)

Click here for additional data file.

List of AIVs used to analyse phylogenetic congruence.

(XLSX)

Click here for additional data file.

a-h Molecular analysis of the eight genomic segments for the Cambodian A(H9N2) isolates identified between 2015 and 2016. Data for each segment, as well as any associated accessory proteins is listed: a) PB2, b) pan class="Gene">PB1, c) n>n class="Chemical">PA, d) HA, e) NP, f) NA, g) MP and h) NS

(XLSX)

Click here for additional data file.

Analysis of HA antigenic sites in Cambodian A(H9N2) AIVs compared to human and poultry vaccine viruses.

(XLSX)

Click here for additional data file.

N-glycosylation sites in the HA and NA proteins predicted for Cambodian A(H9N2) isolates identified between 2015 and 2016.

(XLSX)

Click here for additional data file.

Selection pressure analysis of the Cambodian A(H9N2) genes using FEL, FUBAR, MEME and SLAC.

(XLSX)

Click here for additional data file.

Susceptibility of Cambodia A(H9N2) isolates to a panel of four neuraminidase inhibitors: oseltamivir, zanamivir, laninamivir and peramivir.

(XLSX)

Click here for additional data file.