Twelve H6N2 viruses were isolated from various farming systems between September 2015 and February 2019, in three major chicken-producing regions in South Africa (Fig. 1), namely Gauteng (n = 7), North West (n = 4) and Kwa-Zulu Natal (n = 1) Provinces (Table 1). No viruses were isolated from the Western Cape or other provinces during this period, although serological test results suggest that infection is widespread (unpublished national veterinary laboratory records). Viruses isolated in 2015 showed normal patterns of agglutination of chicken red blood cells (RBCs) on the HI tests. Viruses isolated in 2016 showed either weak agglutination or no agglutination on the first passage but dwarfing of the embryos and deaths was observed. Similarly, isolate 432/2019 did not agglutinate RBCs on the first passage, but the allantoic fluid tested positive for IAV on rRT-PCR and H6 conventional RT-PCR (results not shown).
Genome assembly, phylogenetic and reassortment analysis
Between 4.9 and 6.2 million Ion Torrent reads with an average length of 177 base pairs (bp) were generated for each of the 12 sequencing libraries. Maximum likelihood phylogenetic trees (Figs. 2 to 9) were reconstructed for each assembled genome segment to determine the evolutionary relationships with H6N2 viruses isolated from 2002 to 2013 and closest reference sequences retrieved from the Genbank sequence database [22–23]. Genetic clusters A to C were assigned based on topology to evaluate genome reassortment events.
Separation into sub-lineages I and II was evident in the phylogenetic trees for the surface glycoproteins encoded by segment 4 (HA, Fig. 2) and segment 6 (NA, Fig. 3). Each sub-lineage contained the early isolates from 2002 at the root. Progenitor strains A/ostrich/South Africa/KK98/1998 (H6N8) and A/ostrich/South Africa/9508103/1995 (H9N2) were more closely rooted to sub-lineage II in Figs 2 and 3, respectively. None of the viruses isolated from 2015 onward were classified as sub-lineage II.
In segment 1 (PB2 gene; Fig 4), segment 2 (PB1 and PB1-F2 genes, Fig. 5), segment 3 (PA and PA-X genes, Fig. 6), segment 5 (NP gene, Fig. 7) and segment 8 (NS1 and NEP genes, Fig. 9), no separation into sub-lineages I and II was observed, instead grouping was chronological with the late viruses clustering separately from the early strains. Segment 7 that encodes the M1 and M2 genes was the exception where strains separated into sub-lineages I and II as with HA and NA (Fig. 8).
Genome reassortments were assessed by arranging the assigned genetic clusters in Table 2, where the combination of eight segments formed the genotype. Isolates 338087/2015, 344378/2015 and 344579/2015 that cluster closely together in all eight phylogenetic trees (Figs 2 to 9) shared the same genotype (a). The latter two are from the same site. Both locations were in the Gauteng province, but 338087/2015 was isolated from broilers whereas the source of 344378/2015 and 344579/2015 was a layer farm (Table 1). Isolates 341797/2015, 339678/2015, 398997/2016, N2826/2016, 401156/2016 and 402385/2016 also shared the same genotype (b). The first two viruses were isolated from the same site in the North West province. Apart from 402385/2016 and 401156/2016, also from the North West Province but isolated a year after 341797/2015 and 339678/2015, the other viruses in this genotype originated from the Gauteng province. Strains 404573/2016 and 398997/2016 were isolated from the same site, one month apart. Reassortment in segments 6 and 8 of 404573/2016 (genotype c) indicated that this was not a single outbreak, but instead the farm experienced a re-introduction of H6N2. Strain H44954/2016 from KwaZulu-Natal was distinctive; its genotype (d) and position in the phylogenetic trees suggests that the KwaZulu-Natal Province might have variants that are antigenically distinct from the Gauteng and North West Provinces. The most recent isolate, 432/2019 (Gauteng Province) also had a unique genotype (e), its long branches indicated progressive genetic drift from the 2015/2016 isolates, but it shared recent common ancestors in all segments with the Gauteng and North West Province strains apart from segment 4 (HA). Here the HA gene was more closely related to that of H44954/2016 from the KwaZulu-Natal Province (Fig 2).
Genetic drift and antigenic diversity in the HA protein
At the nucleotide sequence level in the HA gene, identities between the H6N2 viruses isolated in 2015/2016 ranged from 99.9% to 93.8%, but the similarity of the 2015/2016 viruses with the single virus isolated in 2019 ranged from 96.6% to as low as 92.6% (Additional file 1: Table S1a). When translated to amino acids which reflect the antigenic variation, identities between 2015/2016 viruses ranged from 100% to 94.8%, whereas the 2015/2016 viruses had a range of 96.7% to 93.8% sequence identity with the single virus from 2019 (Additional file 2: Table S1b).
To further evaluate variation at the amino acid/antigenic level over time, all viruses in sub-lineage I were grouped according to year of isolation, and the distances between groups were calculated in MEGA 5 (Additional file 3: Table S1c). Between the 2002 virus group and the 2012/2013 virus group homology was 96.1%, between the 2002 virus group and the 2015/2016 virus group homology was 95.6% and between the 2002 virus group and the single 2019 virus homology was only 90.5%. Between the 2012/2013 virus group and the 2015/2016 virus group the homology was 96.7%, and over the subsequent 3-year period the sequence homology between the 2012/2013 group and the single 2019 virus had dropped to 95.7%. Amino acid sequence homology between the 2015/2016 group and the 2019 virus was 94.7%.
Molecular markers of host adaptation and virulence across the genome
The HA receptor binding site is surrounded by antigenic sites that are recognised by the most potent neutralizing antibodies. The linear motif RYVRMGTESMN (residues 199 to 205, Additional file 4: Fig. S1, H3 numbering), located on the surface on the globular head of the H6 HA protein, was identified as a highly conserved H6 neutralizing B-cell epitope in monoclonal antibody-binding studies . All sub-lineage I viruses isolated from 2015 onwards contained a V201I mutation in this epitope, and 338087/2015, 344378/2015 and 344579/2015 had a M204L mutation.
Key mutations near the receptor binding site in the HA protein implicated in the switch between H6 viruses binding avian-type receptors to human-type receptors include A138S, I55V, P186L, V187D, E190V, A193N, G225D, Q226L and G228S (H3 numbering; [26–28]). Of the aforementioned, the A138S mutation was absent in the H6 ostrich progenitor but appeared in South African H6N2 sub-lineage I viruses as early as 2002, and have been maintained up until 2019, with the exception of 401156/2016 and 402385/2015 in which an 138L mutation was present (Additional file 4: Fig S1). The V187D mutation was already present in the H6N8 ostrich progenitor virus and has persisted through both sub-lineages I and II to the present time. The A193N mutation was present in all sub-lineage I strains, but absent in the progenitor and sub-lineage II viruses apart from AL25/2002. The remaining six mutations in the H6 HA, I155V, P186L, E190V, G225D, Q226L and G228S were not detected.
The HA0 motif of PQVETRGIF was present in all viruses sequenced here, except for H44954/2016 that acquired an additional E343G mutation in its motif of PQVGTRGIF (Additional file 4: Fig. S1).
N-glycosylation patterns in the HA across all viruses isolated since 2015 were conserved, apart from N2826/2016 that lost N-glycosylation sites at residues 11 and 298 (H3 numbering, Additional file 5: Table S2, Additional file 4: Fig S1). Two predicted O-glycosylation sites were uniform across isolates apart from N2826/2016 that lacked the predicted site at K314.
Like the HA, the NA glycoprotein consists of a globular head and a stalk region. Early studies identified seven antigenic sites on the N2 NA globular head region . Of these, the R344K mutation first emerged in sub-lineage I viruses in 2012/2013 (Additional file 6: Fig. S2) and was also found in four of the twelve viruses from 2015 onwards, namely 338087/2015, 344579/2015, 390997/2016 and 404573/2016.
Most isolates from 2015 onward contained seven predicted N-glycosylation sites across the NA protein (Additional file 5: Table S2; Additional file 4: Fig. S2). As with the HA protein, N2826/2016 was the exception, lacking four sites at residues 200, 309, 313 and 402 (full-length N2 numbering). Five other viruses, 339678/2015, 341797/2015, 398997/2016, 401156/2016 and 402385/2016 also lacked a single the N-glycosylated site at position 309. These five viruses plus 432/2019 also lacked the predicted O-glycosylation site at 351. Isolate 432/2019 lacked the predicted O-glycosylated site at residue 334 that was common to all others, but all viruses contained the predicted O-glycosylation site at 145. H44954/2016 contained two additional O-glycosylation sites at residues 332 and 336.
The neuraminic acid binding site of NA is situated in the globular head, and the residues responsible for the catalytic function are constant for all NA subtypes in IAVs. However, in response to mutations in the HA that decrease its receptor-binding affinity (for example glycosylation), NA is also capable of evolving a second sialic acid binding site, called the HB site, that functionally replaces the HA binding site [30–31]. The HB site is characterised by three loops. One loop is formed by a serine triplet, S367, S370 and S372, a second loop formed from residues 400 to 403 by chain interactions with W403, and K432 is involved in the formation of the third loop from 430 to 433 (Additional file 6: Fig. S2). The serine triplet is present in a conserved region of the N2 sequences of all the isolates we studied. The tryptophan of the second loop is also present in all viruses, but adjacent to this an N-glycosylation site is present in all viruses except for N2826/2016. In the third loop of the HBS, all sub-lineage I strains (but not sub-lineage II) except for H44934/2016 contain a Q432R mutation. Human influenza viruses usually lack the HB site amino acid sequence .
The NA stalk region of sub-lineage I H6N2 viruses is characterised by a 21 amino acid deletion [22–23] which was present in all the viruses analysed here (Additional file 6: Fig. S2), but an additional single codon (serine) was deleted at residue 46 in isolate 432/2019.
The NA protein is a target for antiviral therapy. Six mutations are associated with the acquisition of resistance to neuramindase inhibitors, namely E119V, R152K, D198N, H274Y, R292K and N294S . V96A is specifically associated with resistance to zanamivir and oseltamivir , however none of these mutations were present in the neuraminidase proteins of the South African H6N2 viruses (Additional file 6: Fig. S2).
In the PB2 protein, K73R, A199S, L475M, D567N, E627K, A661T and K702R mutations are associated with IAVs capacity to infect humans [34–36]. All were absent from the South African H6N2 viruses apart from K702R (Additional file 7: Fig. S3). The K702R mutation was present in the H6N8 progenitor virus and in one early sub-lineage II isolate AL25/2002 but was absent from the later three sub-lineage II strains in 2012. Similarly, the mutation was present in the early sub-lineage I viruses AL19/2002 and W–04/2002, but absent from the two strains isolated in 2012/2013. It however was present in all the sub-lineage I viruses post 2015 apart from three Gauteng isolates: 338087/2015, 344378/2015 and 344579/2015. D701N in the PB2 protein was associated with high pathogenicity in mice , this mutation was present in early sub-lineage I strain AL19/2002 but was not detected in any other South African H6N2 strains.
H99Y and I368V mutations in the PB1 protein are associated with airborne transmission among ferrets  but were absent from the H6N2 viruses. Neither of the human-associated markers in the PB1-F2 protein, R79Q and L82S  were detected, nor was N66S that is associated with increased virulence, replication efficiency and antiviral response in mice. No human-associated markers in the PA protein, namely P28L, D55N, V100A, K356R, Q/T/S400L or T552S [34–36] were detected (data not shown). Human-associated markers in the NP protein (V33I, I109V), M2 protein (L55F) NS1 protein (E227R/K) and NS2 protein (S70G) [34–36] were similarly absent from the H6N2 viruses (data not shown).
Comparison of assays and antigens used for serology tests
No clinical signs were observed in any of the Specific Pathogen Free (SPF) chickens throughout the animal experiment. A serial titration experiment with H6N2 antisera was used to demonstrate the relative detection limits of the approved methods used for H6 antibody testing in South Africa, especially where specific antibodies in a serum sample are low. In Table 3, columns represent groups of chickens that were either vaccinated with AVIVAC AI whole virus inactivated H6N2 vaccine (A and B) or not (C), not exposed to field challenge (A) or challenged with H44954/2016 (B and C). Undiluted sera and sera titrated up to 1:256 were tested with two ELISA tests and three HI antigens. The BioChek ELISA detected IAV-specific antibodies in only the undiluted serum of the chicken that was vaccinated with AVIVAC AI but not exposed to field virus (group A). IDEXX ELISA results were negative. H6N2-specific antibodies were detected by HI up to the 1:8 titration using the 2002 H6N2 antigen (homologous with the AVIVAC vaccine) and the 2016 H6N2 antigen (homologous with the field challenge virus), but HI using the H6N8 antigen only detected H6N2 antibodies at the positive threshold cutoff of 4 log2.
In the vaccinated chicken exposed to field virus (group B), both ELISAs detected influenza virus-specific antibodies, BioChek assay up to the 1:16 titration and IDEXX assay up to the 1:64 titration. The 2002 H6N2 antigen was the most sensitive on HI, detecting serotype-specific antibodies up to the 1:64 titration, with the H6N8 antigen detecting antibodies up to the 1:16 titration and the 2016 H6N2 antigen faring slightly better at the 1:32 titration.
The non-vaccinated chicken exposed to the field virus elicited a humoral response that was strongly detectable by the IDEXX ELISA (up to 1:64) and less robustly by the BioChek assay (up to 1:8). Whereas the 2002 H6N2 antigen and the 1998 H6N8 antigen only detected serotype-specific antibodies up to the 1:2 titration and undiluted sera, respectively, the 2016 H6N2 antigen detected H6N2-specific antibodies up to the 1:64 dilution.