Epidemiology and Evolutionary Dynamics of High Pathogenicity Avian Influenza (HPAI) H5N1 in Bangladesh

. Understanding disease clustering and transmission patterns improves the prevention and control of disease. Herein


Introduction
Since the virus was frst discovered in 1996 in southern China, the global spread of High Pathogenicity Avian Infuenza (HPAI) H5N1 in poultry and wild birds has posed a signifcant panzootic risk and a concern to public health [1]. Its spread outside Asia and the emergence of the HPAI H5N1 virus in Europe in the second half of 2005 underlined the magnitude of the pandemic threat [2]. In 2020, a new clade 2.3.4.4b of H5N1 emerged [3], increasingly spreading through the migratory birds and causing outbreaks in many parts of Africa, Asia, and Europe [4,5]. By April 2023, the virus was discovered in a wide spatial area that includes more than 100 countries and islands across the world, with 8344 outbreaks in domestic and 5573 outbreaks in wild birds [6,7]. In addition, there has been a spread to nonavian species [3]. Since its emergence, the virus has resulted in at least 873 laboratory-confrmed human cases in 21 countries [8].
HPAI H5N1 virus emerged in Bangladesh in 2007. Since then, 585 outbreaks have been recorded in 54 of Bangladesh's 64 districts, making it one of the countries with the highest number of reported cases worldwide [6,9], challenging both the country's commercial and backyard poultry production systems [10][11][12]. Fortunately, there have been just eight human cases of H5N1 in Bangladesh to date. After the implementation of vaccination of poultry against HPAI H5N1 in Bangladesh in 2012, the virus continued to cause sporadic outbreaks in poultry in various locations of the country [13,14]. In 2008 alone, due to HPAI H5N1 outbreaks in Bangladesh, more than 1.8 million chickens were culled, resulting in an estimated economic loss of US$ 40 million [9,15].
Consequently, the control of HPAI in poultry is of paramount importance for Bangladesh. To this end, it is crucial to understand the dynamics of H5N1 outbreaks for future planning, decision-making, and acting. Studying spatiotemporal dynamics can, for instance, identify highrisk regions where monitoring and early detection eforts should be concentrated [16] and support the development of mitigation strategies such as culling or vaccination campaigns [17,18]. Te information is not only of relevance to Bangladesh, but also in building an understanding of spatiotemporal dynamics globally and how the virus is spreading across national boundaries and continents [19,20]. Investigations into the spatiotemporal dynamics of HPAI H5N1 outbreaks may also include identifying disease clusters, which can also lead to better disease control and prevention eforts [21]. Te space-time pattern of HPAI H5N1 outbreak clustering in Bangladesh could designate places where more research should be conducted to learn about the virus's entrenchment and propagation [22][23][24]. In zoonotic disease outbreaks worldwide, including HPAI H5N1, the signifcance of using space-time clusters has been demonstrated previously [23,[25][26][27]. Troughout Asia, signifcant space-time efects on the spread and severity of HPAI H5N1 outbreaks were found in China [28], South Korea [29], Vietnam [30][31][32], and Tailand [33]. Terefore, we conducted an analysis of the spatiotemporal dynamics of HPAI H5N1 outbreaks, including a range of geospatial tools to identify clusters and hotspots of HPAI H5N1 in Bangladesh. While several studies have previously reported on the epidemiology of HPAI H5N1 in Bangladesh, these were limited in the extent to which spatiotemporal dynamics were investigated and were also limited to the period prior to the onset of poultry vaccination against HPAI H5N1 in 2012 [10,11], which was a paramount initiative to reduce the burden of HPAI H5N1 on Bangladesh's poultry industry.
Also, our understanding of the evolutionary dynamics of HPAI H5N1 in Bangladesh is still limited. What [36][37][38]. Here, we combined these previous works and present a holistic overview of the phylodynamic of HPAI H5N1 sequences collected from the domestic chickens, waterfowl, and other wild birds in Bangladesh between 2007 and 2020 and explore the genetic diversity of HPAI H5N1 over time and across species in Bangladesh.  [39]. Most of the HPAI outbreaks reported in Bangladesh have occurred in poultry farms [10]. HPAI outbreak reporting in Bangladesh relies on passive and active reporting. We tabulated all 585 reported HPAI outbreaks in Bangladesh from the Global Animal Disease Information System (EMPRES-i) database between 2007 and 2020, available at the Food and Agriculture Organization (FAO), which also includes data from the World Organization for Animal Health (WOAH) websites, between 2007-2020 [6,9]. Te EMPRES-i data was accessed on January 31, 2022. No outbreak data have been ofcially reported to FAO and WOAH by Bangladeshi authorities after 2018. FAO reports outbreak data, including subtype, pathogenetic status, detection date, and geographical location, from which we determined the district (n � 64) in Bangladesh in which the outbreak took place using digital maps from the open-source DIVA-GIS website [40].

H5N1 Sequence Data.
On May 23, 2022, we retrieved all accessible HA gene sequences of H5N1 with respective epidemiological metadata (n � 413) from the GISAID Epifu database [41]. We also accessed the national center for biotechnology information (NCBI) infuenza virus database [42] for cross-checking but found no additional HA gene sequences for Bangladesh. For all sequences, we collected corresponding metadata from GISAID supplemented with other published information from published articles and other publicly available sources. Tese data included geospatial data and information on whether the host was from the wild or domestic and the type of production system (backyard or commercial) from which it came. Tese data were used to supplement the outbreak data from EMPRES-i. For the building of phylogenetic trees, all sequences from 2007 to 2020 were downloaded.

Description of the Temporal and Spatial Distribution of HPAI H5N1
Incidences. Te monthly number of HPAI H5N1 outbreaks across Bangladesh, as well as the number of afected districts across the nation and incidences based on metadata associated with published sequences, was presented as a time plot to show the temporal dynamics of the HPAI H5N1 outbreaks. Te distribution of outbreaks across species and production systems (i.e., backyard or commercial farm) was demonstrated using prevalence, including a 95% binomial exact confdence interval of the prevalence [43]. During epizootic waves, the frequency of disease outbreaks rises quickly and then progressively falls towards the end of the epizootic [17]. To detect HPAI H5N1 outbreak epizootic waves, we considered the endpoint for the epizootic curve to be reached if there were no reported outbreaks in two consecutive months [44].
For spatial presentations and data analyses, we consid- Adaptive kernel density estimation is a method for estimating the frst-order features of a spatial stochastic process refecting a global or large-scale trend [45], which we used to create a smoothed presentation of how the distribution and frequency of outbreaks varied across the nation. Te kernel density is estimated as follows. Let (x 1 , x 2 , . . . , x n ) be independent and identically distributed samples drawn from some univariate distribution with an unknown density f at any given point x. We are interested in estimating the shape of this function f. Its kernel density estimator is as follows: where K is the kernel, a non-negative function, and h > 0 is a smoothing parameter called the bandwidth.

Statistical Analyses of Global Spatial Patterns in HPAI H5N1
Outbreaks. We used Global Moran's I index [24], Global Geary's C [46], and Global Getis-Ord Gi * [24] to analyze the spatial patterns in HPAI H5N1 outbreaks. A Moran's I value close to zero indicates a random distribution, while a negative Moran's I value indicates more dispersed H5N1 outbreak sites from random and positive values which indicate a clustering of H5N1 outbreaks. Moran's I value is calculated using the following equation: where x i and x j are the deviations of an attribute for feature i and j from their respective means, w i,j is the spatial weight between features i and j, n is equal to the total number of features, and S o is the aggregation of all spatial weights.
Geary's C measures spatial autocorrelation with values between 0 and 1 indicating positive and values higher than 1 indicating negative autocorrelation is calculated as follows: where N is the number of spatial units indexed by i and j, x is the variable of interest, x is the mean of x, w ij is a matrix of spatial weights with zeroes on the diagonal (i.e., w ii � 0), and W is the sum of all w ij . Finally, Global Getis-Ord Gi * values that are negative identify cold spots and values that are positive identify hot spots or areas where there are particularly many outbreaks. Using the parameters used in the calculation presented above, this statistic is calculated using the following equation: (4)

Statistical Analyses of Local Spatial Patterns in HPAI H5N1
Outbreaks. Aside from global spatial statistics, local spatial statistics were also calculated to identify specifc areas with particularly high or low numbers of outbreaks. Te local Moran's I statistic is calculated as follows: where x i is the attribute feature of i, x is the mean of the corresponding attribute, ω i,j is the spatial weight between features i and j, with n being the total number of functions, and Te value of local Moran's I range from +1, indicating high-high (HH) or low-low (LL) clusters, through 0 (� a random pattern) to −1, indicating high-low (HL) or lowhigh (LH) outliers. Tus, a HH (LL) cluster denotes several neighboring areas with a relatively high (low) incidence of outbreaks, while a HL (LH) outlier is a high (low) value surrounded by a region with a predominantly low (high) outbreak frequency.
Te local Getis-Ord Gi * statistic is calculated using the following equation: where x j is the attribute value for feature j, w i,j is the spatial weight between feature i and j, n is equal to the total number of features and the interpretation of the values is similar to that of the Global Getis-Ord Gi * statistic.  2011-2014, and 2015-2018. In this analysis, we largely followed the procedure as outlined in [27] using a temporal scanning window of a maximum of 2 years within each four-year period and a spatial scanning window that allowed a maximum of 50% of the outbreaks.

Bayesian Phylogenetic Analysis.
We used MEGA 11 to align the sequences (https://www.megasoftware.net) and identify sequences with artifacts that were eliminated. Ultimately, 274 sequences were retained for analyses. We constructed the Bayesian phylogenetic tree across all sequences using BEAST v.1.10.4 (https://beast.community) with an uncorrelated lognormal relaxed molecular clock. Te gamma-distributed rate variation among sites with four rate categories (HKYþG) [47] was used. Te Markov Chain Monte Carlo (MCMC) was run for 50 million steps, and the trees and the parameters were sampled every 5,000 steps. Tree Annotator (https://beast.community/treeannotator) was used to generate time-scaled maximum clade credibility (MCC) tree in BEAST and was visualized using Figtree 1.4.3 (http://tree.bio.ed.ac.uk/software/fgtree/). We also constructed a phylogenetic tree for 289 sequences from 2011-2020 belonging to clade 2.3.2.1a following the same methodology as described above for all sequences combined.

Description of Temporal and Spatial Distribution of HPAI H5N1 Incidences.
Tere have been 585 H5N1 outbreaks in poultry and wild birds since the frst incidence was identifed in 2007. Te vast majority, i.e., 98.8% (95% CI: 97.6-99.5) of the outbreaks, concerned chickens, with only 1.2% (95% CI: 0.5-2.5) of the outbreaks in house crows (Corvus splendens). Furthermore, most outbreaks in chickens were reported among commercial chickens (89%), followed by backyard chickens (9.7%) ( Table 1). Te temporal dynamics of HPAI H5N1 outbreaks in Bangladesh (Figure 1), show two major peaks in February 2008 and January 2011, with 58 and 49 outbreaks, respectively. After vaccination against HPAI H5N1 was introduced at the start of 2012, the number of cases dropped dramatically. Although the ofcial reporting of outbreaks stopped in 2018, H5N1 cases were still detected in 2019 and 2020 based on metadata related to published sequences (green line in Figure 1). Nine waves of HPAI H5N1 infections can be decerned between 2007 and 2020. Table 2 enumerates all waves, along with the number of documented outbreaks. All waves were centered around wintertime in Bangladesh. Te frst 5 waves were the largest, concerned chickens exclusively, and occurred annually between the frst occurrence of HPAI H5N1 in Bangladesh and the moment vaccination against HPAI H5N1 was introduced. Figure 2(a) demonstrates the temporal distribution of these 5 epizootic waves in more detail, while Figure 2(b) depicts the spatial scale as the frequency of outbreaks across the impacted districts. Te latter shows Dhaka to be the most afected district by wave 1 and 2, followed by Chittagong, Cox's Bazar, and Narayanganj in waves 3, 4, and 5. Te remaining four small waves, while not occurring annually, remained centered around wintertime and foremostly concerned outbreaks amongst house crows.
Te spatial distribution of the 585 HPAI H5N1 outbreaks in Bangladesh and how it changes over time are presented in more detail in Figure 3. Te maps depict that H5N1 outbreaks concentrated in and around Dhaka district, also when the number of outbreaks decreased over time throughout the country.
Te adaptive kernel density estimation for the HPAI H5N1 outbreaks in each of the three-time periods (Figure 4) changes little in locations over time, with a major concentration around the Dhaka district. Tis confrms the impression already provided by Figure 3, including the decrease in intensity of the kernel density in the 2015 to 2018 period. Overall, outbreaks occurred along an oblique line that ran from southeast to northwest through the country.

Statistical Analyses of Spatial Patterns in HPAI H5N1
Outbreaks. Confrming the impression of this nonrandom distribution and clustering of outbreaks across the country, global Moran's I was signifcant and positive across the entire period 2007-2018 as well as the four-year interval ( Table 3). In line with this, global Geary's C values tended to be lower than 1 across most periods, suggesting spatial autocorrelation in the outbreak, but in none of the cases this was signifcant (Table 3). Finally, again in line with the abovementioned, global Getis-Ord Gi * values are signifcantly positive, suggesting the existence of outbreak hot spots for all investigated periods (Table 3).
Te abovementioned signifcant Global Moran's I and Getis-Ord Gi * spatial statistics call for a more in-depth, local analysis identifying the districts that stand out and ultimately cause the global, nonrandom distribution of outbreaks. In doing so, we found that from 2007 to 2010, the local Moran's I had greater values in Dhaka, Gazipur, Manikganj, Narayanganj, and Narsingdi. From 2011 to 2014, a similar pattern was observed, with a greater value in Dhaka, while Madaripur was added to the cluster. Finally, from 2014 to 2018, in Gazipur, the signifcance of the H5N1 cluster decreased and the local Moran's I values decreased dramatically ( Figure 5).

Phylodynamic and Clade Diversity of H5N1 Virus in
Bangladesh. Following a comprehensive phylogenetic analysis of HA gene sequences of H5N1 HPAI virus isolates from Bangladesh, four diferent clades of H5N1 viruses were discovered (Figures 8 and 9 Te time-scaled phylogenetic tree of HA genes indicated that H5N1 clade 2.3.2.1a viruses in Bangladeshi avian hosts (domestic Anseriformes, domestic Galliformes, and wild birds) diverged into at least nine genetic subgroups (R1-R9), as supported by a high posterior probability (>99%) in the MCC tree ( Figure 10). After 2016, most subgroups (R1-R8) were not observed and replaced by viruses belonging to the R9 subgroup.

Discussion
We investigated the spatiotemporal patterns and evolutionary dynamics of HPAI H5N1 outbreaks in poultry and wild birds in Bangladesh between 2007 and 2020. Most outbreaks occurred in the January-March period when Bangladesh's ambient temperatures are typically at their lowest. Tis broadly corresponds with the period also   designated as high-risk in other Southeast Asian countries [30,48,49] and might relate to higher survival of avian infuenza viruses in colder environments when outside a host [50,51]. However, our phylodynamic analyses clearly indicate that HPAI H5N1 viruses continued circulating during the summer periods and can be considered endemic in Bangladesh. Te nine identifed waves difered in length, the timing of peak incidence and number of outbreaks involved. Aside from these wave variations being possibly due to variations in surveillance intensity and sensitivity, they also may have been due to diferences in response to each outbreak wave and diferences in disease dynamics [15]. Regarding the latter, wave 5 in 2011 coincided with the introduction of a novel H5 clade 2.3.2.1 in Bangladesh. Te advent of this virus may have changed transmission dynamics, such as the severity of infection as a percentage of fock size, in contrast to viruses prevalent in 2008 [44]. While vaccination of chickens in commercial poultry farms, which commenced in early 2012 towards the end of the ffth wave, resulted in a dramatic decrease in HPAI H5N1 outbreaks across Bangladesh, based on published sequence data, we concluded that HPAI H5N1 continued to be present all over the country. Tis is yet another reason to consider HPAI H5N1 to be endemic in Bangladesh. Although most outbreaks occurred in poultry and mainly in chicken, our study revealed that approximately 1% of outbreaks occurred in wild birds, especially house crows. Tat wild birds may get infected with avian infuenza   but subsequently recover has previously been concluded based on the occurrence of antibodies against avian infuenza in blood samples of wild birds in Bangladesh [52]. Waves 7-9 were dominated by house crow mortalities, and these could be connected to live bird markets (LBMs). Several subtypes of avian infuenza virus, including HPAI viruses, are circulating in LBMs, and scavenging house crows may get infected with H5N1 during feeding on ofal at LBMs [53][54][55].

Geospatial Modeling of H5N1 Outbreaks in Bangladesh.
We used exploratory and analytical geospatial statistics to characterize the spatiotemporal pattern in HPAI H5N1 outbreaks in chicken and wild birds in Bangladesh from 2007 to 2018. Tese analyses showed that outbreaks were centered in and around the Dhaka district, supporting observations from the previous studies [10,11]. Dhaka has the highest human population density in the country. Human activity can substantially infuence the spread of infectious diseases [15,56], with human population density typically being associated with high poultry production and trading, increasing cross-infection risk [57].
Tere is strong demand for chicken products in the metropolitan city of Dhaka, where poultry populations from diferent parts of the country are congregated in LBMs [58], which are known to act as hubs for AIV transmission [53]. Road networks have previously been considered to be important in explaining the spread of HPAI H5N1 [33,59,60]. Using space-time permutation modeling, we found clusters in the Dhaka, Bogura, and Barisal region from 2007 to 2010, with large cities (Ramu, Cox's Bazar, and Araihazar, Narayanganj) apparently increasing the risk of viral propagation in the surrounding areas [27]. Te intensity of H5N1 outbreak occurrence was highest in 2007-2010, followed by 2011-2014, and then the lowest intensity was observed in 2015-2018, due to vaccination programs against the H5N1 virus in commercial poultry starting in 2012 [61].   [74].

Evolution and Transmission Dynamics of H5N1 Subtypes and Teir Clade
Te new reassorted 2.3.2.1a H5N1 viruses present in group R9 contain the hemagglutinin, neuraminidase, and matrix genes of the ancestral R8 group and other internal genes of LPAI Eurasian-lineage of AIV that are closely related to LPAI viruses isolated from nomadic domestic ducks and wild birds in a wetland region of northeastern Bangladesh. In this area, domestic ducks have frequent contact with migratory waterfowl, and this interface may have played a crucial role in the emergence of this novel genotype of highly pathogenic H5N1 viruses [63,71,75]. Since 2011, in Bangladesh, clade 2.3.2.1a has been detected and reassorted in multiple host species, including chickens, ducks, geese, migratory wild waterfowl, quail, pigeons, and house crows [63,74,76]. Te fact that the current HPAI H5N1 clade 2.3.2.1a (new) is continuously evolving and widening its host reservoir, including spillover into mammalian hosts, including humans, in Bangladesh, is a worrying development. It underscored the need for continuing genomic surveillance and improved biosecurity measures to prevent new viral introductions and the spread between farms and other entities of the Bangladeshi poultry industry.
For our overview presented here, we are confdent that our data collation includes the most outbreak data. We cannot rule out some underreporting of outbreaks due to fear of compulsory culling and insufcient fnancial compensation. Despite these possibilities, it is unlikely that any  Figure 4).
missing data would lead to biases in the patterns observed.
Our study therewith provides an in-depth description of the spatiotemporal and evolutionary dynamics of H5N1 outbreaks in Bangladesh, which we hope will serve as a foundation for further studies aimed at improving our understanding of the drivers of HPAI H5N1 outbreaks, which can be used in the development of efective disease control strategies in Bangladesh.   Figure 8), highlighting the temporal pattern of the circulating clades. Each clade is denoted using a diferent color.

Conclusion
Our research explored the spatiotemporal dynamics of HPAI H5N1 outbreaks and the genetic evolution and clade diversity of the virus in Bangladesh. Our study identifed nine HPAI H5N1 epizootic waves between 2007 and 2020. We observed a signifcant decrease in the number of outbreak reports after the onset of HPAI H5N1 vaccination for commercial poultry in 2012. Nevertheless, post-2012 surveillance data reveals year-round presence of the H5N1 virus in farms and LBMs, indicating the virus is endemic in Bangladesh. Te application of spatiotemporal analyses showed the presence of HPAI H5N1 outbreak hotspots, with densely populated areas being most vulnerable to outbreaks, and providing crucial information for the development of efective disease control and prevention strategies. Since 2012, and despite a national HPAI H5N1 vaccination campaign, clade 2.3.2.1a dominates the avian infuenza landscape in Bangladesh, with strains from this clade being found in an increasing number of species, including mammals. In the face of this worrisome proliferation, we recommend continued monitoring of the evolutionary dynamics of HPAI H5N1, implementing stringent biosecurity measures in both LBMs and farms, and the monitoring of vaccine efcacy to put a stop to this rise of HPAI H5N1 and prevent incursions of novel strains of H5Nx viruses into Bangladesh.

Data Availability
Te data used in the analyses can be obtained from the corresponding author upon reasonable request.

Ethical Approval
No ethical consideration was required as the study uses open-source datasets.

Disclosure
Te content is solely the authors' responsibility and does not necessarily represent the ofcial views of the USA Government. An earlier version of the abstract to our manuscript was submitted to the 7 th World One Health Congress in Singapore, November 7-11, 2022, and published at https://www.epicentro.iss.it/globale/pdf/abstractbook_v2-1. pdf.

Conflicts of Interest
Te authors declare that they have no conficts of interest.