Ecology and spread of the North American H5N1 epizootic

Dataset collection and processing

Information on case detections in North America

In this study, a detection is defined as a positive PCR test from a collected sample. In Canada, year-round surveillance in wild and domestic populations is coordinated by the Canadian Food Inspection Agency, Environment Canada, the Public Health Agency of Canada and the Canadian Wildlife Health Centre⁶⁴. In the USA, the United State Department of Agriculture Animal and Plant Health Inspection Service (APHIS) manages HPAI surveillance and testing in wild birds through investigation of reported morbidity and mortality events, hunter-collected game birds/waterfowl, sentinel species/live bird collection, and environmental sampling of water bodies and surfaces^43,65. USDA APHIS also surveilles domestic birds using several reporting methods: mandatory testing through the National Poultry Improvement Plan, coordination with state agricultural agencies, routine testing in high-risk areas and backyard flock surveillance⁶⁶.

Data on detections of HPAI in the USA used in analyses for this study were collected from USDA APHIS. Reports for mammals, wild birds and domestic poultry were all downloaded in November 2023 (download date: 25 November 2023)⁴⁰. During the time period analysed in this study (November 2021-September 2023), most HPAI detections in the USA were reported in wild birds (Supplementary Fig. 1a). Data on domestic bird detections are reported with information on poultry type (such as duck, chicken) and by whether the farm is classified as a commercial operation or backyard flock. Backyard flocks are categorized by the USDA as operations with fewer than 1,000 birds^47,67 and by the World Organization for Animal Health (WOAH) as any birds kept in captivity for reasons other than for commercial production⁶⁸. Among domestic birds, detections (1,177 total) came predominantly from commercial chickens (9.3%), commercial turkeys (28.5%), commercial breeding operations (species unspecified) (15.3%) and birds designated WOAH non-poultry, which refers to backyard birds (42.3%) (Supplementary Fig. 1b). Other domestic bird detections occurred in game bird raising operations (2.5%) and commercial ducks (2.0%). The North American epizootic has impacted a broad range of mammalian hosts, with detections (399) reported in red foxes (24.3%), mice (24.1%), skunks (12.2%) and domestic cats (13.2%). Other mammalian hosts (26.2%) represent a wide range of species including harbour seals, bobcats, fishers and bears (Supplementary Fig. 1c).

Genomic data processing and initial phylogenetics

We downloaded all available nucleotide sequencing data and associated metadata for the haemagglutinin protein of all HPAI clade 2.3.4.4b H5Nx viruses from the GISAID database on 25 November 2023 (ref. ⁶⁹). For each subset of the data described for further phylodynamic modelling, the following process was followed. We first aligned sequences using MAFFT v.7.5.20, sequence alignments were visually inspected using Geneious and sequences causing significant gaps were removed and nucleotides before the start codon and after the stop codon were removed^70,71. We deduplicated identical sequences collected on the same day (retaining identical sequences that occurred on different days). We identified and removed temporal outliers for all genomic datasets by performing initial phylogenetic reconstruction in a maximum-likelihood framework using IQtree v.1.6.12 and the program TimeTree v.0.11.2 was used to remove temporal outliers and to assess the clockliness of the dataset before Bayesian phylogenetic reconstruction^72,73. This resulted in a dataset of 1,824 sequences that were used in further analyses (Supplementary Fig. 17).

Biases in genomic data and N
_e inference

Sequencing data sampled in North America are heavily skewed toward sequences from the USA (USA, 1,590; Canada, 224; Central America, 8), and from the first 6 months of the outbreak, with 74% of all available sequences sampled from January to July 2022 (Supplementary Fig. 2). To evaluate whether sequencing data reflect case detections, we inferred the viral N_e—a measure of viral genetic diversity shown to be mathematically related to disease prevalence and the disease transmission rate²⁶. We inferred N_e using a nonparametric population model (Skygrid), which captures relative changes in genetic diversity and the variability of growth rate in the virus population over time, providing a proxy for epidemic dynamics as previously described. N_e is modestly correlated with detections (highest Spearman rank correlation: 0.65, P = 4.4 × 10⁻¹¹) (Fig. 1c and Supplementary Figs. 3 and 4), with peaks in N_e preceding peaks in detections by about 1 week (Supplementary Fig. 5), probably reflecting the lag between viral transmission and case detection. We interpret these results to suggest that, despite uneven sequence acquisition across time, the diversity of sampled sequences roughly reflect the amplitude of H5N1 cases. Given these results, we opted to use sequencing data for the entire sampling period for broad inferences on introductions and geographical spread, but supplement these analyses with controls for sampling differences between groups. For more-intensive reconstructions of transmission patterns between wild birds, commercial poultry and backyard birds, we focus on the initial 6-month period with the most densely sampled data, coupled with experiments to assess the impacts of sampling on results. Finally, although we retained data from Canada and Central America for all subsequent analyses, our results are probably most informative about transmission within the USA due to the heavy skewing of data towards the USA.

AVONET database

We downloaded the AVONET database for avian ecology data and merged it to available host metadata from GISAID for each sequence⁷⁴. We used the species if provided to match the species indicated in the AVONET database. If host metadata in GISAID was defined using common name for a bird, we determined the taxonomic species name and used that for further merging with the AVONET data (for example, ‘mallard’ was replaced with Anas platyrhynchos) for the given region to match the species to its respective ecological data. Domesticity status (whether a sequence was isolated from a wild host or a domestic host) was determined using available metadata downloaded from GISAID using the ‘Note’ and ‘Domestic_Status’ fields in sequence associated metadata. Moreover, if a given sequence strain name (in the field ‘Isolate_Name’) indicated domestic status (for example, A/domestic_duck/2022) these sequences were labelled as belonging to domestic hosts.

Phylodynamic analysis

The following Bayesian phylogenetic reconstructions and analyses were performed using BEAST (v.1.10.4)⁷⁵.

Empirical tree set estimation and coalescent analysis

We performed Bayesian phylogenetic reconstruction for each dataset before discrete trait diffusion modelling to estimate a posterior set of empirical trees. The following priors and settings were used for each subset of the sequencing data. We used the HKY nucleotide substitution model with gamma-distributed rate variation among sites and log-normal relaxed molecular clock model^76,77. The Bayesian SkyGrid coalescent was used with the number of grid points corresponding to the number of weeks between the earliest and latest collected sample (for example, for a dataset collected between 4 November 2021 and 11 August 2023, we would set 92 grid points)⁷⁸. We initially ran four independent MCMC chains with a chain-length of 100 million states logging every 10,000 states. We diagnosed the combined results of the independent runs diagnosed Tracer v1.7.2. to ensure an adequate effective sample size (ESS > 200) and reasonable estimates for parameters⁷⁵. If ESS was inadequate additional independent MCMC runs were run increasing chain length to 150 million states, sampling every 15,000 states were performed. We combined the tree files from each independent MCMC run removing 10–30% burn-in and resampling to get a tree file with between 9,000 and 10,000 posterior trees using Logcombiner v.1.10.4. A posterior sample of 500 trees was extracted and used as empirical tree sets in discrete trait diffusion modelling.

Discrete trait modelling framework

For each discrete trait dataset, we used an asymmetric continuous time Markov chain discrete trait diffusion model and implemented the Bayesian stochastic search variable selection (BSSVS) to determine the most parsimonious diffusion network⁷⁹. We inferred the history of changes from a given trait to another across branches of the phylogeny, providing a rate of transitions from A to B per year for each pair of trait states. When reporting these results, we refer to state A as the source population/state and B as the sink population/state. We implemented the BSSVS, which enables us to determine which rates have the highest posterior support by using a stochastic binary operator which turns on and off rates to determine their contribution to the diffusion network. In addition to the discrete trait diffusion rate, we used a Markov Jump analysis to observe the number of jumps between discrete states across the posterior set of trees and estimated the Markov rewards to determine the waiting time for a given discrete trait state in the phylogeny^49,50. The Markov reward proportion is calculated as the proportion of the phylogeny at a given time being a given discrete state. By looking at the proportion of a given state over time across the phylogeny, we can provide a proxy for how long transmission has occurred in each group between transition events. We calculate the transition rate as a realization of the CTMC process by dividing the number of Markov jumps by the tree height (branch length from the earliest tip to the root of the tree), and separately, by tree length (sum of all branch lengths). For each pairwise transition rate, we calculate the level of BF support that the given rate has. The BF represents the support of a given rate, and is calculated as the ratio of the posterior odds of the given rate being non-zero divided by the equivalent prior odds, which is set as a Poisson prior with a 50% prior probability on the minimal number of rates possible⁷⁹. We use the support definitions by Kass and Rafferty to interpret the BF support where BF > 3 indicates little support, a BF between 3 and 10 indicates substantial support, a BF between 10 and 100 indicates strong support, and a BF of greater than 100 indicates very strong support⁸⁰.

Empirical tree sets were used with the discrete traits defined for each sequence to perform discrete trait diffusion modelling. Each discrete trait model was implemented using three independent MCMC chains with a chain length of 10 million states, logging every 1,000 states. Runs were combined using LogCombiner v.1.10.4, subsampling a posterior sample of 10,000 trees/states. The BF support for transition rates were calculated using the program SPREAD3 (ref. ⁸¹). Maximum clade credibility trees were constructed using TreeAnnotator v.1.10.4.

Extraction of phylogenetic metrics

We calculated the transitions between states across branches of phylogenies estimated from ancestral state reconstructions using the Baltic python package⁸². To calculate the persistence of a given discrete trait, we used the program PACT v.0.9.5, which calculates the persistence of a trait by traversing the phylogenetic tree backwards and measuring the amount of time that a tip takes to leave its sampled state⁸³.

Dataset subsampling and definition of discrete traits

Geographical introductions analysis

We characterized the geographical introduction of HPAI into North America by randomly sampling 100 sequences from Europe and Asia for each year between 2021 and 2023 (total, 300 non-North American) and all available North American sequences across the study period. After removal of temporal outliers, this resulted in a dataset of n = 1,927 sequences annotated by continent of origin. The sequencing data available from North America broken down by country are as follows: USA (1,590), Canada (224), Honduras (2), Costa Rica (5) and Panama (1).

Migratory flyways analysis

To characterize geographical transmission within North America after introduction, we constructed a dataset of sequences subsampled based on migratory flyway. We used place-of-isolation data to match the US state or Canadian province that the sequence was collected from with the respective US Fish and Wildlife Service Migratory Bird Program Administrative Flyway³⁰. We subsampled 250 sequences for each flyway (Atlantic, Mississippi, Central and Pacific) to create a dataset of 1,000 sequences collected between November 2021 and August 2023. In addition to USFWS flyways, we defined four geographical regions going north to south based on latitude lines, with the following delineations for each group. We divided North America into four regions segregated by latitude, with the northernmost group above the 49° N parallel and the southernmost group below the 36° N parallel. We then sampled 916 sequences uniformly across these categories and inferred transitions between these regions.

Host order analysis

We classified sequences by host taxonomic order, inferring the host species using designations in the strain name and/or metadata to match species records in AVONET⁷⁴. To ensure that each discrete trait had an adequate number of samples for the discrete trait analysis of host orders, we combined orders in two instances based on taxonomic and behavioural similarity. The order Falconiformes (n = 14), representing falcons, was added to Accipitriformes (n = 363), which includes other raptors such as eagles, hawks and vultures. Pelecaniformes (n = 34), including pelicans, were grouped with Charadriiformes (n = 74, shorebirds and waders) due to their similar aquatic lifestyles and behaviours. Mammals were kept as a broad non-human classification as most samples were of the order carnivora (foxes, skunks, bobcats), apart from samples of dolphins (Artiodactyla) and Virginia opossum (Didelphimorphia). The following orders were omitted due to a low number of sequences: Rheaforimes (n = 2), Casuariiformes (n = 1), Apodiformes (n = 2), Suliformes (n = 7), Gaviiformes (n = 1), Gruiformes (n = 1) and Podicipediformes (n = 1).

Discrete trait approaches assume that the number of sequences in a dataset are representative of the underlying distribution of cases in an outbreak, resulting in faulty inference when this assumption is violated^27,39,53 and bias when groups are unevenly sampled^82,84,85. To account for differential sampling among these host order groups, we considered two distinct subsampling approaches. The first is a proportional sampling regime in which sequences are sampled proportional to the detections in each host group each month. This common sampling regime assumes that case detections in each group are the closest proxy for the case distribution in the outbreak, and attempts to align sampling with underlying model assumptions. However, this approach may not be appropriate if case detection is heavily biased between groups. For HPAI H5N1 in North America, detections in wild birds are primarily identified when humans report sick or dead birds to wildlife health authorities or wildlife rescues (Supplementary Fig. 1a), which may skew detections towards birds with dedicated rescue services or birds that reside in closer proximity to humans. For example, Anseriformes and raptors comprised 50.2% and 20.3% of all sequences, respectively, which could arise from high case intensity or a higher rate of case acquisition. A second, complementary subsampling approach is to sample sequences equally, meaning that sequences are sampled from each group in perfectly equal numbers. By forcing the number of sequences from each group to be equal, the transmission inference must be driven by the underlying sequence diversity in each group rather than by sampling differences. Given the high variation among detections within each host group, we opted to pursue both sampling regimes and focus on results that were concordant in both. We first performed an AI test to confirm that clustering was sufficient for discrete trait inference (Extended Data Table 1). Next, for each regime, we performed three independent subsamples, where the dataset was sampled either proportional to cases or equally. For the equal sampling regime, each dataset included 100 randomly sampled sequences per host group, except for Passeriformes, for which only 57 sequences were available. To account for variation across subsampled datasets, we combined the results for the three independent subsamples to summarize statistical support (Supplementary Fig. 12 and Supplementary Tables 4 and 5). Owing to similar tree topologies across replicates, we visualized the phylogeny of the dataset with the highest posterior support (equal order subsample 1) in the main text and make the results of all analyses available in supplement (Supplementary Fig. 8 and Supplementary Tables 4–11). Finally, to measure the effects of potential sampling bias on the inferred transition rates, we performed a modified tip-shuffle analysis. We generated 100 datasets in which the host tip assignments were randomly shuffled, re-inferred the host group at internal nodes and infer a mean root state probability for each host across the 100 shuffled datasets. We then compared the root state probability in the empirical data to that inferred in the shuffled data as a measure of the impact of sampling on the results as previously described (see the ‘Tip-shuffle analysis’ section for further details)³⁸.

For the equal sampling regime, we randomly subsampled 100 sequences for each host order between 4 November 2021 and 11 August 2023, resulting in a dataset of n = 655 sequences whereby all isolates for host orders with less than 100 samples, Passeriformes (n = 57) and Strigiformes (n = 99) (removing one temporal outlier), were used (Supplementary Fig. 18). We repeated this random subsampling three times, resulting in three separate datasets. For the proportional sampling regime, we performed three subsamples of sequences based on the proportion of detections in each host order group, which were collected between 4 November 2021 and 11 August 2023. Three random proportional samples were taken each with the following number of sequences for each group: Accipitriformes (133), Anseriformes (342), Passeriformes (12), non-human-mammal (16), Galliformes (83), Charadriiformes (40), Strigiformes (29) (total n = 655 sequences).

Migratory behaviour analysis

We defined discrete traits for use in discrete trait diffusion modelling based on the available sequence metadata and merged AVONET data. In addition to taxonomic order, we defined migratory behaviour. Birds were classified as sedentary (staying in each location and not showing any major migration behaviour), partially migratory (for example, small proportion of population migrates long distances, or population undergoes short-distance migration, nomadic movements, distinct altitudinal migration) or migratory (the majority of population undertakes long-distance migration). We subsampled sequences based on migratory behaviour including non-human-mammals and domestic birds to create a subsample of 500 sequences with equal sampling across behaviour groups.

Rationale for inclusion of turkeys as domestic birds

While commercial turkey operations represent 53.7% of all detections on commercial farms⁴⁴, the presence of wild turkeys throughout North America makes categorizing turkey sequences as domestic or wild status ambiguous. 98% of all turkey sequences are not associated with metadata on domestic/wild status, and were therefore excluded from the first analysis of domestic/wild bird diffusion. However, epidemiological data suggest that most deposited turkey sequences probably stem from domestic outbreaks. Among case detections during the study period, only 139 were reported in wild turkeys, representing 1.5% of all wild bird detections. By contrast, commercial turkey outbreaks comprised 28.5% of all domestic detections in the study period, suggesting that unlabelled turkey sequences are most likely to have come from domestic birds. While these data are not conclusive, we opted to perform an additional analysis to determine whether our exclusion of turkey sequences (that are most likely domestic) may have biased our results. In the analyses detailed below, turkeys are assumed to be domestic.

Domestic/wild titration analysis

To study the impact of sampling of wild birds on the estimation of rates between domestic and wild birds, we created five separate datasets with varying numbers of wild birds for sequences collected between 2021 and 2023. We randomly sampled 270 domestic sequences and 270 wild sequences as the initial 1:1 ratio dataset. We then made four more datasets increasing the number of wild sequences by a factor of 0.5 (adding 135 wild sequences), resulting in a final titration of 1:3 domestic to wild sequences (n = 1,080). We applied a two-state asymmetric CTMC discrete trait diffusion model in which sequences were labelled as domestic or wild. All priors and model parameters selected are the same as those described in the empirical tree set description above. To study the impact of the inclusion of turkeys in the transmission between domestic and wild populations, we annotated all unannotated sequences collected from turkeys as domestic (see the rationale in section above). We then created three datasets starting with 525 domestic and 525 wild bird sequences, adding 263 sequences to successive titrations resulting in 1:1, 1:1.5 and 1:2 (domestic:wild) sequencing datasets with a final titration size of 1,575 sequences. We again applied a two-state asymmetric CTMC discrete trait diffusion model in which sequences were labelled as domestic or wild, and all priors and model parameters selected are the same as those described in the empirical tree set description above. To determine whether the proportion of turkeys to other domestic birds would impact the results of the previously described titration analysis we built a dataset in which the domestic bird group had equal numbers of turkey and domestic (non-turkey) sequences. This dataset included 173 turkey, 173 domestic bird and 692 wild bird sequences, totalling 1,038 sequences. Given that turkeys comprised 53.7% of commercial poultry outbreaks in the study period, this sampling regime conforms to both equal and proportional sampling regimes. We applied an asymmetric CTMC discrete trait diffusion model using a BSSVS for a three-trait model with the following states: wild birds, domestic birds (not turkey) and turkey. We performed three independent runs of this analysis using the models and parameters described in the empirical tree analysis section above. All titration replicates were performed using an MCMC chain length of 100 million states sampling every 10,000 states.

Commercial, backyard, wild-bird titration analysis

Metadata and annotated sequences were made available describing sequences as being from backyard birds for sequences collected in early 2022 which distinguished them from commercial poultry (previously all sequences being determined domestic)¹⁵. We used these metadata to create a dataset with equally sampled backyard birds and commercial birds (n = 85 for each bird type) and then added all available wild birds (n = 722) in 25% increments creating four separate datasets for sequences collected between Jan 2022 and June 2022. This resulted in a final dataset of n = 942 sequences. We performed discrete trait diffusion modelling using an asymmetric CTMC diffusion model described in the previous section for sequences labelled as backyard bird, commercial bird and wild bird. Calculation of the lag time between the cumulative Markov Jumps for backyard birds and commercial birds was calculated as the average length of time between points where cumulative Markov jumps are equal between backyard birds and commercial birds. This was calculated for each tree in the posterior.

Assessment of sampling bias

BaTs analysis

To determine whether the discrete traits analysed correlated with shared ancestry in the phylogeny, we employed tip trait association tests implemented in the Bayesian Tip-Association Significance (BaTs) program (v.1.0)³¹. This program assesses the phylogenetic structure of discrete traits across viral lineages using three metrics: the AI, parsimony score (PS) and maximum monophyletic clade size (MC). The AI measures the imbalance of internal nodes of a phylogeny for a given set of traits. The PS calculates the number of state changes in the phylogeny. The MC measures the maximum number of tips belonging to a monophyletic clade for each discrete trait of interest. These metrics are calculated for the phylogeny as tips are randomly swapped to create a null distribution to compare against. Taken together, these metrics quantify the degree of clustering within the phylogeny, with lower AI and PS values indicating stronger phylogenetic structure, suggesting that closely related taxa tend to share the same trait, whereas higher values indicate weaker structure and more-frequent transitions between trait states. Statistical significance was assessed by comparing observed values against a null distribution generated through randomization, with P values reported for each test. All discrete trait groupings showed evidence for clustering by trait, supporting the use of trait modelling across the tree. The results of BaTS analyses for each discrete trait in this study are provided in Extended Data Table 1.

Tip-shuffle analysis

To assess the sensitivity of each of our discrete trait reconstructions to differences in sampling between groups, we implemented a modified version of a tip-swap analysis³⁸. As originally developed, a tip swap analysis attempts to assess the impact of trait sampling on discrete trait measurements. An operator is implemented within the MCMC chain that randomly picks pairs of tips and swaps their trait values, thus generating a posterior set of trees among which pairs of trait assignments have been randomly swapped. The probability of each state at the root is then computed, and compared to the inferred root state probabilities in the empirical data. As the root state probabilities in randomized datasets should primarily reflect the frequency of each trait in the analysis, empirical results that differ substantially from this null distribution are interpreted as evidence that the sequencing data are informing the analysis beyond what is expected based on trait frequencies alone. Thus, traits for which the root state probability differs considerably from the root state probability in the null data are frequently interpreted as being informed by the data, rather than sampling bias. While this approach has been shown to perform well on small phylogenies^1,86, the strategy of swapping single pairs of tips poses challenges for larger trees. In our flyways dataset, which includes around 1,000 tips, we found that, even with extremely high operator values (4,000), the traditional tip-swap analysis resulted in a posterior set of trees in which the majority of tips (93.3%) remained assigned to their true state at least 50% of the time, resulting in a null dataset that was only partially randomized. We believe that this is due to the high number of tips in our analysis, resulting in only an extremely small fraction of tips randomized at any given step in the MCMC chain. To overcome this limitation, we instead performed a randomized tip-shuffle analysis. Using the empirical set of trees inferred for each discrete trait analysis, we generated 100 null datasets in which we shuffled the trait assignments randomly across the tips. In this approach, we preserve the phylogenetic tree topology and the ratio of samples from each group, but shuffle their assignments at the tips. For each discrete trait analysis, we generated 100 distinct shuffled versions of the empirical trees, reran the analysis and summarized the resulting posterior distribution by inferring a maximum clade credibility tree. We then computed the root state probabilities for each trait for each MCC tree, and computed the mean root state probability across all 100 replicates. This computed mean is reported in Extended Data Table 2, in the column labelled ‘mean root state probability across 100 datasets with randomly shuffled tip states’. We then compared the root state probabilities in the empirical data (reported as ‘root state probability in empirical data’ in Extended Data Table 2) to the shuffled data as a measure of the impact of sampling on the results. As expected, the root state probabilities inferred in the shuffled datasets are proportional to the number of sequences included for each group. For the analyses using an equal sampling regime (migration, flyway, host orders equal and initial titration tests), this leads to approximately equal expected root state probabilities across groups. By contrast, the root state probabilities in the empirical data generally differ significantly from expectation, suggesting that the phylogenetic results are informed by the genetic data rather than from sampling alone. The results of tip shuffling analyses for each discrete trait in this study are provided in Supplementary Figs. 19–26 and Extended Data Table 2.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.