Epidemiology and molecular characterization of avian influenza A viruses H5N1 and H3N8 subtypes in poultry farms and live bird markets in Bangladesh

Avian influenza virus (AIV) remains a global threat, with waterfowl serving as the primary reservoir from which viruses spread to other hosts. Highly pathogenic avian influenza (HPAI) H5 viruses continue to be a devastating threat to the poultry industry and an incipient threat to humans. A cross-sectional study was conducted in seven districts of Bangladesh to estimate the prevalence and subtypes (H3, H5, and H9) of AIV in poultry and identify underlying risk factors and phylogenetic analysis of AIVs subtypes H5N1 and H3N8. Cloacal and oropharyngeal swab samples were collected from 500 birds in live bird markets (LBMs) and poultry farms. Each bird was sampled by cloacal and oropharyngeal swabbing, and swabs were pooled for further analysis. Pooled samples were analyzed for the influenza A virus (IAV) matrix (M) gene, followed by H5 and H9 molecular subtyping using real-time reverse transcription-polymerase chain reaction (rRT-PCR). Non-H5 and Non-H9 influenza A virus positive samples were sequenced to identify possible subtypes. Selected H5 positive samples were subjected to hemagglutinin (HA) and neuraminidase (NA) gene sequencing. Multivariable logistic regression was used for risk factor analysis. We found that IAV M gene prevalence was 40.20% (95% CI 35.98–44.57), with 52.38%, 46.96%, and 31.11% detected in chicken, waterfowl, and turkey, respectively. Prevalence of H5, H3, and H9 reached 22%, 3.4%, and 6.9%, respectively. Waterfowl had a higher risk of having AIV (AOR: 4.75), and H5 (AOR: 5.71) compared to chicken; more virus was detected in the winter season than in the summer season (AOR: 4.93); dead birds had a higher risk of AIVs and H5 detection than healthy birds, and the odds of H5 detection increased in LBM. All six H5N1 viruses sequenced were clade 2.3.2.1a-R1 viruses circulating since 2015 in poultry and wild birds in Bangladesh. The 12 H3N8 viruses in our study formed two genetic groups that had more similarity to influenza viruses from wild birds in Mongolia and China than to previous H3N8 viruses from Bangladesh. The findings of this study may be used to modify guidelines on AIV control and prevention to account for the identified risk factors that impact their spread.

www.nature.com/scientificreports/ Influenza A virus is a negative-strand RNA virus belonging to Orthomyxoviridae family, and the virion carries surface proteins known as hemagglutinin (HA) and neuraminidase (NA) 1 . Based on their potential to cause disease in chickens, avian influenza viruses (AIVs) are grouped into two categories: highly pathogenic avian influenza virus (HPAIV) and low pathogenic avian influenza virus (LPAIV) [2][3][4] . AIVs infect chickens, turkeys, and other gallinaceous birds and inflict significant economic losses worldwide [5][6][7] . In terms of humans and poultry, Bangladesh is one of the world's most densely populated countries. The poultry industry in Bangladesh supports economic growth and poverty reduction in rural and urban areas by creating employment opportunities and food products 8 . LPAIVs and HPAIVs, including the highly pathogenic H5N1 viruses, have been found in waterfowl, pet birds, wild birds, and chickens in Bangladesh [9][10][11][12][13] . In Bangladesh, over 580 outbreaks of HPAI H5N1 have been reported in poultry and wild birds since 2007 [14][15][16] . Because the poultry industry accounts for 20% of the livestock sector in Bangladesh, the culling of an estimated 250 million diseased animals to date in response to these outbreaks causes food insecurity and negatively impacted the economic growth 17 . In addition, eight human cases of H5N1 have been reported in Bangladesh, with one fatality 18 . Human H5N1 infections have been reported in Vietnam, Thailand, Indonesia, Hong Kong, China, and Cambodia, all of which have a history of poultry exposure in LBMs and commercial and free-range farms, implying that both LBMs and farms can contribute to the spread of AIVs among poultry and from poultry to humans [19][20][21][22] . Furthermore, the co-circulation of LPAIV is a source of concern. In Bangladesh, the average number of reported poultry outbreaks caused by AIV subtype H5 per year has declined, dropping from 83 and 10 outbreaks in commercial and backyard poultry, respectively, in 2007-2012 to two and zero outbreaks in 2013-2019. As compensation policies were phased out, underreporting and AIV subtype H5 vaccination in commercial poultry might be among the causes for the reduction in the number of avian influenza outbreaks 23 . LBMs are the backbone of poultry trade in Bangladesh. Birds of different species and geographical origins are introduced into LBMs on daily basis and these birds may be housed together, allowing for local transmission of several virus subtypes and possible reassortment 20,22,24 . On the other hand, commercial and backyard poultry farming contribute to the development of the poultry industry of Bangladesh 25 . Waterfowl, either wild or raised in various settings, such as backyard, nomadic, and free-range, play an additional role in virus transmission among commercial and wild bird populations 26,27 . Most studies aimed at determining the level of viral circulation in poultry are conducted in LBM, with only a few conducted-on poultry farms 24,28-30 . To address this gap and inform AIV control in poultry in both LBMs and farms, a thorough knowledge of infection patterns and risk management is essential. The study presented here quantifies the extent of H5, H9, and H3 virus circulation, factors influencing the spread of AIVs, and phylogenetic analysis of viruses in duck and turkey farms and LBMs in seven representative districts of Bangladesh. Overall prevalence of IAV subtypes (H3, H5, and H9). Figure 2 presents the prevalence of influenza A virus subtypes. IAV subtype H5 had the highest prevalence (22%, 95% CI 18.58-25.85). The prevalence of IAV subtype H3 and IAV subtype H9 was 3.4% and 6.9%, respectively.

Prevalence of influenza A virus M-gene.
Prevalence of IAV subtypes (H3, H5, H9) in host taxa. IAV subtype H3 was detected only in waterfowl (6.88%, 95% CI 4.31-10.80). On the other hand, IAV subtype H5 was detected in all the species. The prevalence of subtype H5 was highest in chickens (38.09%, 95% CI 20.29-59.81), while 25.51% (95% CI 20.45-31.32) of the waterfowl tested positive for IAV subtype H5. Further, 25% of pigeons and 16% of turkeys were found to have IAV subtype H5. However, IAV subtype H9 was not detected in waterfowl. Only 12.71% of the turkey and 4.76% of the chicken samples were positive for IAV subtype H9, as shown in Fig. 3.
Choropleth map of the prevalence (IAV subtypes H5, H3, and H9) by the district. The prevalence of IAV M-gene was highest in Thakurgaon. All the samples collected from Thakurgaon were positive for IAV M and subtype H9 (Fig. 5). The prevalence of IAV M-gene was 55% and 49% in Dhaka and Sirajganj, respectively. IAV subtype H5 was detected in 39% and 38% of the samples from Dhaka and Sirajganj, respectively. While the prevalence of IAV subtype H9 was 7% in Dhaka, 0.99% and 1.41% of the samples were IAV subtype H3 positive in Dhaka and Sirajganj, respectively. The prevalence of the IAV M-gene was 38% and was 17.57% for IAV subtype H3 in Kushtia. However, a very low prevalence was recorded for IAV subtypes H5 (4%) and H9 (1%). The prevalence of IAV M-gene was very low in Meherpur (4%) and Rajshahi (5%). In addition, no samples from Rajshahi and Meherpur were detected as being H5 positive. All samples from Naogaon tested negative for the M-gene.  health status significantly affect M-gene prevalence. However, the season is not a significant factor for subtype H5, but host taxa, health status, and interface are significant factors. The results of the multivariable logistic regression analysis are shown in Table 2. In winter, the chances of being positive for M-gene were 4.93 times higher than in summer, and waterfowl had 4.75 times the odds of being positive for M-gene and 5.71 times the odds of being positive for subtype H5 than chicken. In terms of health, dead birds were substantially more likely than healthy birds to be both M-gene positive and subtype H5 positive. Furthermore, a bird from LBM had a 3.28 times higher chance of being subtype H5 positive than a bird from a farm (  Fig. 6 and Supplementary Fig. 1). Four substitutions known to play a role in increased alpha-2,6 sialic acid receptor binding have been identified  www.nature.com/scientificreports/ in the HA sequences of the H5N1 strains isolated in the current study: D94N, S155N, T156A and K189R (H5 numbering). No classical neuraminidase inhibitor resistance markers were detected 31 .
H3N8. Phylogenetic analyses of HA and NA segments were also performed to gain better understanding of the evolutionary relationships between the characterized 12 H3N8 viruses (eight from Kushtia, two from Dhaka, one from Rajshahi, and one from Sirajganj). The HA segments of the twelve H3 viruses were grouped into the Eurasian lineage (Fig. 7). Ten H3 viruses were identical and clustered in one group; the other two viruses were identical and clustered in another group. The HA sequences of the twelve H3 viruses were genetically more similar to sequences of viruses from Mongolia, China, and Japan than the previously reported H3 viruses from Bangladesh. The twelve N8 neuraminidase genes grouped into the Eurasian lineage. The N8 sequences were closely related to those of viruses from Mongolia and China ( Supplementary Fig. 2).

Discussion
AIVs including highly pathogenic strains, can cause zoonotic disease and are often associated with severe economic, animal, and public health consequences 32 . Influenza A viruses (IAVs) subtype H5 have raised significant concern worldwide due to their high pathogenicity and zoonotic potential. More concerningly, three human  www.nature.com/scientificreports/ infections with avian influenza A (H3N8) virus were reported from China as of April 2023 33 . As a result, determination of the variables that impact influenza virus prevalence is essential. This study investigated the risk factors associated with influenza A, specifically H5 subtype prevalence. The present study detected influenza A viruses in waterfowl, chicken, turkey, and pigeon, similar to previously reported data [34][35][36] . We screened avian samples and found evidence of H3, H5, and H9 viruses in the tested samples, albeit with specific patterns. No H7 subtype viruses were detected in the tested samples. IAVs subtype H9 were detected only in turkey and chicken samples, and were not detected in waterfowl samples. A higher prevalence of H9 in chickens than waterfowl in LBM was previously reported 37,38 similar to the present study. The prevalence of IAV subtype H3 in the present study was lower than that reported in other South Asian countries 39 . We detected IAV H3 subtype only in waterfowl, which is consistent with the finding that H3 viruses are among the most commonly identified subtypes in ducks (proportion up to 91.76%) 38,39 . In addition, previous IAVs H3 subtype from Bangladesh were detected primarily in ducks 41,42 . IAV subtype H3 was detected in Kushtia, Sirajganj, Rajshahi, and Dhaka in the present study and more often in farms than LBMs, suggesting wild bird involvement through interactions with free-range farmed ducks.
In the present study, IAV subtype H5 had the highest prevalence. IAV subtype H5 has been spreading in Bangladesh since 2007, with more than 500 outbreaks reported in chickens. Analysis showed that all IAV subtype H5 in the tested samples belonged to the 2.3.2.1a clade. It was also found that the prevalence of influenza A M-gene and influenza A subtypeH5 were greater in chickens than in waterfowls, but the logistic regression model showed an increased risk of the IAV M-gene and IAV subtype H5 in waterfowls, similar to a previous report 43 .
Results of the current study showed that the risk of AIVs was higher in the winter than in the summer, consistent with many outbreaks in Bangladesh and other countries occurring during the same season 35,[44][45][46] . This may be explained due to low temperatures and low humidity, which favor AIVs persistence 47 . The odds of IAV H5 subtype detection were higher in LBM than in farms in the present study. The LBM biosecurity in Bangladesh is widely seen as falling short of minimum acceptable standards, and their capacity to act as viral evolution drivers facilitating the genesis of new emerging strains is significant 48 .
Molecular markers of influenza HA and NA are associated with increased virulence, adaptation to mammals, or resistance to antiviral agents. Unfortunately, no information on avian H3N8 virus molecular markers is available in literature, except for mutation W222L in HA was shown to allow for equine to dog adaptation for influenza A H3N8 49 . However, many molecular markers HA and NA are well characterized for H5N1. In the present  www.nature.com/scientificreports/ study, four substitutions D94N, S155N, T156A and K189R (H5 numbering) which were previously reported to play a role in the increased α -2,6 SA receptor binding have been identified in the HA of influenza A H5N1 in Bangladesh and no classical neuraminidase inhibitor resistance markers were detected in the NA sequences. Among the limitations of the present study is that only six H5N1 and twelve H3N8 viruses were selected for sequencing and only HA and NA were sequenced. This limited our ability to conduct accurate spatial and temporal analyses of the data. Detailed sequence analysis of the internal gene segments of the tested H5 viruses may add significant data regarding whether these viruses are similar to previously circulating AIVs in the sub-continent or whether these viruses are further evolving by reassortment with other subtypes. However, sequencing of internal gene segments of the tested influenza viruses in the present study was not possible due to lack of resources.
The findings of the present study could also be extended to include biosecurity practices such as disposal of dead bird, handling of sick bird, and other factors influencing AIV transmissibility and prevalence.
In conclusion, the findings of the present study highlight that multiple subtypes of AIVs (H5, H3, and H9) are circulating in both LBM and farms in Bangladesh. We identified the season, host type, and health status of birds as risk factors influencing AIV prevalence in the selected study areas. Farmers and workers in poultry farms are highly encouraged to get trained regularly to prevent possible AIV zoonotic transmission. The present study shows that avian influenza viruses are circulating in poultry populations in LBMs in Bangladesh, emphasizing the need for more extensive genomic surveillance studies to determine the avian influenza gene pool, monitor the possibility of emergence of new viruses of zoonotic potential, and to track their transmission. Routine and programed surveillance will help the early detection and rapid responses to prevent possible AI outbreaks. Many developing countries have similar poultry rearing practices in LBMs and farms as in Bangladesh. As a result, such recommendations based on our data and findings have global significance and implications.

Materials and methods
Sample collection. We conducted a cross-sectional study at the LBM and farm interfaces. In the instance of LBM, only LBMs in Dhaka, the capital and largest city in Bangladesh, were sampled due to the large numbers and concentration of LBMs and the presence of an AIV gene pool. On the other hand, for farms, the majority of live birds in Dhaka originate from the country's northeast part. In light of this fact, we collected duck and turkey samples from farms in seven major districts in the northeastern part of the country. Cloacal and oropharyngeal swab samples were collected from chickens (Gallus gallus domesticus), domestic ducks (Anas platyrhynchos domesticus) pigeons (Columba livia domestica) and turkey (Meleagris gallopavo) birds between September 2018 and November 2019 from seven districts, as shown in Fig. 1, using sterile swabs which were placed in a 1.8 mL sterile cryo-vial containing 1 mL of viral transport media (VTM) as previously described 49 .

RNA extraction and AIV rRT-PCR.
Pooled cloacal and oropharyngeal samples (n = 500) were shipped to the Division of Virology, Department of Infectious Diseases, St. Jude Children's Research Hospital (Memphis, TN, USA) for virus isolation and characterization. We extracted the viral RNA using the RNeasy Mini Kit (Qiagen, USA), and cDNA was synthesized using the SuperScript™ III Reverse Transcriptase (Invitrogen, USA). We tested all swab samples using real-time reverse transcription PCR (rRT-PCR) using universal M gene and H5, H3, H7 and H9 HA-specific primers, as previously described 51,52 . We considered a positive result for any gene tested if the Ct (cycle threshold) was less than 35 and with a characteristic amplification curve.

Virus isolation. Influenza A virus M-gene positive samples were cultured for virus isolation by inoculation
in embryonated chicken eggs (ECEs) as previously described 32 . A volume of 100 μl from each swab sample was injected into the allantoic fluid (AF) of three 10-day-old ECEs per sample and incubated at 35°C. After 72 hr, the AF was collected and tested for hemagglutination using 0.5% turkey erythrocytes. All the first passage in ECE (E1) AF samples were tested using the Flu Detect® (Zoetis Inc.) and RNA was extracted from AF samples and was retested for AIV by M gene rRT-PCR. The AIV-positive RNA were submitted for sequencing. M-gene positive, H5/H9 negative samples were subtyped using molecular methods using DNA sequencing. A multisegment RT-PCR-was performed using gene-specific primers according to the previously published protocol to amplify the whole influenza viral genome 53 . PCR products were then gel-extracted and purified using GE Healthcare illustra™ GFX PCR DNA and Gel Band Purification Kit (Sigma Aldrich, MO, USA). Library preparation of samples was performed using Illumina's Nextera XT DNA Sample Preparation kit (Illumina, CA, USA) according to the manufacturer's protocol. Amplicons were sequenced on the Illumina's MiSeq platform using the paired-end approach. Sequencing reads were demultiplexed, quality-trimmed and filtered prior to consensus sequence generation using the Pallas pipeline, developed by the Hartwell Genomics Centre at St. Jude Children's Research Hospital. Final analysis and the generation of consensus sequences were completed using CLC Genomics Workbench (v.11.0.1).

Phylogenetic analysis and molecular characteristics.
Multiple sequence alignment with the MUS-CLE algorithm and phylogenetic and evolutionary analyses were conducted using IQ-Tree 54 and MEGA11 55 . Briefly, phylogenetic trees of full-length HA and NA gene segments of six H5N1 viruses and 12 H3N8 viruses were generated using the maximum likelihood method with Tamura-Nei model and 1000 bootstraps replicate. All HA and NA sequences of gene segments of AIVs subtype H5N1 from Bangladesh available in the https:// www. ncbi. nlm. nih. gov/ genba nk/ and in the EpiFlu database of the Global Initiative on Sharing All Influenza Data https:// gisaid. org/ were retrieved and used as references, together with the first 50 blast hits for each gene www.nature.com/scientificreports/ segment and subtype, and for H3 and N8 trees together with reference sequences from the known genotypes as previously described 56 . In addition to H5N1 viruses from Bangladesh, the phylogenetic analyses also included H5N1 viruses from other countries with the closest genetic similarity, as indicated by the BLAST search. Bivariate analysis and statistical modeling. Bivariate analysis was conducted for different factors related to IAV M-gene and H5 in the study area. Seasons were classified into temporal patterns (winter and summer). Host taxa were categorized into pigeon, turkey, chicken, and waterfowl. Health status was classified based on healthy and dead birds. The interface was divided into LBM and Farm. "Chi-Square Test" was performed using R to identify the association of factors with the prevalence of the IAV M-gene and H5. Variables that had a p-value < 0.05 were considered for further statistical modeling. For statistical modeling, we utilized significant variables in the bivariate analysis. As our outcome variables on IAV M-gene, and IAV subtype H5 were binary, we considered multivariable logistic regression for each outcome.
Ethical approval. All procedures and methods were carried out in accordance with relevant guidelines and regulations and were reported in accordance with ARRIVE guidelines. All procedures were approved by the Ethics Committee of the Chattogram Veterinary and Animal Sciences University, Bangladesh, under reference number CVASU/Dir (R&E) EC/2019/126/(1).

Data availability
Sequence data generated in the present study were deposited in GenBank,