Evolutionary Dynamics of H5 Highly Pathogenic Avian Influenza Viruses (Clade 2.3.4.4B) Circulating in Bulgaria in 2019–2021

The first detection of a Highly Pathogenic Avian Influenza (HPAI) H5N8 virus in Bulgaria dates back to December 2016. Since then, many outbreaks caused by HPAI H5 viruses from clade 2.3.4.4B have been reported in both domestic and wild birds in different regions of the country. In this study, we characterized the complete genome of sixteen H5 viruses collected in Bulgaria between 2019 and 2021. Phylogenetic analyses revealed a persistent circulation of the H5N8 strain for four consecutive years (December 2016–June 2020) and the emergence in 2020 of a novel reassortant H5N2 subtype, likely in a duck farm. Estimation of the time to the most recent common ancestor indicates that this reassortment event may have occurred between May 2019 and January 2020. At the beginning of 2021, Bulgaria experienced a new virus introduction in the poultry sector, namely a HPAI H5N8 that had been circulating in Europe since October 2020. The periodical identification in domestic birds of H5 viruses related to the 2016 epidemic as well as a reassortant strain might indicate undetected circulation of the virus in resident wild birds or in the poultry sector. To avoid the concealed circulation and evolution of viruses, and the risk of emergence of strains with pandemic potential, the implementation of control measures is of utmost importance, particularly in duck farms where birds display no clinical signs.


Introduction
In 2016, HPAI H5 viruses of clade 2.3.4.4B descending from the Gs/GD lineage caused one of the biggest epidemics in Europe by the number of countries affected and the species diversity of infected wild birds [1]. Different HPAI H5 subtypes were detected after the occurrence of reassortment events with LPAI viruses circulating in wild birds. Within the European territory, alongside the HPAI H5N8, which represented the main circulating subtype, HPAI H5N2, H5N5 and H5N6 were also identified [1]. The latter was subsequently responsible for a new epidemic wave in north-eastern Europe between December 2017 and January 2019 [2,3]. In December 2019 a new incursion of a HPAI H5N8 subtype was reported in Poland [4]. Since then, the virus has spread to 8 central-eastern European countries (Bulgaria, Czech Republic, Germany, Hungary, Poland, Romania, Slovakia, Ukraine), affecting mainly the poultry population, with more than 300 outbreaks notified in Hungary and Poland [5][6][7][8].
The last outbreak was reported in June 2020.
In October 2020, with the beginning of the migration of wild birds to their wintering areas, a new HPAI H5 incursion occurred in northern Europe. Similar to the 2016 epidemic wave, the virus had first been detected further east and arrived across a wide area causing variable mortality in the affected wild bird populations and rapidly spread to most of the European countries [9][10][11][12].
In Bulgaria, the first detection of a HPAI H5N8 virus dates back to December 2016. Since then, many outbreaks have been reported in domestic and wild birds in different administrative regions [13]. A recent work by Venkatesh et al. (2020) describes the 2017/2018 HPAI H5N8 epidemic in Bulgaria. The study reports the identification of two separate viral introductions, one in the north-eastern region of Dobrich and another in central and eastern Bulgaria, and suggests that the duck sector played a critical role in HPAI H5N8 virus spread and maintenance in the country [9]. Here, we describe the genetic diversity of HPAI H5 viruses of clade 2.3.4.4B in Bulgaria, analyzing sixteen newly sequenced H5 viruses identified in the country between 2019 and 2021 and the emergence of a new H5N2 subtype, likely from a local reassortment event.

Sample Collection
Between March 2019 and February 2021, sixteen outbreaks caused by the HPAI H5N8 and H5N2 subtypes were registered in chicken and duck farms in Plovdiv, Lovech and Pleven regions ( Table 1).
The National Diagnostic and Research Veterinary Medical Institute of Sofia (Bulgaria) sent clinical samples (intestines, tracheal and cloacal swabs collected from H5Nx infected mule ducks, layer hens and domestic ducks) as well as infectious allantoic fluids (harvested after inoculating 9-11 day old embryonated specific pathogen free (SPF) hen's eggs for virus isolation) to the EU/OIE/National Reference Laboratory for Avian Influenza and Newcastle Disease at the Istituto Zooprofilattico Sperimentale delle Venezie in Legnaro (Italy) for genetic characterization.

Results
In this study, we analyzed the genetic characteristics of the HPAI H5N8 and Czech Republic and Slovakia during the first half of 2020, the outbreaks reported in Bulgaria in the same time period were not caused by the strain newly introduced in Europe through wintering migrations of wild birds (Figure 1). viruses which caused sixteen outbreaks notified in Bulgaria since the beginning of 2019, and specifically in the periods March-April 2019, February-June 2020 and February 2021. The phylogenetic analysis of the hemagglutinin (HA) gene showed that the HPAI H5N8 viruses collected in March-April 2019 (outbreaks 19/1-19/5) and the H5N8 (outbreaks 20/2, 20/5, 20/7) and H5N2 (outbreaks 20/3, 20/6, 20/8) viruses collected in the period February-June 2020 are related to HPAI H5N8 viruses of clade 2.3.4.4B, which have been identified in Bulgaria since the end of 2016. Despite the introduction in December 2019 in east-central Europe of a novel HPAI H5N8, which extensively circulated in Poland, Germany, Hungary, Czech Republic and Slovakia during the first half of 2020, the outbreaks reported in Bulgaria in the same time period were not caused by the strain newly introduced in Europe through wintering migrations of wild birds (Figure 1).  The first outbreak of 2020 (20/1) occurred in a mule duck fattening farm. The clinical sample and the respective isolate were both confirmed to be co-infected with H5N2 and H5N8 viral subtypes by Real Time PCR (using specific primers both for the N2 and for the N8 subtypes), as well as sequencing. Specifically, next generation sequencing (Illumina) of the clinical sample only determined an H5N2 subtype, which was the prevalent one, while a partial sequence of about 300 nucleotides of the N8 segment was obtained using Sanger sequencing. Differently, next generation sequencing of the isolate clearly showed a mixed H5N2/H5N8 virus population and the complete sequence of both N2 and N8 genes was obtained. The analyses of the complete genomes confirm that the 2019-2020 H5N8 belong to the genotype responsible of the 2016 European epidemic (Figure 2), clustering with HPAI H5N8 viruses circulating in Bulgaria in 2017-2018, for all the eight gene segments. Specifically, they can be grouped into three different sub-clusters: the first one includes the HPAI H5N8 viruses collected in 2019 (outbreaks 19/1 to 19/5); the second one is composed of HPAI H5N8 viruses from outbreaks 20/2 (ducks) and 20/5 (mule ducks) from Plovdiv region; the third one contains HPAI H5N8 viruses from outbreaks 20/7 (hens) and 20/9 (hens), respectively from Kurdzhali and Plovdiv regions.
sample and the respective isolate were both confirmed to be co-infected with H5N2 and H5N8 viral subtypes by Real Time PCR (using specific primers both for the N2 and for the N8 subtypes), as well as sequencing. Specifically, next generation sequencing (Illumina) of the clinical sample only determined an H5N2 subtype, which was the prevalent one, while a partial sequence of about 300 nucleotides of the N8 segment was obtained using Sanger sequencing. Differently, next generation sequencing of the isolate clearly showed a mixed H5N2/H5N8 virus population and the complete sequence of both N2 and N8 genes was obtained.
The analyses of the complete genomes confirm that the 2019-2020 H5N8 belong to the genotype responsible of the 2016 European epidemic (Figure 2), clustering with HPAI H5N8 viruses circulating in Bulgaria in 2017-2018, for all the eight gene segments. Specifically, they can be grouped into three different sub-clusters: the first one includes the HPAI H5N8 viruses collected in 2019 (outbreaks 19/1 to 19/5); the second one is composed of HPAI H5N8 viruses from outbreaks 20/2 (ducks) and 20/5 (mule ducks) from Plovdiv region; the third one contains HPAI H5N8 viruses from outbreaks 20/7 (hens) and 20/9 (hens), respectively from Kurdzhali and Plovdiv regions.  The H5N2 viruses, all identified in Plovdiv region from outbreaks 20/1 (H5N2/N8 coinfected mule ducks), 20/3 (hens), 20/6 (hens) and 20/8 (hens), derive from reassortment events between the HPAI H5N8 and LPAI viruses. Specifically, the HA and M genes (Figures 1 and S7) cluster with HPAI H5N8 viruses circulating in Bulgaria in 2017-2020, showing the highest identity (99.1-99.3% for the HA gene) with viruses of sub-cluster 3 (Figures 1 and S1-S8); while the N2 gene clusters with a LPAI H4N2 virus identified in ducks in Bulgaria in 2018 (95.5-95.6% identity) (Supplementary Figure S5). The internal gene segments of the H5N2 subtype show the highest similarity with different LPAI viruses identified in wild and domestic birds in Eurasia (Supplementary Figures S1-S4 and S8). In  Figure S8). Of note, a sample co-infected with H5N8 and H5N2 viral subtypes had already been identified in ducks in Plovdiv region in 2018 (accession numbers EPI1743963-EPI1743964, GISAID EpiFlu database) but its N2 gene is very different from that of the 2020 H5N2 viruses (Supplementary Figure S5), indicating an independent reassortment event.
To estimate the time of emergence of the HPAI H5N8 and H5N2 viruses in Bulgaria we performed a Bayesian evolutionary analysis of all the eight gene segments, and we estimated the time to the most common recent ancestor (tMRCA). The HPAI H5N8 viruses belonging to the cluster of the 2018-2020 Bulgarian viruses seem to have emerged in the first half of 2017 (December 2016-August 2017) (Figures 3 and S9-S18). Conversely, the reassortant H5N2 viruses collected in Bulgaria in 2020 appear to have been generated between May 2019 and January 2020, with the exception of NS gene segment, for which the limited number of allele B viruses widens the tMRCA estimate to between October 2017 and July 2019 (Figures 3 and S9-S18). The 2020 reassortant H5N2 strain was identified only in duck and chicken farms located in Plovdiv region, suggesting that it likely emerged in this area, which was hardly affected by the H5N8 subtype outbreaks.
In February 2021 a new introduction of HPAI H5N8 viruses of clade 2.3.4.4B was recorded in the northern region of Pleven. The 2021 HPAI H5N8 viruses (outbreaks 21/1-21/3, Table 1) cluster separately from the viruses previously detected in the country and show the highest identity with the HPAI H5N8 viruses, which have been circulating in Europe since October 2020 [12]; in particular the three Bulgarian viruses are almost identical and cluster with HPAI H5N8 viruses identified in Poland, England and Russia (Krasnodar Region) in November 2020 and January 2021 (Figures 1 and 2).
Molecular analysis of the HA protein revealed that the Bulgarian H5N2 viruses obtained from outbreaks 20/1 and 20/3 possess a potential additional glycosylation site at position 70-72 (D70N) (numbering starting from the first methionine), which likely indicates an adaptation of the virus to domestic birds [24]. For the NA protein, potential additional glycosylation sites were identified in A/chicken/Bulgaria/217_20VIR1724-1/2020 (H5N8) from outbreak 20/7 at position 411-413 (E411N). No evidence of known mutations associated with mammalian adaptation were recorded in any of the analysed viruses.   Table 1) cluster separately from the viruses previously detected in the country and show the highest identity with the HPAI H5N8 viruses, which have been circulating in Europe since October 2020 [12]; in particular the three Bulgarian viruses are almost identical and cluster with HPAI H5N8 viruses identified in Poland, England and Russia (Krasnodar Region) in November 2020 and January 2021 (Figures 1 and 2).
Molecular analysis of the HA protein revealed that the Bulgarian H5N2 viruses obtained from outbreaks 20/1 and 20/3 possess a potential additional glycosylation site at position 70-72 (D70N) (numbering starting from the first methionine), which likely indi-  Figure 1 and Figures S1-S8) and the H5N2 (light green box in Figure 1 and Figures S1-S8) subgroups.

Discussion
In Europe, the last cases of HPAI H5N8 viruses related to the 2016-2017 epidemic were detected in March and July 2018 in Italy and Russia, respectively [25,26]. After the detection of the HPAI H5N8 in 2016, Bulgaria has continued to identify up to June 2020 H5 viruses related to the 2016-2017 European epidemic in domestic birds. This may indicate an undetected circulation of the virus in resident wild birds, or more likely in poultry sector, in particular in duck farms where the HPAI H5 viruses have shown to be able to circulate causing no or few clinical signs [27,28]. Extremely low biosecurity measures in foie gras farms and persistent circulation of a single strain in this poultry sector have been previously reported in the country [29]. Several studies based on epidemiological and phylogenetic analyses have recognized domestic Anseriformes to play a critical role in HPAI virus persistence, spread and evolution as well as being an important source for virus infection in chickens [30][31][32][33][34].
After France, Bulgaria is one of the main European foie gras producers, with a highly developed commercial poultry sector which includes mule ducks and hybrid ducks raised for foie gras and meat production. Following recovery from the 2016-2017 avian influenza crisis, production underwent a sharp increase in 2018. In 2019, the Bulgarian industry produced about 2600 tonnes of foie gras intended for export mainly to France (71%), Belgium (14%) and Spain (14%) (https://www.itavi.asso.fr/download/10930, accessed on 10 October 2021, Situation de la production et du marché du foie gras en avril 2020, accessed on 20 September 2021). This rapid expansion, not always accompanied by the implementation of risk-appropriate biosecurity measures with only very basic sanitation and hygiene, could represent the ideal environment for avian influenza viruses to be introduced, maintained undetected and reassort, as demonstrated by the identification of two independently originated H5N2 reassortant strains in 2018 and 2020 in one of the regions with the highest duck population density [35].
Of note, the 2020 H5N2 was firstly identified in a farm rearing mule duck, where the H5N8 was also detected. Additional testing would have been useful to elucidate the extent of circulation of the two subtypes in the farm. The subsequent identification of a genetically related H5N2 strain in hens (outbreaks 20/3, 20/6 and 20/8) may indicate a virus transmission between the two poultry sectors, which may have been fostered by the use of vehicles interchangeably employed to transport eggs, fodder and animals from/to different production holdings. However, we cannot exclude that the virus might have circulated undetected also in other poultry sectors (i.e., backyards). The estimation of tMRCA indicated that this subtype may have been generated between the end of 2019 and the beginning of 2020, which clearly indicates that this genotype was detected soon after its emergence, about two years after the introduction in Bulgaria of the H5N8 progenitor strain (first half of 2017 based on tMRCA estimation). However, given the lack of identification of the LPAI progenitor virus, it is difficult to determine whether the reassortment event took place within the mule ducks farm or in another setting (i.e., in another farm or in wild birds). Increasing active surveillance in wild birds may help to shed light on the ecological context which may favor the emergence of reassortant viruses.
The new introduction in February 2021 of a novel H5N8 strain, which has been showing a high propensity to reassort and the ability to acquire mutations associated with mammalian adaptation, is of concern [36]. The implementation of adequate control measures as well as a regular and capillary surveillance in the poultry industry, in particular in the duck farms, also in absence of evident clinical signs is of utmost importance to stop the circulation of the virus in the country. This could help preventing the emergence of novel reassortant viruses, or of variants with unexpected biological properties, which may represent a serious risk for the operators in the poultry sectors, but also for the general public health.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/v13102086/s1, Figure S1: Maximum Likelihood phylogenetic tree of the PB2 gene, Figure S2: Maximum Likelihood phylogenetic tree of the PB1 gene, Figure S3: Maximum Likelihood phylogenetic tree of the PA gene, Figure S4: Maximum Likelihood phylogenetic tree of the NP gene, Figure S5: Maximum Likelihood phylogenetic tree of the NA (N2 subtype) gene, Figure S6: Maximum Likelihood phylogenetic tree of the NA (N8 subtype) gene, Figure S7: Maximum Likelihood phylogenetic tree of the M gene, Figure S8: Maximum Likelihood phylogenetic tree of the NS gene, Figure S9: Maximum Clade Credibility Tree of the PB2 gene, Figure S10: Maximum Clade Credibility Tree of the PB1 gene, Figure S11: Maximum Clade Credibility Tree of the PA gene, Figure S12: Maximum Clade Credibility Tree of the HA gene, Figure S13: Maximum Clade Credibility Tree of the NP gene, Figure S14: Maximum Clade Credibility Tree of the NA-N2 subtype gene, Figure S15

Data Availability Statement:
The consensus sequences of the viruses analysed in this study were submitted to the GISAID EpiFlu™ database under the accession numbers reported in Table 1.