Transmission dynamics of low pathogenicity avian influenza (H2N2) viruses in live bird markets of the Northeast United States of America, 2013–2019

Abstract Live bird market (LBM) surveillance was conducted in the Northeast United States (US) to monitor for the presence of avian influenza viruses (AIV) in domestic poultry and market environments. A total of 384 H2N2 low pathogenicity AIV (LPAIV) isolated from active surveillance efforts in the LBM system of New York, Connecticut, Rhode Island, New Jersey, Pennsylvania, and Maryland during 2013–2019 were included in this analysis. Comparative phylogenetic analysis showed that a wild-bird-origin H2N2 virus may have been introduced into the LBMs in Pennsylvania and independently evolved since March 2012 followed by spread to LBMs in New York City during late 2012–early 2013. LBMs in New York state played a key role in the maintenance and dissemination of the virus to LBMs in the Northeast US including reverse spread to Pennsylvania LBMs. The frequent detections in the domestic ducks and market environment with viral transmissions between birds and environment possibly led to viral adaptation and circulation in domestic gallinaceous poultry in LBMs, suggesting significant roles of domestic ducks and contaminated LBM environment as reservoirs in maintenance and dissemination of H2N2 LPAIV.


Introduction
Live bird markets (LBMs) are implicated as an optimal environment for avian influenza virus (AIV) maintenance, replication, and spread among susceptible poultry species co-housed with daily introduced naive poultry leading to viral maintenance within the market. The Northeast (NE) United States (US) has the largest number of LBMs in the USA with a complex system of production flocks, dealers/haulers, and markets (Garber et al. 2007;Jagne, Bennett, and Collins 2021). The LBMs in the NE US, consisting of more than 140 markets across 8 states [Pennsylvania (PA), New York (NY), New Jersey (NJ), Connecticut (CT), Rhode Island (RI), Maryland (MD), Massachusetts (MA), and New Hampshire (NH)], serve as a major source of fresh poultry meat products for mainly ethnic immigrant populations in large cities (Senne, Pedersen, and Panigrahy 2005). Avian influenza surveillance in these markets is voluntary under the 'Prevention and Control of H5 and H7 Avian Influenza in the Live Bird Marketing System Uniform Standards for a State Federal-Industry Cooperative Program (USDA-APHIS 2020)'. As of 2021, the States of NY, PA, and NJ have the most active systems with NY State having the most markets presently at 87, followed by 36 in NJ, 11 in PA, and a small number (1-3) of LBMs in each of other states (Jagne, Bennett, and Collins 2021).
The frequent detections of AIV from LBMs raised public health concerns. The NE US LBMs were identified as the probable source of the virus responsible for the Pennsylvania H5N2 low pathogenic AIVs (LPAIVs) outbreak in the summer of 1983 which subsequently mutated to high pathogenicity (HP) AIV in October 1983 resulting in the HPAIV outbreak in 1983-4 in PA, VA, and MD (Senne, Pedersen, and Panigrahy 2005). During 1996-2006, the H7N2 LPAIVs spread from LBMs of the NE US causing five consecutive outbreaks in commercial poultry despite the efforts to eradicate the virus by market closures followed by extensive cleaning and disinfection (Akey 2003;Davison, Eckroade, and Ziegler 2003;Dunn et al. 2003). The voluntary LBM surveillance program was established in NE US LBMs to monitor AIVs after the H5N2 virus outbreak in PA during the mid-1980s (Senne, Pedersen, and Panigrahy 2005). The NE US LBMs have been inspected and tested at least once per quarter, including testing of AIV from five to eleven randomly selected birds of each bird type and environmental samples from each market (Jagne, Bennett, and Collins 2021).
The H2N2 subtype has been sporadically detected in NE US LBMs (Schafer et al. 1993;Joseph et al. 2015). Previous studies demonstrated that H2N2 viruses detected in the 1990s formed a monophyletic clade, indicating that the viruses were closely related and shared a common ancestry. This current LBM H2N2 lineage is neither related to the clade from 1990 nor are they related to the human pandemic 1957 H2N2 influenza A virus, which belong to a Eurasian lineage and demonstrated a low or intermediate risk to mammalian transmissions (Schafer et al. 1993;Makarova et al. 1999;Liu et al. 2004;Piaggio et al. 2012;Jones et al. 2014). In this study, 384 H2N2 LPAIV submitted to the United States Department of Agriculture's (USDA) Animal and Plant Health Inspection Services (APHIS) from active voluntary surveillance efforts at the NE US LBMs, including NY, CT, RI, NJ, PA, and MD during 2013-2019 were sequenced and analyzed to aid in tracing the origin and investigate the transmission dynamics by incorporating host species and sampling locations into statistical Bayesian phylogenetic analysis models.

Genome sequencing
A total of 431 H2N2 LPAIV collected from LBMs in the NE US and confirmed by the USDA APHIS during 2013-2019 were sequenced in this study. We excluded low coverage/poor quality data and obtained complete genome sequences of 381 viruses, including 288 from NY, 57 from NJ, 17 from CT, 14 from PA, 3 from RI, and 2 from MD. All eight segments of the isolates were amplified by multiplex RT-PCR and whole-genome sequenced by using the MiSeq system (Illumina, San Diego, CA, USA) as previously described (Lee 2020). Briefly, the Nextera XT DNA Sample Preparation Kit (Illumina) and the 500 cycle MiSeq Reagent Kit v2 (Illumina) were used according to the manufacturer's instructions to generate and sequence multiplexed paired-end sequencing libraries. Genome sequence assembly was carried out using IRMA v.0.6.7 (Shepard et al. 2016) and verified manually using reference-based assembly in SeqMan NGen v14 program (http://www.dnastar.com/tnextgen-seqman-ngen.aspx). Genome sequences were deposited in NCBI BioProject (BioProject accession code: PRJNA750416) and GenBank (accession nos. KY272856-KY272863, KY272971-KY 272986, MW727341-MW727380, and MN998501-MN998508).

Intravenous pathogenicity index in chickens
The in vivo lethality testing for 34 representative viruses was conducted according to OIE (http://www.oie.int/en/internationalstandard-setting/terrestrial-manual). The selected representative samples include the index detection from each state and an additional virus newly detected every 6 months, or virus possessing any change in the HA cleavage site sequences. Briefly, 0.2 ml of a 1/10 dilution of infectious allantoic fluid was inoculated intravenously into eight 4-to 8-week-old specific-pathogen-free chickens and chickens monitored for 10 days for mortality; death in 6 to 8 birds is considered highly pathogenic. The in vivo pathogenicity assays with live viruses were conducted at the National Veterinary Services Laboratories, Animal and Plant Health Inspection Service, US Department of Agriculture in Ames, Iowa, USA, and in accordance with approved institutional animal care and use protocols.

Phylogenetic analysis
Comparative phylogenetic analysis was conducted to trace the evolution and spread of the H2N2 viruses in NE US. All closely related HA sequences of wild-bird origin H2 viruses searched by the Basic Local Alignment Search Tool (BLAST) search on Global Initiative on Sharing All Influenza Data (GISAID) database were added to the dataset as reference sequences which were exclusively wild-bird-origin H2 virus sequences. Reference sequences for phylogenetic analyses of PB2, PB1, PA, NP, NA, M, and NS gene segments were selected by BLAST search on GISAID database and subsampled based on nucleotide sequence identity of 99-99.5 per cent using CD-HIT (Fu et al. 2012) in February 2020. We generated maximum-likelihood phylogenies of each gene segment (PB2, ntax = 419; PB1, ntax = 429; PA, ntax = 435; HA, ntax = 411; NP, ntax = 440; NA, ntax = 423; M, ntax = 439; NS, ntax = 408) using RAxML v8 and the GTR nucleotide substitution model, with among-site rate variation modeled by using a discrete gamma distribution (Stamatakis 2014). Bootstrap support values were generated using 1,000 rapid bootstrap replicates. We performed a root-to-tip regression analysis against sampling dates on the ML tree using Tem-pEst (http://tree.bio.ed.ac.uk/software/tempest/) to determine the temporal signal of the HA sequences. As this revealed significant clock-likeness, we performed two distinct Bayesian discrete trait phylogeographic analyses (DTA) using BEAST version 1.10.4 (Drummond and Rambaut 2007). (i) To investigate virus transmission between locations, we reconstructed the virus transmission history between states geographically using an ancestral state reconstruction approach with a Bayesian stochastic search variable selection to determine the most probable spatial transmission history. The closely related HA sequences of wild-bird origin H2 viruses from BLAST search added to the dataset were labeled as 'Wild Bird'. In this analysis, we defined geographic region as discrete nominal categories, including Wild Bird (13 sequences), 'New York' (NY) (35 sequences), 'Pennsylvania' (PA) (14 sequences), 'New Jersey' (NJ) (29 sequences), and 'Connecticut and Rhode Island' (CT_RI) (20 sequences). To account for potential sampling biases, sequences were subsampled to preserve diversity with respect to location (state) which resulted in 111 sequences. Particularly, the NY and NJ dataset was subsampled based on nucleotide identity at 99.5 per cent and 99.8 per cent using CD-HIT (Huang et al. 2010), respectively. (ii) The phylogenetic data were also used to infer the transmission dynamics between host species in NY category, including 'domestic Anseriformes in New York' (Dm_Ans_NY) (103 sequences), 'domestic Galliformes in New York' (Dm_Gal_NY) (105 sequences), 'environment in New York' (Env_NY) (80 sequences), thus providing evidence for which host type plays a major role in the spread of H2N2 LPAIV in LBMs. The environmental samples include viruses detected in swabs from floors, cages, drains, walls, scales, door handles, etc. at LBMs. For both analyses, Bayesian relaxed clock phylogeny of HA gene was reconstructed by using BEAST version 1.10.4 (Drummond and Rambaut 2007). The HKY + G nucleotide substitution model with an uncorrelated lognormal relaxed molecular clock was used along with a Gaussian Markov Random Field (GMRF) Bayesian skyride coalescent tree prior (Minin, Bloomquist, and Suchard 2008). Bayes factor (BF) was calculated using the SPREAD v. 1.0.7 (Bielejec et al. 2011) to identify the best supported viral transitions between discrete categories. We identified a transition as significant when posterior probability >0.5 and BF >3. The rate and number of viral transitions between discrete categories (Markov jump) and the time spent in a state between transitions (Markov reward) were estimated using stochastic mapping . The Markov Chain Monte Carlo (MCMC) was run in parallel for 3 chains, each with 40-100 million steps and samples across chains combined after 10 per cent burn-in. The parameters were analyzed with TRACER v1.5 (http://tree.bio.ed.ac.uk/software/tracer/) and all parameters had an effective sample size greater than 200. A maximum clade credibility (MCC) tree was generated using TreeAnnotator and visualized using FigTree 1.4.4 (http://tree.bio.ed.ac.uk/software/ figtree/). We analyzed posterior trees using the program PACT (http://www.trevorbedford.com/pact) to compute the number of transition events through time (0.25 years per section) and discrete state proportion through time.
To assess the presence of sampling bias in our analysis, we conducted the generalized linear model (GLM) analysis as an extension of DTA to identify the impact of viral source/sink sample sizes on our Bayesian discrete inference datasets using BEAUti v1.10.4 software (Drummond and Rambaut 2007). The GLM result indicates virus migration rates as a linear combination of coefficients and coefficient indicators, and predictors. The coefficient represents the effect size of predictors affecting the migration rate of the virus and the coefficient indicator describes if the predictor was included in the model. The sample size of each discrete state was used as a predictor to inform the viral transition rates of two separate DTA.

Results
Lethality testing of representative LBM isolates confirmed that the H2N2 isolates were of low pathogenicity (no mortalities). ML phylogenetic analysis for each of the eight gene segments was used to infer the reassortment events of H2N2 LPAIV (Fig.  S1). Each gene segment exhibited a well-supported monophyletic clade of H2N2 LPAIV, suggesting the H2N2 LPAIV evolved in the absence of reassortment during 2013-2019. The H2N2 viruses detected in the 1990s were not phylogenetically related to the currently circulating H2N2 viruses in the NE US (data not shown). Molecular dating analysis of HA gene demonstrated that the time to most recent common ancestor (TMRCA) of H2N2 LPAIV in NE US LBMs was estimated to be 1 March 2012 [95 per cent Bayesian credible interval (BCI): 26 September 2011-27 June 2012] ( Fig. 1 and S2). The mean substitution rate of HA gene was 6.617 × 10 −3 substitutions/site/year (95 per cent BCI: 4.67-8.80 × 10 −3 ), which is in range of global AI virus  substitution rates as described elsewhere (Chen and Holmes 2006). Source sink dynamics of the H2N2 LPAIV identified from the LBMs in NE US during 2013-2019 was analyzed by estimating the transition rate (TR), Markov jump count, and Markov rewards. The GLM analysis revealed the potential impact of sample size has negligible influence to cause a potential bias toward the origin or destination state in this analysis (source indicator: 0.391, source coefficient: 0.419, sink indicator: 0.416, and sink coefficient: 0.251) ( Table 1) (Fig. 3). The total Markov reward time was the highest in LBMs in NY (37.550,, followed by NJ (15.751, 95 per cent HPD: 11.902-20.508), PA (12.1960, 95 per cent HPD: 4.947-18.607), and CT_RI (5.190, 95 per cent HPD: 2.535-9.889) (Fig. 4A). Collectively, these results indicate that the LBMs in NY play a central role in the maintenance and spread of H2N2 in NE US.
As a secondary analysis to explore the extent of viral maintenance and spread between host species and environment within LBMs in NY, DTA was conducted with three discrete nominal categories: Dm_Ans_NY, Dm_Gal_NY, and Env_NY. The frequent detections of virus transitions from Env_NY to Dm_Ans_NY (TR: 0.941, BF: 2616.411, posterior probability: 1.000), Env_NY to Dm_Gal_NY (TR: 1.471, BF: 49733.920, posterior probability: 1.000), Dm_Ans_NY to Env_NY (TR: 0.980, BF: 4972.287, posterior probability: 1.000), Dm_Ans_NY to Dm_Gal_NY (TR: 0.873, BF: 4520.149, posterior probability: 1.000), and Dm_Gal_NY to Dm_Ans_NY (TR: 1.023, BF: 49733.920, posterior probability: 1.000), suggest the H2N2 LPAIV is frequently transmitted between poultry and environment, and maintained in Anseriformes poultry and environment of NY LBMs ( Fig. 5 and Table 3). The GLM analysis revealed the potential impact of sample size is negligible in this analysis (source indicator: 0.163, source coefficient: −0.088, sink indicator: 0.033, and sink coefficient: 0.005; Table 1 Fig. 4B). The Markov reward time data are also consistent with previous findings on the significant role of domestic ducks in LBMs to facilitate the introduction and transmission of AIVs (Kwon et al. 2020;Youk et al. 2020). Frequent viral transmissions from Dm_Ans_NY to Dm_Gal_NY and Env_NY support the idea that domestic ducks may have acted as a seeding population in the LBMs in NY that is concordant with the highest trunk proportion being identified in Dm_Ans_NY during 2015-2016 (Fig. 6). An underlying mechanism of dynamics of circulating H2N2 LPAIV could be due to the cycle of deposition of viral pathogens in contaminated environment and reintroduction of the virus into newly introduced naïve poultry in the market system (Aguero-Rosenfeld et al. 2005). Our data highlight the key role of domestic ducks and environmental virus depositions in the persistence of H2N2 LPAIV in LBMs in NY, followed by dissemination to LBMs in neighboring states.
We identified the NA stalk region deletions in 3 different clades designated as subgroup A which contains 3 NJ viruses collected in 2016, subgroup B which contains 30 viruses collected in NJ 2018-2019, and subgroup C which contains 4 NY viruses collected in 2019 (Fig. S1). The sequences in subgroup A

Discussion
The risk of AIV influx and spread within the LBM system is constantly present due to the large avian population from various sources residing in the markets (Cardona, Yee, and Carpenter 2009;Jagne, Bennett, and Collins 2021). In the NE US live bird marketing system, poultry are raised on special production farms or commercial farms that sell wholesale to dealers and haulers who in turn provide to retail LBMs (Jagne, Bennett, and Collins 2021). In particular, poultry producers from PA provide the majority of the birds going into the LBMs for all the NE States. Our phylogeography data coincided with this poultry production and supply flow supporting that the initially introduced H2N2 LPAIV to PA subsequently spread to NY which has the largest number of LBMs in the US. After the multiple virus introductions from PA to NY LBMs in late 2012 and early 2013, the H2N2 viruses appear to have been maintained in NY LBMs and frequently spread to neighboring states including reverse spread to PA LBM. We were not able to investigate the association between interstate poultry trade and its role on the spread of H2N2 viruses due to limited information for poultry trade and supply networks between LBMs. The H2N2 virus in the LBMs has independently evolved since March 2012 with the same gene constellation until mid-2019. During the preparation of this article, we detected reassortant viruses in NY LBM possessing PB2, PB1, PA, and NS genes from other North American lineage non-H2N2 wild-bird-origin viruses in late 2019 (data not shown). Ongoing surveillance and sequencing will aid in determining whether these represent transient or may become established in LBMs.
In the present study, two separate phylodynamic analyses were used to estimate the most probable transmission dynamics between the LBMs in the NE states of USA, and between host species and environment within the NY LBMs. Our analysis between locations provides strong evidence for a central role of NY LBMs in the maintenance, amplification, and spread of H2N2 viruses via LBM chains from late 2013 to early 2019. It should be noted that we used only HA gene for this analysis, but transmission dynamics of other genes were not investigated. Although all eight genes of viruses in our dataset were monophyletic in ML analysis, there is a possibility that minor reassortment events occurred between the H2N2 viruses which could possibly cause alterations in the tree topology of other genes and influence the posterior probability of NY state as a source location. The secondary DTA was conducted to further investigate the transmission history of H2N2 between host species and the environment within NY LBMs. This analysis highlights the key role of domestic ducks and environmental virus depositions in persistence of H2N2 LPAIV in LBMs in NY, followed by dissemination to LBMs in neighboring states. We used all viruses identified from NY LBMs in the analysis since the NY LBMs acted only as a source population after the initial virus introductions from PA LBM during 2012-2013. It is most likely that there was no virus influx from other states to NY LBM except only one case from PA to NY (A/Muscovy_duck/NY/17-031928-003/2017) during the first quarter of 2017 as shown in Fig. 3. The transmitted A/Muscovy_duck/NY/17-031928-003/2017 virus was negligible in the DTA model because it disappeared without further circulations in NY LBMs.
The dissemination of H2N2 LPAIV between LBMs was most likely promoted by the trading of birds via market chains. In the trading process, flocks are collected at the farms or LBMs and distributed via ground transportation to distribution centers or directly to LBMs (Senne, Pedersen, and Panigrahy 2005). Newly entering birds in the LBMs originate from various sources including commercial poultry farms and backyard flocks in neighboring states. New birds are introduced to the LBMs daily and co-housed or caged closely with precedent flocks providing optimal conditions for AIVs to replicate and adapt to the new hosts which may result in long-term viral maintenance. The infectious period of AIV in infected birds is well-established to range from 6 to 10 days (Stallknecht et al. 1990;Webster et al. 1992;Brown et al. 2006). During the period, infected avian hosts undergo continuous shedding of virions in their feces and respiratory secretions which may transmit the virus directly to other hosts by ingestion or inhalation or may lead to limited viral maintenance in the environment. The viral shedding of infected poultry can contaminate the LBM environment samples, including swabs from floors, cages, drains, walls, scales, and door handles, from which the virus may persist for days or weeks depending on the temperature and other environmental factors, and reinfect host avian species (Breban et al. 2009;Jagne, Bennett, and Collins 2021). Based on our analysis, the frequent H2N2 AIV detections in the poultry and waterfowl, as well as the market environments emphasize the potential role of H2N2 LPAIV contamination in the LBM environment as an ongoing source of infection for incoming naïve domestic poultry. The H2N2 LPAIV Markov reward time and transmission data between discrete host categories in the NY LBMs indicate extensive viral maintenance within domestic Anseriformes and in the LBM environments which could facilitate viral spread between the susceptible hosts. Consistent with this finding, housing susceptible poultry species with ducks in proximity was considered as a potential risk in NE US LBMs due to characteristics of ducks as a natural AIV reservoir (Webster et al. 1992). Previous surveillance studies in East Asian LBMs also suggested the domestic duck as a major contributor to the spread and maintenance of various AIV strains (Mase et al. 2005;Wee et al. 2006;Uchida et al. 2008;Kim et al. 2010Kim et al. , 2012Sakoda et al. 2012;Lee et al. 2014;Kanehira et al. 2015;Kwon et al. 2016Kwon et al. , 2017Kwon et al. , 2020. Long-term circulation of wild-bird-origin AIVs in LBMs allows viruses to adapt among domestic poultry affecting viral host specificity and/or virulence via genetic evolution. Due to the confined spaces in LBM settings, different avian species could be caged in the same market with close proximity and it was commonly practiced in LBMs before regulations, such as Uniform Standard Program, were enacted (Jagne, Bennett, and Collins 2021). Since waterfowl are the natural reservoir for AIVs, viral transmission risks are present from co-housing domestic Gallinaceous species and domestic waterfowl population from various sources. In previous studies, changes in the viral host range were related to the length of NA stalk region of AIVs (Castrucci and Kawaoka 1993;Guangxiang, Jeffrey, and Palese 1993). Truncation in NA stalk region of AIVs has been associated with viral host adaptation in gallinaceous species (Banks et al. 2001;Hossain, Hickman, and Perez 2008;Sorrell et al. 2010). The majority of the domestic host species affected by viruses with NA truncation have been domestic gallinaceous poultry. In the current study, the presence of multiple NA stalk region deletions supports adaptation, extensive transmission, and long-term circulation of H2N2 LPAIV, not only in domestic waterfowl but also in gallinaceous poultry in the LBMs.
The complete genome sequences produced in this study greatly expanded the existing data set of AIV genome sequences from LBMs in the USA and allowed us to determine the evolutionary history and molecular epidemiology of H2N2 LPAIVs detected in LBMs across spatial and temporal scales. Molecular analysis demonstrated that the current LBM lineage H2N2 LPAIVs that has circulated in the NE US LBMs since 2012 continues to evolve, has acquired markers of gallinaceous poultry adaptation as early as 2019, and has reassorted with wild bird lineage viruses on more than one occasion. In addition, transmission and maintenance of the virus in domestic ducks and the LBM environment were important epidemiological features of the H2N2 LPAIV. These results highlight the potential for domestic ducks and environmental contamination in the LBMs to serve as reservoirs to maintain and disseminate AIV, subsequently contributing to viral adaptation and circulation in domestic poultry. Coordinated response measures such as trace outs, and heightened biosecurity measures such as market sell-downs, closures for cleaning and disinfection have been recommended by a subcommittee of the Live Bird Marketing Working Group that was created to help address this issue. Continued genomic surveillance is one tool in aiding the eradication of this virus and is crucial for ongoing monitoring of AIV incursions into the NE US LBMs.

Supplementary data
Supplementary data is available at Virus Evolution online.

Funding
This study is supported by the US Department of Agriculture, Agricultural Research Service project no. 6040-32000-066-51S.
Dong-Hun Lee is partially supported by the Bio and Medical Technology Development Program of the National Research Foundation, funded by the government of South Korea (grant no. NRF-2018M3A9H4056535). The findings and conclusions in this publication are those of the authors and should not be construed to represent any official USDA or U.S. Government determination or policy.

Conflict of interest:
None declared.