Highly pathogenic avian influenza virus of the A/H5N8 subtype, clade 2.3.4.4b, caused outbreaks in Kazakhstan in 2020

Collecting samples for the study

The main reference laboratory in the Kazakhstan’s state system for veterinary control and surveillance, the National Reference Veterinary Center (NRVC), organized expeditions to provinces in Northern and Southern Kazakhstan, where outbreaks of avian influenza were officially registered. The NRVC veterinarians monitored areas within a 10 km-radius from poultry farms or disease-affected private households, as well as around settlements officially designated as “disease reservoirs”. Birds’ carcasses or visibly diseased birds were examined. In total, 2212 birds were examined and used for samples collection. Samples were brains, lungs, trachea, spleen, intestinal tract and feces. The vast majority of the samples were from chickens, domestic ducks and geese. Other samples were from mute swan (Cygnus olor), gray crow (Corvus cornix), dove (Columba livia). Veterinarians diagnosed the disease in live birds based on physical signs such as coma or prostration, unnatural position of head and neck, difficulty breathing, refusal to feed, hyperemia and swelling of visible mucous membranes, flows from beak and nasal openings, diarrhea. When examining carcasses, attention was paid to hemorrhages in the mucous membranes, under serous membranes of internal organs, in the subcutaneous tissue or in the muscles, the presence of enteritis, conjunctivitis or inflammation of other internal organs, pulmonary edema.

Molecular diagnostics and virus subtyping

Frozen samples were transported to NRVC where PCR-diagnostics was performed. Pieces of organs were homogenized using TissueLyser bead mill (Qiagen, Hilden, Germany) to produce 10% suspensions in saline (0.9% NaCl). The suspensions were transferred in 1.5 ml tubes and centrifuged at 10,000 rpm for 30 s. The supernatants were used for RNA extraction. Samples of feces were suspended to obtain 10% suspensions in saline. The suspensions were centrifuged and the supernatants were used for RNA extraction. Total RNA was purified using QIAamp Viral RNA Kits (Cat# 52906, Qiagen, Germany) as per the manufacturer’s instructions.

Viral RNA was detected using a test-system for real-time PCR (RT-PCR), namely PCR-INFLUENZA-A-FACTOR (Cat# R10515-VET, LLC “VET FACTOR”, Troitsk, Russia). This test-system targets for amplification conserved (subtype-independent) sequences in the MP segment in the viral genome. Another kit, PCR-FLU-TYPE-H5/H7/H9-FACTOR (Cat# R12817-VET, LLC “VET FACTOR”, Troitsk, Russia) was used for HA typing. The kit allows assigning HA subtypes H5, H7, H9 based on results of one multiplexed reaction. Thermal cycler Rotor Gene 6000 was used to conduct the assay. Neuraminidase (NA) subtype was determined using sequencing and identification of the closest prototypes.

Whole genome sequencing

Ten samples collected by NRVC veterinarians during their visits to the most disease-affected settlements were selected for sequencing the full AIV genome. Larger utilization of the full-genome sequencing was not possible due to the high costs. The samples were transferred to the National Center for Biotechnology (NCB). Reverse transcription was performed using SuperScript IV Reverse Transcriptase (Cat# 18090010, Invitrogen, CA). The first strand-cDNA was primed using a universal primer MBTuni-12 (5′-ACGCGTGATCAGCAAAAGCAGG-3′). Eight segments in the AIV genome were PCR-amplified using a pair of MBTuni-12 and MBTuni-13 (5′-ACGCGTGATCAGTAGAAACAAGG-3′) (Zhou et al., 2009). These primers target conserved sequences at the 5′- and 3′-termini of genomic segments. The high-fidelity DNA polymerase Phusion (Cat# F-530XL, Invitrogen, CA) was used in PCR. Amplification products were purified using AMPure magnetic particles (Cat# A63880, Beckman Coulter, USA) according to the manufacturer’s instructions.

A fragment library was prepared using Nextera DNA Flex Library Prep Kit (20015829, Illumina, USA). Then, the MiSeq Reagent Kit v3, 600 cycles (MS-102-3003, Illumina, USA) and MiSeq platform (Illumina) were used to conduct the sequencing.

In addition to ten samples sequenced by the authors of this work, six different samples from the year-2020-outbreaks were sequenced in foreign laboratories, namely the Animal and Plant Health Agency (APHA) (Addlestone, United Kingdom) and WHO National Influenza Centre (WHO NIC) in St. Petersburg (Russia). In order to increase the representation of viruses from Kazakhstan in the study, we downloaded the APHA- and WHO NIC-generated sequences from GISAID (GISAID, 2022; Elbe & Buckland-Merrett, 2017; Shu & McCauley, 2017) and included in the analysis in this work. The contribution from the originating laboratories is acknowledged in the Supplementary Information.

Bioinformatics and phylogenetic analysis

Quality of sequencing reads was assessed using the FastQC v0.11.9 program (Babraham Institute, 2019). Trimming to quality values Q30 was performed using the SeqTK v1.3-r117/Sickle suite (Shen et al., 2016). The sequencing reads were assembled into contigs using the BWA tool (Li & Durbin, 2010). The following reference sequences were used as guides during the assembly of contigs: isolate A/White-fronted Goose/AN/1-15-12/2016 (Genbank entries: MH988775.1 for HA, MH989503.1—NA); and isolate A/chicken/Omsk/0112/2020 (GISAID accessions: EPI1813342—PB2, EPI1813343—PB1, EPI1813341—PA, EPI1813338—NP, EPI1813340—MP, EPI1813339—NS). Consensus sequences were produced using FreeBayes  (FreeBayes, 2021) and BCFtools (Danecek et al., 2021). The complete genomes of Kazakhstan’s viruses determined in our laboratory were deposited in GISAID. Our entries are listed in the Supplementary Information.

Sequences of historic prototypes and isolates in contemporary circulation were downloaded from GISAID, their listing with acknowledgements to originating laboratories is provided in the Supplementary Information. In the presented phylogenetic trees, viruses are identified by the GISAID isolate name (EPI_ISL) and accession numbers (EPI#) for particular segments. The goal of this study was to identify isolates in the global circulation which share the closest similarity to viruses from the recent outbreaks in Kazakhstan, and hence, to point to a possible path of introduction. Accordingly, all available complete-genome sequences from Russia, Europe, China and Asia-Pacific collected between Jan. 2020–Aug. 2021 (i.e., in the time frame overlapping the Kazakhstan’s outbreaks) were used in the initial analysis. For Africa and Middle East, sequences deposited from 2018 to 2021 were used because of the shortage of recent (2020-2021) submissions from these regions. The phylogenetic trees were constructed using the Bayesian probabilistic inference with Markov Chain Monte Carlo (MCMC) sampling algorithms in MrBayes v.3.2.7a (Ronquist et al., 2012). The initial analysis was performed with a number of evolutionary process models, i.e., the Jukes-Cantor (JC69) model, generalized time-reversible (GTR) model, and GTR with rate variation across sites (GTR+G) with four rate categories. The initial analysis was done using chain length of 1,000,000. Models showing good convergence such as an estimated sample size (ESS) of at least 1000, potential scale reduction factor (PSRF) close to one (PSRF = 1), and also computation times <24 h, were used in the further analysis. Phylogenetic trees and support statistics for the branches presented in the paper and Supplementary Information were constructed using the MCMC chain length of 10,000,000. With the final settings, ESS values were universally above 10,000 and PSRF = 1 indicating perfect convergence. Consensus trees (50% majority rule trees) are presented in the paper and the Supplementary Information. Posterior probability percents are shown near the nodes to indicate the statistical support for clustering. The H5 HA clades naming and labeling HA tree-branches as the clades 2.3.4.4a–2.3.4.4 h are in accordance to the recent nomenclature (Anonymus, 2022).

FluSurver (A*STAR Bioinformatics Institute, 2022) was used to identify mutations in the influenza virus surface proteins which possibly influence biological characteristics such as host specificity, drugs resistance, virulence, etc.

Bird flu in Northern Kazakhstan, 2020

During autumn 2020, unusually high numbers of wild birds’ carcasses and also live but dying birds were noticed in many places in provinces of Northern Kazakhstan. Local veterinarians reported that initial die-offs started in wild wetland birds which had been migrating southwards. However, exact losses to wild birds’ populations are unknown because no attempts were made to quantify diseased or dead birds in wild habitats. Soon, die-offs in poultry began and attracted attention of national media and officials. Affected poultry were chickens, ducks, geese and turkeys. First officially registered outbreaks in poultry are dated in September 2020. The situation with the disease deteriorated during autumn of 2020, because the disease propagated across the country and increasing numbers of farms faced poultry die-offs. According to data in the official NRVC report, the headcount of diseased birds in private households and commercial farms totaled 69,828 (National Reference Veterinary Center, 2021, unpublished report: http://www.nrcv.kz/). Destruction of all poultry in affected farms was mandatory. The total loss of poultry was 2 million at the end of 2020 including the culling (Sputnik, 2021). Many farms stopped production for about 3 months resulting in an increase in market prices for poultry products (EurasiaNet, 2020). Eleven out of fourteen Kazakhstan provinces had reported outbreaks before the country was declared avian flu-free. Four provinces of Northern Kazakhstan (Akmola, Kostanay, Pavlodar, North-Kazakhstan) had the largest numbers of outbreaks.

The state surveillance system for avian influenza collected samples from 2212 bird’s carcasses. All the samples were tested for AIV using a commercial RT-PCR test system. Total number of PCR-positives amounted to 1976. Ten samples were taken from bird carcasses with overt pathological signs compatible with avian influenza. The samples were collected by NRVC specialists during on-site visits during different expeditions to North Kazakhstan and the Almaty province (Southern Kazakhstan). Map in Fig. 1 shows sampling locations. The ten samples were utilized for full-genome sequencing of Kazakhstan’s AIVs.

Molecular variability of surface glycoproteins

Surface glycoproteins determine antigenic properties of the virus and undergo driving selection under the pressure of the host immune response. Comparing the HA nucleotide sequences within the group of Kazakhstan’s viruses identifies 32 variable positions (per 1704 sequenced positions). Similar comparison of the NA nucleotide sequences shows 36 polymorphic positions (per 1413 nt length). The majority of the nucleotide mutations are silent, as only four polymorphic sites were found at the amino acid level in HA, and 12 polymorphic sites were found in deduced NA amino acid sequences. Data on the amino acid mutations found in surface glycoproteins of Kazakhstan’s viruses are presented in Tables 1–2 for HA, and Tables 3–4 for NA.

The deduced HA amino acid sequences contain the same cluster of basic amino acid residues (KRRKR-G) in all Kazakhstan viruses. This is the cleavage site separating two HA subunits (HA1 and HA2). The presence of such “polybasic cleavage site” is considered to be a major virulence determinant and associated with the HPAI phenotype (Tscherne & García-Sastre, 2011).

The FluSurver service developed by the A*STAR Bioinformatics Institute in Singapore is a tool allowing identifying and annotating mutations in influenza virus proteins which may change virus’ properties such as host range, drug resistance, etc. As compared to the nearest (FluSurver-proposed) reference strain A/Sichuan/26221/2014(H5N6), the HA glycoproteins in Kazakhstan’s viruses share common mutations: K3N, G16S, N110S, T139P, T156A, Q185R, V194I, A201E, N252D, E284G, M285V, I298V, K492E, V538A, I547M and V548M. In addition to the common mutations, two viruses (A/goose/Kazakhstan/4-190-20-B-H5N8-1/2020 and A/chicken/Kazakhstan/220-B-2-H5N8-4/2020) have S12G.

The NA sequences were screened for the presence of biologically relevant mutations using the service FluSurver. As compared to the nearest service-identified prototype A/Baikal teal/Korea Donglim/3/2014(H5N8), Kazakhstan’s viruses share common mutations V8A, V27L, T32M, N46K, V49I, T81A, V106I, S136A, A190T, V201I, V213I, A245S, G263D, R264Q, T265A, T295M, I303V, T329A, V359M, S397L, Y450H, and K469G. Mutations V71I, E72K, N84S, P88S, Q330R, and K469E are present in a small fraction (10–20%) of our collection. Additionally, 40% of our viruses have K141R.

Relation of Kazakhstan’s viruses to isolates in contemporary circulation

The main goal of phylogenetic analysis in this work was to identify isolates in the contemporary circulation which are the closest relatives to Kazakhstan’s viruses. In preliminary experiments we compared different approaches for building phylogenetic trees and calculating the statistical support for branches. The Maximum Likelihood (ML) method with the Bayesian probabilistic selection of the best model-fitting phylogenies, produced the trees with more robust clustering as compared to simple distance-based methods (Minimum Evolution (ME), etc.). Statistical support values for branches in Bayesian-MCMC-phylogenetic trees are the posterior probabilities, and these universally have been higher than bootstrap values commonly used with distance methods. This justifies the use of the Bayesian MCMC approach. The trees presented in Figs. 2 and 3 and figures in the Supplementary Information have been produced using the Bayesian MCMC modeling.

Phylogenetic trees for the segments HA and NA show that all sixteen viruses from Kazakhstan belong to the A/H5N8 subtype. The HA H5 clade 2.3.4.4b comprises all viruses from Kazakhstan (Fig. 2). The H5 clade 2.3.4.4b had evolved at around 2016 and currently this clade predominates worldwide (World Health Organization, 2018; Lee et al., 2017b; Lee et al., 2017a; Chen et al., 2019; Le & Nguyen, 2014; Yamaji et al., 2020; Gu et al., 2013; Tarek et al., 2021; Zhang et al., 2020; Pohlmann et al., 2019; Laleye et al., 2021; Lewis et al., 2021; Globig et al., 2017; Liang et al., 2020; Baek et al., 2021).

In phylogenetic trees (Figs. 2 and 3 and the Supplementary Information) Kazakhstan’s viruses cluster with isolates from Southern Russia (e.g., Rostov-on-Don, Krasnodar), Russian Caucasus (North Ossetia-Alania), Ural region (Chelyabinsk), South-Western Siberia (Omsk). Other closely related prototypes are from Eastern Europe.

In the HA tree, the representatives from Kazakhstan group into five clusters with high statistical support. These clusters are provisionally labeled as branches 1–5 (B1–B5) in Fig. 2. The given labeling is used hereafter to describe the findings. The posterior probability percents for branches B1–B5 are quite high (99–100%). Also, the branch B5 may be divided into lower-order clusters and a single-tip branch. Some of the lower-order branches (labeled B5-1, B5-2, B5-3 in Fig. 2) also have high statistical support (100%). With this regard, the clusterization of the NA sequences follows a very similar pattern. More specifically, branches B1–B4 in the NA tree (Fig. 3) have high statistical support. The latter comprise the same representatives from Kazakhstan as the eponymous branches in the HA tree. Sequences from the cluster B5 in the HA tree appear not to form a separate cluster in the NA tree. However, all studied NA sequences from Kazakhstan do form a mixed group with sequences from South-Western Siberia and Eastern Europe, and this group has significant statistical support (73%).

Some of the clusters defined in the HA and NA trees are actually present in phylogenetic trees for all eight segments. The Kazakhstan’s representatives in the branches B2, B4, B5-2, B5-3 show the same pattern of clusterization in phylogenetic trees for the segments MP, NP, NS, PA, PB1 and PB2 (these trees are presented in the Supplementary Information).

Whereas the described relationships between the HPAI A/H5N8 viruses from Kazakhstan, Siberia and Europe was easy to discern, it has been difficult to discover epidemiological correlates for molecular data such as the dependence between phylogenetic clustering and sampling locations, timing, host species or case pathology. For example, the branches B1 and B2 both comprise viruses from the North and South of Kazakhstan. The branches B3, B4, B5-3 represent viruses from samples collected only in the North, and B4 (the largest cluster of Kazakhstan’s viruses) represents closely related sequences in all trees. In contrast to B5-3, the branch B5-2 members are from the South. Finally, one “stray” variant in the single-tip branch B5-1 is actually from Central Kazakhstan (this sample was sequenced in the APHA laboratory). Nevertheless, the presented data confirm that the HPAI viruses in the circulation involving the territory of Kazakhstan are genetically diverse. The current knowledge on the epidemiological correlates for the Kazakhstan’s HPAI viruses suffers from the real scarcity of both molecular and observational data, and the way to establish such correlations is to enhance the surveillance.