H9N2 avian influenza virus dispersal along Bangladeshi poultry trading networks

Abstract Avian influenza virus subtype H9N2 is endemic in Bangladesh’s poultry population. The subtype affects poultry production and poses a potential zoonotic risk. Insufficient understanding of how the poultry trading network shapes the dissemination of avian influenza viruses has hindered the design of targeted interventions to reduce their spread. Here, we use phylodynamic analyses of haemagglutinin sequences to investigate the spatial spread and dispersal patterns of H9N2 viruses in Bangladesh’s poultry population, focusing on its two largest cities (Dhaka and Chattogram) and their poultry production and distribution networks. Our analyses suggest that H9N2 subtype avian influenza virus lineage movement occurs relatively less frequently between Bangladesh’s two largest cities than within each city. H9N2 viruses detected in single markets are often more closely related to viruses from other markets in the same city than to each other, consistent with close epidemiological connectivity between markets. Our analyses also suggest that H9N2 viruses may spread more frequently between chickens of the three most commonly sold types (sunali—a cross-bred of Fayoumi hen and Rhode Island Red cock, deshi—local indigenous, and exotic broiler) in Dhaka than in Chattogram. Overall, this study improves our understanding of how Bangladesh’s poultry trading system impacts avian influenza virus spread and should contribute to the design of tailored surveillance that accommodates local heterogeneity in virus dispersal patterns.


Introduction
The endemic circulation of avian influenza viruses (AIVs) in Bangladesh's poultry populations poses a threat to animal and human health (Marinova-Petkova et al. 2016;Turner et al. 2017). Low pathogenic AIV (LPAIV) subtype H9N2 was first detected in Bangladesh in 2006 (Shanmuganatham et al. 2014;Marinova-Petkova et al. 2016;Rimi et al. 2019) and is known to cause reduced egg-laying and hatching (Kariithi et al. 2020;Ripa et al. 2021). Viruses from this subtype are now common within Bangladesh at live bird markets (LBMs) where most consumers purchase poultry and are often detected at relatively lower prevalence in poultry farms or backyard rearing sites (Gerloff et al. 2016;Turner et al. 2017;Kim et al. 2018;Parvin et al. 2020;Gupta et al. 2021; Moyen et al. 2021). Close and frequent contact between birds and humans in these locations increases the risk of zoonotic spillover (El-Shesheny et al. 2020;Parvin et al. 2020;Bi, Li, and Shi 2022). Reducing virus spread is therefore particularly important for minimising pandemic emergence risk and protecting avian health (El-Shesheny et al. 2020;Parvin et al. 2020;Bi, Li, and Shi 2022).
Live poultry trading is known to contribute to AIV spread (Wu and Perrings 2018;Yang et al. 2020). Previous studies have identified broad associations between growth in poultry trade volumes and the likelihood of establishment and persistence of several infectious diseases (Wu and Perrings 2018). However, we often lack a more detailed understanding of the complexity of bird production and distribution networks and how these networks may influence AIV maintenance and spread (Gerloff et al. 2016;Parvin et al. 2020;Moyen et al. 2021). In many countries where AIVs are endemic in poultry, knowledge of viral prevalence is too limited to adequately explore how trading practices might impact circulation from infection data alone (Chattopadhyay et al. 2018;Parvin et al. 2020;Gupta et al. 2021;Moyen et al. 2021;Ripa et al. 2021). Even when surveillance is routinely conducted, LPAIV H9N2 can easily be missed because it rarely causes severe disease (Parvin et al. 2020).
Phylodynamics can reveal in-depth information about virus dynamics from virus genome sequences, including how an outbreak changes size over time, how a virus lineage spreads spatially, and what factors may influence viral dispersal patterns (Lu, Leigh Brown, and Lycett 2017;Yang et al. 2019;Kwon et al. 2020). For example, phylodynamic studies have demonstrated that live poultry trade networks shape AIV movement over large spatial scales in movement in China (H5N1, H7N9, and H5N6) (Yang et al. 2020) and that restricting duck transport and culling can suppress highly pathogenic AIV (HPAIV) H5N1 movement between regions in France (Chakraborty et al. 2022). Phylodynamic approaches, however, have rarely been applied to studying AIV in a diseaseendemic poultry production system at high resolution (Yang et al. 2020). Within Bangladesh, specifically, the few published studies on the molecular epidemiology of AIVs have either relied on small numbers of samples (e.g. (Gerloff et al. 2016;Ripa et al. 2021)) or have not incorporated information on the precise sampling location or type of chicken (e.g. Marinova-Petkova et al. 2016;Parvin et al. 2019Parvin et al. , 2020. Hence, they do not permit full exploration of how the virus dispersal patterns vary between different components of the poultry system (Marinova-Petkova et al. 2016;Parvin et al. 2019Parvin et al. , 2020. Within Bangladesh's poultry production and distribution network, mobile poultry traders ('middlemen') collect and transport poultry from farms in rural and peri-urban areas to LBM vendors within urban areas. Here, poultry is either sold directly to endusers or may be traded further between market vendors before sale (Moyen et al. 2018(Moyen et al. , 2021Moyen 2019;Høg et al. 2021). The most commonly sold poultry types are broiler (exotic, industrial chicken breeds), sunali (chicken crossbreed of Rhode Island Red cocks and Fayoumi hens), deshi (indigenous chicken breeds), and ducks (Moyen et al. 2018;Gupta et al. 2021). Broilers and sunalis are raised on commercial farms, whereas deshis and ducks are raised in a traditional scavenging system ('backyard') (Moyen et al. 2018). Recent studies show that poultry trading practices vary substantially across the network (Moyen et al. 2018(Moyen et al. , 2021Moyen 2019). The two largest cities in Bangladesh, Dhaka and Chattogram, contain numerous LBMs where AIV infections have been consistently reported and receive poultry from largely nonoverlapping regions of Bangladesh (henceforth, 'production areas') (Moyen et al. 2018). Trade practices differ greatly between these cities, with more frequent inter-market trading of birds in Dhaka than in Chattogram (Kim et al. 2018;Moyen et al. 2018). Broiler chickens are typically sold in LBMs close to the farm, whereas ducks and deshi chickens are generally moved over longer distances by more intermediaries (Moyen et al. 2018;Moyen 2019). Typical distance travelled from farm to market, level of intramarket trade, and the number of middlemen varies between LBMs (Kim et al. 2018;Moyen et al. 2018).
This study investigates the spatial spread and dispersal patterns of H9N2 viruses among poultry marketed in Dhaka and Chattogram cities and their respective supply networks. Unlike previous studies in Bangladesh, we curated and analysed highresolution data on the location of sampling and type of chicken, enabling the first phylogenetic insights into how the complexities of poultry trading practices influence AIV dispersal there (Marinova-Petkova et al. 2016;Parvin et al. 2019Parvin et al. , 2020. We use phylodynamic analyses of the rapidly evolving haemagglutinin (HA) gene segment to examine the dispersal of H9N2 virus lineages between cities and stages of the production chain. We explore the hypothesis that H9N2 virus genetic diversity is randomly distributed across LBMs. Furthermore, we query how often H9N2 lineages spread between chicken types or locations and whether observed patterns of dispersal are driven by the overlap between chicken type production areas. Finally, to investigate if AIVs with different subtypes in Bangladesh may have similar dispersal patterns to those of H9N2 in the country, we perform similar molecular clock phylogenetic analyses on Bangladeshi HPAIV H5NX sequences.

Bangladeshi AIV sequences
Our study focuses on a dataset of AIV H9 HA segment sequences sampled between 2003 and 2019. The dataset includes eightytwo newly generated H9N2 HA sequences and H9 HA sequences for which only the HA segment was generated, but that are assumed henceforth to be from H9N2-infected birds (justification of this assumption is provided in Results and Discussion) (study sources in Table 1, accessions provided in Table S1). The dataset also includes 216 previously published sequences from both the Global Initiative on Sharing Avian Influenza Data (GISAID) (Shu and McCauley 2017) EpiFlu database (www.gisaid. org) (n = 211, Table S2) and from Ripa et al. (2021) (n = 5). Newly generated sequences were sampled through observational studies and routine surveillance in Dhaka and Chattogram from 2016 to 2018 (details in Table 1), one of which has been previously described (Kim et al. 2018).
Whilst we focus on HA here, sequences from other segments were generated for a subset of our H9N2 samples. NA gene sequences from Ripa et al. (2021) were combined with NA sequences from GISAID (Shu and McCauley 2017) EpiFlu database (H9N2; n = 200) corresponding to those viruses for which HA was available (Table S2) and used for a subset of analyses.
Associated information on species, location (either market-or farm-level), and sampling date was available for all newly generated sequences ( Table 1). The type of chicken (either broiler, deshi, or sunali) was also available for seventy newly reported H9N2 virus genome sequences generated through a cross-sectional study of Dhaka and Chattogram's LBMs during February-March 2016 (henceforth, the 'cross-sectional H9N2 dataset') (Kim et al. 2018) (Table 1). For sequences accessed from public databases, we extracted available corresponding location metadata from the Influenza Research Database (Squires et al. 2012) animal surveillance database (www.fludb.org) (Table S3). We also obtained unpublished LBM characteristics for 184 of the publicly available sequences in the H9N2 HA dataset from Professor Richard Webby (St. Jude Research Center of Excellence for Influenza Research and Surveillance) (Accession numbers in Table S4).
In addition to the H9N2 data described above, H5NX sequences were generated via next-generation sequencing through two of the three surveillance studies in Table 1. HA sequences (n = 29) and available corresponding NA gene segment sequences (n = 15) ( Table 1) were combined into datasets with previously published H5NX sequences from GISAID EpiFlu (HA, n = 175; NA, n = 172) (Table S2) and Ripa et al. (2021).

Sequence alignments and clade selection
We aligned sequences using MAFFT (Katoh et al. 2002) v7.453. We removed sequences from the alignment that either were duplicated, short (<70 per cent of the total sequence length), or indicative of containing sequencing or assembly errors.
We used an alignment of all H9NX sequences available on GISAID and our newly generated sequences to estimated preliminary phylogenies with Fasttree v2.1.11 (Price, Dehal, and Arkin 2010). All but seven HA H9 sequences from Bangladesh, H9N2; 1 HA sequence (accession numbers: Table S1). H5NX; 5 HA sequences (accession numbers: Table S1), 5 NA sequences (accession numbers: Table S1) including all sequences generated in this study, fell within a monophyletic clade containing H9N2 virus sequences sampled from poultry in India and Bangladesh (n = 333). We retained sequences from that monophyletic clade where exact sampling dates and at least district-level location data were known, resulting in 298 HA sequences in the H9N2 virus alignment (Bangladesh: n = 284; India: n = 14). We confirmed the presence of an appropriate temporal signal for both subtypes using TempEst (Rambaut et al. 2016) v1.5.3 ( Figures S1 and S2). Finally, for the respective available NA sequence dataset (H9N2, n = 205; H5N1, n = 191), we repeated all alignment and temporal signal checks as detailed above.

Molecular clock phylogenies
We investigated the introduction date of H9N2 to sampled locations using the Bayesian phylogenetic package BEAST (Drummond and Rambaut 2007;Suchard et al. 2018) v1.10.4. First, we compared four possible models: pairwise combinations of two molecular clock models (uncorrelated lognormal relaxed clock ) and strict clock (Marco and Ferreira 2008)) and two coalescent models (constant size (Kingman 1982;Griffiths and Tavaré 1994) and Bayesian skygrid (Gill et al. 2013)), all using an SRD06 substitution model (Shapiro, Rambaut, and Drummond 2006). We executed multiple Markov Chain Monte Carlo (MCMC) chains comprising 100 million steps and sampling every 10,000 steps. As identified using path-sampling (Lartillot and Philippe 2006) and stepping-stone-sampling analyses (Fan et al. 2011;Xie et al. 2011), the best model used the uncorrelated relaxed clock ) and a Bayesian skygrid coalescent prior (Gill et al. 2013). The posterior tree distributions for H9N2 and H5NX HA alignments were each obtained from two parallel MCMC chains with 250 million steps, sampling every 25,000 steps. We assessed the convergence of each run using Tracer (Rambaut et al. 2018) v1.7.1 (http://tree.bio.ed. ac.uk/software/tracer/) and confirmed the presence of appropriate parameter convergence by visually inspecting and then log combining multiple parallel runs. We summarised the information on maximum clade credibility (MCC) trees using TreeAnnotator (Drummond and Rambaut 2007) v1.10.4, with the first 10 per cent discarded as burn-in. We repeated all BEAST analyses for the NA sequence alignments using similar methods and parameters.
We used a generalised linear model (GLM) extension of the skygrid coalescent model (Gill et al. 2016;Dellicour et al. 2020) (henceforth, 'skygrid-GLM') to determine whether the effective population size (N e ) of the H9N2 HA lineage was associated with the tonnage of chicken or duck meat production in Bangladesh. Data were available yearly for each predictor from FAOSTAT (Food and Agriculture Organization of the United Nations, 2020) statistical database. To provide correspondence with the estimation of virus effective population size at two-month intervals, we used the 'zoo' package (Zeileis and Grothendieck 2005) in R (RStudio Team 2020) v4.1.2 to estimate values for each predictor every two months between 2000 and 2020 using linear interpolation. We then conducted the same analysis for the H5 HA dataset.

AIV transmission dynamics in live bird markets
We investigated whether individual H9N2 virus clades tend to be either randomly distributed across LBMs in a city or associated with specific markets. First, we identified genetically diverse 'clusters' (henceforth known as 'clades') in the HA H9N2 MCC tree using Clusterpicker (Ragonnet-Cronin et al. 2013) v1.2.3. These clades (n = 32) were identified according to a clade support threshold (0.7 posterior probability), which is consistent with values chosen in previous AIV phylogenetic clustering analyses (Gerloff et al. 2013;Lee et al. 2018) and a maximum sequence identity threshold (4.5 per cent genetic distance) that provided a consistent aggregation of closely related sequences (i.e. with most sequences into a specific well-supported clade). Next, to quantify lineage diversity within and between markets, we used the adjusted Rand index (Rand 1971) by means of the 'phyclust' package (Chen 2011) in R (RStudio Team 2020) v4.1.2 to determine the similarity of clustering by market and clade in the cross-sectional H9N2 dataset. To determine if there was a significant difference in the observed clustering pattern to that expected based on random viral movement, we compared the estimated median adjusted Rand index with a median adjusted Rand index calculated by permuting the market location 10,000 times while controlling for the city. We further generated binary adjacency networks in R (RStudio Team 2020) v4.1.2, in which the edges linked market nodes if they shared a genetically defined clade. To further test whether AIV clades are randomly distributed across markets, we determined the density of the observed network using the 'sna' package (Butts 2008) in R (RStudio Team 2020) v4.1.2 and subsequently compared it to the density value from permuting (n = 10,000) market location with respect to city. This density metric is the ratio of observed edges to the number of possible edges for the given network.

H9N2 dispersal between different chicken types
To determine whether H9N2 virus genetic diversity is structured according to Bangladesh's three main chicken types (broiler, sunali, and deshi) and two main cities (Dhaka and Chattogram), we used Bayesian Tip-association Significance Testing (BaTS) (Parker, Rambaut, and Pybus 2008) v1.0. We removed tips in each tree in the HA H9N2 posterior tree distribution when tip information on either chicken type and sampling city was unavailable, using the R (RStudio Team 2020) v4.1.2 library, 'ape' (Paradis and Schliep 2019) package. This generated a distribution of downsampled posterior trees each containing seventy tips, corresponding to those sequences in the cross-sectional dataset ( Table 1). The first BaTS analyses we performed assessed if sequences tended to cluster based on the sampling city (Dhaka and Chattogram). The second set of BaTS analyses determined if sequences tended to cluster by chicken type (broiler, sunali, and deshi); these analyses were conducted separately for each city because significant clustering of sequences by the city was observed in the first BaTS analysis. For each BaTS analysis, we computed median empirical values for each BaTS statistic (association index, AI; the parsimony score, PS) (Parker, Rambaut, and Pybus 2008) using 1,000 subsampled post-burn-in posterior trees from the respective phylogeny. We calculated P values by permuting the market location within each city 1,000 times to determine if empirical estimates significantly differed from the null expectation (based on random viral movement).
We also used a GLM extension of a phylogeographic discrete trait analysis (DTA) in BEAST v1.10.4 (Drummond and Rambaut 2007;Faria et al. 2013;Suchard et al. 2018) to determine whether covariates of the trading network predicted the H9N2 viral lineage movement. We used the six pairwise combinations of chicken types (sunali, deshi, and broiler) sold in Dhaka and Chattogram as the six discrete traits. The GLM predictors included (i) the number of sequences associated with each of the six categories as a proxy for sampling effort (both origin and destination), (ii) the weekly estimated sales at LBMs for all six categories (both origin and destination) as detailed in Moyen (2019) and Moyen et al. (2021), (iii) a binary chicken type similarity index of the six discrete traits (i.e. same chicken type = 1 and different chicken type = 0), (iv) an equivalent city similarity index, and (v) the extent of overlap between production areas from which birds are sourced (1, Pianka index) between each discrete trait. These production areas are the set of Upazilas (sub-district) where farms supplying markets within a city are located, which were previously computed from reconstructed transaction networks based on traders' interviews as described in Moyen et al. (2021). Descriptions for all GLM predictors included are summarised in greater detail in Table S5. We completed four separate MCMC analyses of 100 million steps sampling every 10,000 steps using the HA H9N2 empirical tree distribution downsampled to tips from the cross-sectional dataset (i.e. those tips with all appropriate metadata). The first model, termed model A, used all covariates except the chicken type and city similarity matrices. Models B and C were the same as model A but excluded the number of sequences or the weekly sales number, respectively. Model D used all five predictors. Model D was ran to check whether results obtained regarding the importance of production catchment area from models A-C were sensitive to the inclusion of chicken type and city similarity matrices, as these are somewhat correlated with production area matrices (Moyen et al. 2021). For each run, the presence of appropriate parameter convergence was confirmed visually in Tracer (Rambaut et al. 2018) v1.7.1 and multiple parallel runs for each MCMC analysis compared.

Dynamics of avian influenza virus in Bangladesh
All but seven H9N2 sequences from Bangladesh (of which six were sampled from ducks) fell within a H9NX monophyletic clade containing virus sequences sampled from poultry in India and Bangladesh, and this clade (n = 298 sequences with appropriate metadata) was therefore selected for further analysis. All sequences in this clade that had a respective NA sequence were identified as being H9N2 subtype, thus strongly indicating that any HA-only sequences (i.e. cross-sectional dataset) were also H9N2 (Price, Dehal, and Arkin 2010). The most recent common ancestor of all Bangladeshi HA sequences in this major H9N2 clade was estimated around October 2005 (December 2004-July 2006: 95 per cent highest posterior distribution (HPD)) (Fig. 1A), concurrent with the countries' first reported H9N2 cases in early 2006  (Fig. 1A).
Bayesian skygrid reconstructions from the HA segment ( Figures  S3A and B) indicate that the effective population size of H9N2 may have increased in Bangladesh since it was first detected, but neither chicken meat nor duck meat production was found to be a significant predictor of the effective population size (Figures S3C  and D).
Several clades in the HA H9N2 tree contain sequences from only Dhaka or from only Chattogram, suggesting that LPAIV transmission might occur preferentially between birds sold in each city (Fig. 1A). However, there are seven instances where sequence(s) from one of the cities fall within a clade where the basal sequence(s) to the clade were sampled in the other city (Fig. 1A). This is consistent with recurrent sharing of virus lineages between Dhaka and Chattogram or their introduction from the same source. The NA phylogeny contains very few sequences sampled in Chattogram (H9N2; Chattogram n = 4), making it difficult to assess whether H9N2 virus NA lineages are similarly structured between the two cities (Fig. 1B).
Only around 3 per cent of H9N2 HA sequences (8/298) were recorded as being sampled from farms. Farms have been sampled less intensively than LBMs in Bangladesh, in part because lower AIV prevalence and greater distance from urban research laboratories makes surveillance more challenging than at LBMs (Kim et al. 2018;Parvin et al. 2020;Moyen et al. 2021;Ripa et al. 2021). The farm-sampled HA sequences fall throughout the phylogeny (Fig. 1A), sometimes within clades containing sequences from both farms and markets sampled within similar time periods (Fig. S4). Farm-sampled genomes are too rare here to make robust conclusions regarding the direction of viral dispersal between farms and markets.
To begin to explore whether AIVs of different subtypes may show similar patterns to those of H9N2, we attempted similar analyses on HPAIV H5NX. The 95 per cent HPD interval of root dates for the H5 Bangladeshi sequences in both the HA and NA MCC trees overlap (Fig. S5). H5 HA sequences are proportionally less well sampled in Chattogram than in Dhaka (Chattogram, n = 6; Dhaka, n = 185) (Figures S5A and S6). However, the six H5 sequences sampled in the Chattogram fall into three different clades (Fig. S5A). Only 2 per cent of H5NX HA sequences were sampled from farms, but again these fell throughout the respective trees (Fig. S5A). Neither chicken nor duck meat production was a significant predictor of the effective population size of H5NX (Fig. S3). Whilst these results are extremely limited, there is therefore no clear indication from currently available data that H5 HPAIV exhibits strongly different patterns to H9N2.

Distribution of virus clades across LBMs
We investigated whether individual H9N2 HA clades tended to be randomly distributed across LBMs in a city or associated with specific markets. We quantified the diversity of viruses within and between markets in a city by employing the adjusted Rand index to measure the similarity of clustering by market and virus clade for the cross-sectional H9N2 dataset. We found no significant difference (P = 0.993) between the median adjusted Rand index calculated from the empirical and permuted data (Table S6). Therefore, differences in viral genetic diversity between and within markets were not significantly different. The binary adjacency network analysison the same cross-sectional H9N2 dataset found no significant difference (P = 0.881) between the median network density estimated from the permuted data and the median network density calculated from the empirical data (Table S7 and Fig. 2). Therefore, we found no significant difference in the sharing of clades between markets in a city.

Phylogenetic clustering by city and chicken types
To identify possible subnational differences in AIV dynamics between different major cities, we investigated whether HA H9N2 virus sequences sampled from Dhaka and Chattogram phylogenetically clustered by city and chicken type (deshi, sunali, and broiler) with BaTs (Parker, Rambaut, and Pybus 2008) v 1.0 (Fig. S7). The analyses demonstrate significant clustering of sequences by city based on AI and PS statistics (P < 0.001) ( Table S8). The AI statistic suggested sequences clustered significantly by chicken types in Chattogram (P = 0.036) but not in Dhaka (P = 0.278) (Table S9). In contrast to the respective AI statistic, the PS statistic did not significantly support clustering by chicken type in Chattogram (P = 0.071) (Table S9). This discrepancy could reflect the lower statistical power of the PS compared to the AI test (Parker, Rambaut, and Pybus 2008).
A GLM extension of discrete trait phylogeography was used to determine whether different features of the trading network predicted HA H9N2 virus dispersal between six pairwise combinations of chicken type (sunali, deshi, and broiler) and city (Dhaka and Chattogram). The results are summarised in Fig. 3 and Table S10. Our initial analyses (model A; Fig. 3A and Table S10A) involving all covariates bar the chicken type and city similarity matrices showed very strong support (Bayes factor (BF) = 27.8; following Stefan et al. 2019) for greater virus dispersal as the overlap in production areas from which chickens are sourced before being sold increased. This finding was consistent to subsequent analyses that excluded covariates for the number of sequences (model B; Fig. 3B and Table S10B) or the weekly sales number (model C; Fig.  3C and Table S10C) (BF > 70 and >50, respectively). However, when the city and chicken type similarity matrices were also included as predictors (model D; Fig. 3D and Table S10D), the association of greater production overlap with greater virus dispersal was no longer observed (production area matrices; BF < 1). Instead, the higher diffusion was strongly associated (city similarity index BF > 1,000; following Stefan et al. 2019) with the chicken presence in the same city. Thus, whether the bird was sampled in the same city is likely the primary predictor of virus dispersal of all covariates considered here (Moyen et al. 2018(Moyen et al. , 2021.

Discussion
Controlling AIV transmission in Bangladesh is complicated by our lack of understanding of whether virus spread disproportionately occurs within specific components of the poultry system (e.g., with farms, markets, or different chicken types as more important sources of infection). Here, we undertook a phylodynamic investigation to understand H9N2 viral spread in Bangladesh's poultry system, focusing on Dhaka and Chattogram markets and their respective production areas. For H9N2 analyses, we first identified Figure 3. Predictors of HA H9N2 dispersal among the six pairwise combinations of three main chicken types sold in Bangladesh (sunali, deshi, and broiler) and the two largest cities (Dhaka and Chattogram) in four separate DTA-GLM analyses (A-D). Inclusion probability is an estimate of the posterior expectation for the indicator variable associated with each predictor E( ). It suggests the likelihood that the predictor has a meaningful impact on viral diffusion. Bayes Factor (BF) support values for predictors (when >3 BF) are indicated by black text annotations. Coefficient ( | = 1) represents the contribution of each predictor on a log scale conditional when the predictor is included in the model, with the 95 per cent highest posterior density interval of the log GLM coefficients ( ) represented by horizontal lines from the mean. Grey boxes indicate that the predictor was not included in the specific DTA-GLM analyses. a monophyletic HA H9NX clade containing sequences sampled from poultry in Bangladesh and India, for which all available respective NA segment sequence were H9N2. Our analyses of this clade show that most H9N2 dispersal likely occurs between birds within the same major city, but that virus lineages are shared between Dhaka and Chattogram. H9N2 viruses from different lineages appear to be randomly distributed across a city's markets. We found regional differences in H9N2 virus spread suggesting more frequent viral transmission between chickens of different types in Dhaka compared to in Chattogram.
Our analyses suggest that H9N2 subtype lineage movement occurs less frequently between Bangladesh's two largest cities than within each city. Whilst present, virus genetic structuring by city appears less strong than previously hypothesised based on low levels of overlap in the geographic regions from which markets in Dhaka and Chattogram source birds (Moyen 2019). It is perhaps most likely that small overlaps in production areas allow for shared H9N2 lineage introductions to both cities identified here, given the frequent movement of birds from farm to markets and the relatively high prevalence of infection among traded chickens. However, more extensive and reliable data on mobile poultry traders' movements are required to robustly rule out possible alternative explanations (Moyen 2019;Høg et al. 2021;Moyen et al. 2021;Ripa et al. 2021), as all current estimates for chicken type production areas are based on reporting data that may be subject to memory recall errors (Høg et al. 2019;Moyen 2019;Moyen et al. 2021). Alternative or complementary explanations include rapid, direct viral movement between Dhaka and Chattogram (perhaps associated with direct poultry trading between individuals in each city) and slower indirect virus spread across unsampled intermediate locations (Høg et al. 2019;Moyen 2019;Moyen et al. 2021). Although the lack of wild bird sequences in our datasets prevents any robust conclusions of the contribution of wild birds to H9N2 mobility here, if wild birds played an important role in such viral movements, we may expect Bangladeshi sequences to be more frequently intermixed with samples from locations linked to the country via wild migration routes (e.g. China and Mongolia) (Tian et al. 2015;Lycett, Duchatel, and Digard 2019).
H9N2 viruses detected within one market are often interspersed phylogenetically with those from other markets, and genetically different viruses from this subtype appear randomly distributed across markets. This may indicate strong epidemiological connections between these markets, leading to frequent virus spread between them or from shared sources, as suggested by the highly connected trading network captured by Moyen et al. (2018). Our observations suggest that viral genetic diversity within a given LBM might not be hugely dissimilar to viral genetic diversity across all markets in a city, meaning surveillance of only a few markets within a city may be sufficient to capture viral diversity. This pattern could be generated in several ways, including frequent introduction of viral clades to LBMs (as suggested by Moyen et al. 2021) followed by limited intra-market persistence. However, we could not test this or other hypotheses regarding lineage introduction and persistence here because many markets in our dataset were represented by only one sequence.
Our analyses indicate that the H9N2 viruses generally cluster more frequently by chicken type in Chattogram than in Dhaka, suggesting that it may be necessary to consider subnational variation in production and trading processes when evaluating AIV dispersal. This finding could be consistent with several non-exclusive scenarios. In the first scenario, farms in different geographic regions might harbour geographically distinct viruses, which are then imported regularly to LBMs (Ripa et al. 2021). In that case, more significant overlap in the production areas from which different bird types are typically sourced could lead to greater mixing of viruses between those bird types. Weak evidence that H9N2 dispersal between different bird types in each city is associated with their degree of production area overlap somewhat supports this hypothesis. However, this finding is not robust to the inclusion of additional predictors. In the second scenario, our finding may result from differences in the trading practices between cities. Specifically, in Chattogram, middlemen typically sell poultry directly to stallholders (retailers or wholesaler), whereas in Dhaka, inter-market transactions are relatively more common (Moyen 2019;Moyen et al. 2021). As middlemen rarely supply more than one chicken type, there is less opportunity for direct viral transmission between chicken types during middlemen-facilitated transport to Chattogram's LBMs. In contrast, LBMs often sell multiple chicken types, so higher rates of market-to-market trading in Dhaka may increase opportunities for AIV dispersal between different bird types (Moyen 2019; Moyen et al. 2021).
Patterns of AIV dispersal identified here via H9N2 analyses may be different for non-H9N2 subtypes. For instance, HPAIV lineages may be relatively more detectable by poultry traders as host infections would often be more symptomatic (Parvin et al. 2018;Ripa et al. 2021). This could result in contrasting management practices for birds infected with LPAI and HPAI viruses and subsequently differing viral dispersal patterns (Rimi et al. 2019). Although our limited analyses of H5NX Bangladeshi sequences did suggest that similar patterns of intermixing of viral lineages between cities to those observed for H9N2 may exist, we could not repeat all analyses for H5 as for H9 due to a lack of available sequences and chicken type level metadata for the former subtype. Equally, while our H9N2 analyses focus mainly on sequences obtained from chickens, which are the hosts that H9N2 is more commonly found in and are the poultry type that accounts for more than 97 per cent of poultry moving through Dhaka and Chattogram's market stalls (Moyen, 2021), HPAIV lineages such as H5N1 are relatively more frequently detected in ducks (Kim et al. 2018;Kwon et al. 2020). Differences in trading patterns between chickens and ducks may, therefore, result in possible divergent virus transmission dynamics between such AIV lineages (Moyen et al. 2018(Moyen et al. , 2021Moyen 2019).
Our study has several limitations. First, virus genomic surveillance in Bangladesh is likely biased relative to infection distribution. Few sequences are available before 2013 in Dhaka or before 2016 in Chattogram, and change in genomic surveillance over time may therefore bias our results. There is a lack of AIV sequences from farms. Whilst this may be partly a true reflection of higher AIV prevalence at LBMs than at farms (Gupta et al. 2021;Ripa et al. 2021), this may also reflect surveillance bias towards LBMs, which are often in easier-to-reach locations than farms and where AIV is easier to detect. Finally, genomic sampling of H9N2 in Bangladesh is heavily biased towards Dhaka, and to a lesser extent, Chattogram, and therefore possible important sources and sinks of infection outside of these regions may be missed. Although H9 and H5 subtypes have been repeatedly detected in LBMs in both Dhaka and Chattogram (Sayeed et al. 2017;Biswas et al. 2018;Kim et al. 2018;Rimi et al. 2019;Hassan et al. 2020), the relatively lower number of Chattogram sequences may somewhat reflect both the higher number of live birds traded in Dhaka than in Chattogram (Moyen et al. 2021) as well as the greater quantity of mixed bird-type markets in Dhaka than Chattogram, which are generally associated with higher prevalence of infection (Kim et al. 2018). Biased sampling over space and time is less problematic for our analyses based on the cross-sectional dataset than those based on all publicly available data, as the cross-sectional dataset was generated through the same observational study conducted in both Chattogram and Dhaka. Our results show that subnational variation in production and trading processes may affect AIV dispersal between chicken types in Dhaka and Chattogram, and hence our results should be considered geographically specific and should not be extrapolated to other regions with production differences.
Second, for some analyses, we were limited to using only seventy sequences associated with known chicken types (Kim et al. 2018). This likely limited statistical power of our analyses of both market viral movement analyses and of city and chicken type clustering, relative to our analyses that did not use chicken type. Likewise, as these seventy sequences were generated in a single cross-sectional study (Kim et al. 2018) spanning only two months, we could not describe and explore seasonal variations in H9N2 infection patterns in such analyses. The dataset used for this analysis also contained only HA H9 sequences, meaning that we could not determine whether other segments or AIV subtypes may be similarly shaped by poultry trading practices in Bangladesh.
Improving insight into how viruses spread at a range of spatial scales could help guide improvements in AIV control. Our study shows the importance of recording accurate information on chicken type and highlights the need for greater surveillance on farms to understand viral epidemiology in Bangladesh. Our results suggest that nationally uniform interventions to reduce AIV prevalence may be unlikely to provide optimal effectiveness. Instead, actions should be tailored to the specific local structural characteristics of the poultry trading network and AIV dispersal patterns but could be made more efficient through targeted surveillance of a small number of key sites (Moyen et al. 2021). Any recommendations to improve AIV control in Bangladesh should adopt a multi-sector One Health approach to ensure proper consideration of health, social, and economic impacts (Chattopadhyay et al. 2018;Mackenzie and Jeggo 2019).

Data availability
The datasets and BEAST XML files used in this study can be found at https://github.com/lorcancarnegie/H9N2_Bangladesh.git; Gen-Bank Accession numbers for newly reported genetic sequence data are available under the accession numbers detailed in Table S1. Supplementary Data are available as separate files.

Supplementary data
Supplementary data is available at Virus Evolution online. the supplementary electronic material (Tables S2 and S3, respectively).
Funding L.C. was supported by a Biotechnology and Biological Sciences Research Council (BBSRC) London Interdisciplinary Biosciences Consortium (LIDo) studentship (BB/T008709/1). S.C.H. is supported by Sir Henry Wellcome Postdoctoral Fellowship from the Wellcome Trust (220414/Z/20/Z) (https://wellcome.org/). The work was funded from two sources: (1) the BALZAC research programme 'Behavioural Adaptations in Live Poultry Trading and Farming Systems and Zoonoses Control in Bangladesh' (BB/L018993/1), one of the eleven programmes funded under Zoonoses & Emerging Livestock Systems, ZELS, a joint research initiative between BBSRC, DSTL, DFID, ESRC, MRC, and NERC, and (2) the UKRI GCRF One Health Poultry Hub (BB/S011269/1), one of the twelve interdisciplinary research hubs funded under the UK government's Grand Challenge Research Fund Interdisciplinary Research Hub initiative.