These results indicate that the HPAI H7N3 virus that caused the outbreaks in Mexico is not related to any of the previous H7N3 HPAI outbreaks in North America, nor related to other AI outbreaks (HP H5N2 outbreaks in Mexico in 1994–1995, LP H7 outbreaks in Canada in 2009) in domestic birds in recent years [14] , [15] . In addition, no clear pattern of association among the segments of the Mexico H7N3 strains was observed, indicating multiple segment exchange events occurred among North American influenza strains to give rise to it.

In contrast, the three polymerase encoding gene segments PB2, PB1 and PA of the Mexico outbreak strain belong to lineages composed mainly of AIV found in wild waterfowl in California from the beginning of 2012 ( Figure S3 , S4 , S5 ), with segments PB1 and PA having the same most closely related strain: A/American green-winged teal/California/123/2012 (H1N1). The other internal segments of the Mexico H7N3 strains have a different origin. From the Bayesian phylogenetic tree of the NP segment, the closest AIV strain to H7N3 Mexico is an H11N9 strain isolated from Mississippi in 2012 ( Table 1 ). The NP segment of the H7N3 Mexico outbreak strain is also in the same lineage as a small number of H7N7 AIV strains carried by northern teal in Illinois/Missouri in the fall of 2010; these strains belong to the same lineage in the HA segment as well (see above and Figures 2 and S6 ). The NS segment was derived from an H10N7 AIV which was also circulating in the same region at that time ( Figure S8 ). In the M segment, however, Mexico H7N3 strains are distinct from all currently available AIV in North America, suggesting a surveillance gap ( Figure S7 ).

Branches are coloured according to the 4 discrete traits (host order, host species, flyway and location) on internal nodes. Mexican outbreak strains are highlighted with pink. A: Host order. Five host orders are labelled on HA tree: wild birds of the order Anseriformes (ans-wild); wild birds of the order Charadriiformes (cha-wild); wild birds of the order Passeriformes (pas-wild); domestic birds of the order Galliformes and Mexico H7N3 outbreak in the order Galliformes (gal-domestic-Mexico). B: Host species. Wild Anseriformes are classified into the five main species and a group comprising the other rarer species of Anseriformes in this study: mallard (Anas platyrhynchos), northern pintail (Anasacuta), northern shoveller, blue-winged teal, green-winged teal and other Anseriformes (other ans); The order Galliformes are shown as “outbreak” (the H7N3 Mexico outbreak) and “other_gal”; The other orders are shown as: Charadriiformes (cha) and Galliformes (gal), Gruiformes (gru) and Passeriformes (pas). C: Flyway. Four specific North American flyways are labelled on the HA tree: the Atlantic, Mississippi, Central, and Pacific. D: State. 22 states and provinces of the viral sample locations are labelled on the HA tree. The original MCC tree files with all taxa names are deposited in Dryad (doi:10.5061/dryad.j5bf8), and trees for the other 7 segments without taxa names can be found in Figure S2 , S3 , S4 , S5 , S6 , S7 , S8 .

Co-circulation of multiple H7 clades was observed in HA across North America. Interestingly, the HPAI H7N3 Mexico strains are not related in HA to the HP H7N3 outbreak in British Columbia in 2004 and 2007, but instead are closely related to a subgroup of H7 AIV (H7N3, H7N8 and H7N9) from wild waterfowl isolated from Nebraska, Illinois, Missouri and Mississippi in 2010 and 2011. The mean estimate of the date of the common ancestor is February 2010 ( Figure 2 ). On the other hand, the picture in NA is different: the closest related strain to that of the Mexico outbreak is a subtype H2N3 AIV isolated from a green winged teal in Illinois in 2010 ( Figure S2 ).

Sequences for time-scaled phylogenetic analysis were selected from the closest clades to the novel H7N3 HPAI viruses on the maximum likelihood tree of each segment. This dataset comprised 427 AIV strains collected over a 12 year period (2001 to 2012). The time to the most recent common ancestor (TMRCA) for each segment of the novel HPAI H7N3 in Mexico was estimated from the time-scaled phylogenetic trees. The HPAI H7N3 strains sampled in Mexico shared similar common ancestors among different genes between October 2011 and March 2012, i.e. during the winter of 2011–2012 ( Table 1 ). The common ancestor of the HPAI H7N3 Mexico outbreak and the closest related avian influenza strains existed between 1.1 to 3.9 years ago, which varied among their different genomic segments ( Table 1 ). The difference in unsampled diversity among gene segments suggested that the reassortment of North American AIV lineages which led to the H7N3 Mexico outbreak may have involved several events spread over this time period. This can be seen by comparing the closest related strains in the phylogenetic tree for any segment: they can be quite distant from the H7N3 Mexico strain in the other segments. This result supports our hypothesis of the occurrence of multiple reassortment events.

Diverse reassortment events involving the six internal segments can be inferred from the maximum likelihood phylogenies of 2343 North American AIV. Clades identified in the phylogeny for one segment (e.g., PB2) are not maintained in the phylogenies of other internal segments ( Figure S1 ). In addition, internal segments of AIV viruses isolated in distant locations can be closely related to each other within the same time period, which suggests not just frequent reassortment but also rapid movement of influenza viruses across North America.

A: HA. Sequences in grey are from before 1990; the clade colored blue is composed of H7N2 AIV isolated from a single surveillance in poultry in New York; the clade colored pink was selected for time-scaled phylogenetic analysis. B: NA tree. The uncoloured sequences are from before 2000; the AIV clade in pink was selected for time-scaled phylogenetic analysis.

To investigate the origin of the AIV causing the HPAI H7N3 outbreak in Mexico in 2012, an initial phylogenetic analysis using Maximum likelihood was performed for each segment of both the outbreak sequences and a background dataset which comprised all available AIV of North American AIV lineages ( Figure 1 ). The phylogenetic trees of all available H7 segments in north American showed that AIV isolated in recent years have diverged from those before 1990 ( Figure 1A ). In addition, in the HA segment a sub-lineage mainly composed of H7N2 AIV from domestic birds in New York state is clearly separate from the recent lineage composed of AIV from wild birds, which indicates extensive diversity of LP AIV in wild and domestic birds. Since 2000, the N3 NA segment of North American AIV has split into two separate lineages ( Figure 1B ). The mechanism for maintenance of this divergence remains unknown as viruses from both lineages co-circulate in geographically overlapping host populations, mainly wild waterfowl.

Gene flow of the precursor of the HPAI H7N3 outbreak in Mexico

To further explore the origin of the Mexico outbreak strain, a joint analysis of discrete trait models was performed to estimate the overall genetic transmission process. In this the phylogenetic tree space was sampled independently for each segment, while the transition pattern was jointly estimated in a single analysis as the diffusion parameters being applied in the discrete trait models were assumed to be the same (see methods). Four major factors including seven specific traits were tested by implementing Bayesian stochastic search variable selection (BSSVS): i) host population of AIV (order/species); ii) geographic location of sampled AIV (bird migration flyways/provinces and states of North America); iii) subtype of AIV and iv) virulence (pathogenicity/cleavage sites). For each trait, the evolving process of the HPAI H7N3 in Mexico and closely related AIV can be seen from the reconstructed time-scaled phylogeny of each segment independently (with exception of pathogenicity and cleavage sites which only applied to the HA segment), with the branches colored by the specific trait according to the ancestor trait in the internal nodes (Figure 2 and Figure S2, S3, S4, S5, S6, S7, S8).

Host populations of AIV in this study were first analysed by Order: wild Charadriiformes (cha) such as gulls, wild Gruiformes (gru) such as cranes, wild Passeriformes (pas) such as sparrows, domestic Galliformes (gal) such chickens and wild Anseriformes (ans) such as mallards, with Anseriformes comprising the majority of our AIV data set (n = 366/427), see table S1. Among all four strongly supported transitions with Bayes Factor (BF) >100 with mean diffusion rate (R) between 0.01 and 0.07, the highest diffusion rate was found between Charadriiformes and Passeriformes. The other three were found between Anseriformes and other bird orders, and the HPAI H7N3 outbreak in poultry (labelled as Galliformes Mexico) is linked to Anseriformes with strong support (R = 0.02, BF>100) (Figure 3A and Table S5). The results confirm that there has been extensive mixing of influenza A virus between different orders of birds, both wild and domestic.

PPT PowerPoint slide

PowerPoint slide PNG larger image

larger image TIFF original image Download: Figure 3. Inferred host transmission networks of Mexican outbreak AIV. A: Host order. Node labels in the nodes are host orders identified following the abbreviations used in the colored phylogenetic trees (Figure 2 and Figures S2, S3, S4, S5, S6, S7, S8): wild birds of the order Anseriformes (ans-wild); wild birds of the order Charadriiformes (cha-wild); wild birds of the order Passeriformes (pas-wild); domestic birds of the order Galliformes and Mexico H7N3 outbreak in the order Galliformes (gal-domestic-Mexico). Arrows show the direction of transmission between two host orders; the arrow weight and the number above each arrow indicates the per capita transmission rate. Node size reflects the number of AIV for each host order (Table S1). Line colours indicate the overall Bayes Factor test support for epidemiological linkage between host orders, Red lines indicate statistical support with BF>100 (very strong support), dark pink lines indicate support with 30<BF<100(strong support), pink lines indicate support with 3<BF<30. B: Host species (Anseriformes only). Wild Anseriformes are further classified into the five main species and a group comprising the other rarer species of Anseriformes in this study: mallard (Anas platyrhynchos), northern pintail (Anasacuta), northern shoveller, blue-winged teal, green-winged teal and other Anseriformes (other ans); As above, arrows show the direction of transmission between two host species; the arrow weight and the number above each arrow indicates the per capita transmission rate. Node size reflects the number of AIV for each host species (Table S3). Line colours indicate the overall Bayes Factor test support for epidemiological linkage between host species, Red lines indicate statistical support with BF>100 (very strong support), dark pink lines indicate support with 30<BF<100(strong support), pink lines indicate support with 3<BF<30. https://doi.org/10.1371/journal.pone.0107330.g003

To explore which host species might have been the direct donor of the Mexico outbreak strains, AIV belonging to Anseriformes were further divided into the five predominant species: mallard (Anas platyrhynchos), northern pintail (Anasacuta), northern shoveller (Anas clypeata), blue-winged teal (Anas discors) and green-winged teal (Anas carolinensis) (Table S2). Species that were sampled at relatively low levels were combined as “other Anseriformes” (Table S2). Multiple statistically supported transitions (with R from 0.15 to 2.13, BF from 6 to over 100) were identified among different host species within this Order, and both mallard and green winged teal are linked to 3 other host species (Figure 3B and Table S6). This analysis indicates the Mexican outbreak strains were most likely to have been transmitted from green winged teals (R = 0.15 and BF = 6).

The phylogeographic analysis for each segment of the Mexico HPAI H7N3 strain was summarised by a MCC tree in a geographic context. However, to visualize the evolution process in a spatiotemporal mode we converted the spatial annotated time-scaled phylogeny to an annotated map (Figure 4 and Figure S9). Five segments (HA, NA, NP, M, and NS) of the Mexico HPAI H7N3 strain were introduced directly from different states in central US, while PB2, PB1 and PA were introduced from states in the western region. The introductions of segments from several different geographic locations indicate multiple reassortment events were likely to have been involved in the generation of the novel H7N3 Mexico AIV.

PPT PowerPoint slide

PowerPoint slide PNG larger image

larger image TIFF original image Download: Figure 4. Spatial diffusion of AIV segments of the Mexico outbreak AIV. The first three panels represent three segments separately (A: HA, B: NA, C: PB2) and D represents the spatial transmission of all 8 segments jointly. The plotted lines represent the branches of the MCC trees for different segments, distinguished by color; the size of each circle represents the number of lineages with that location state. The map source for this figure was OpenStreetMap (http://www.openstreetmap.org/). The spatial diffusions of other five segments (PB1, PA NP, M and NS) on the map are shown in Figure S9. https://doi.org/10.1371/journal.pone.0107330.g004

Joint discrete trait analysis of all eight segments indicated frequent gene transfer among locations (states and provinces) where the background AIV sequences were isolated. However, in this initial analysis, no significant support was found between Jalisco (the outbreak state) and any other location, probably due to the large number of possible transitions (26 states, 325 irreversible transition pairs) and the limited number of outbreak strains (3; Table S3). Previous studies have shown that incorporating location greatly improves phylogeographic descriptions of the pattern of virus gene flow [16], [17]. Therefore, we enhanced the statistical power of the analysis in two ways: first by decreasing the number of locations by combining states and provinces into major regions (flyways) and secondly by reducing the number of pairwise transitions possible.

To aggregate locations we used the known migration routes, or “flyways”: Atlantic, Mississippi, Central, or Pacific (Figure S10). The distributions of avian influenza virus for each flyway are summarized in Table S4, showing a wide range in rate and statistical support (R = 0.02 to 0.25; BF = 3 to >100). Highly significant links [BF>100, Indicator (I) = 1] were found between major North American flyways, particularly between Atlantic and Mississippi (R = 0.17 exchange/year); Mississippi and Pacific (R = 0.22 exchange/year), Central and Atlantic (R = 0.25 exchange/year) and Central and Pacific (R = 0.18 exchange/year) (Figure 5A and Table S7). Linkages between Atlantic and Pacific flyways, Mississippi and Central flyways were also identified although with weaker support (3<BF<6). The transitions between flyways and the H7N3 HPAI in Mexico are not that strongly supported compared to the between flyway transitions, probably due to the limited number of sequences available from the outbreak AIV, but three direct donors among these flyways to the predecessor of the Mexico outbreak were identified: the Pacific, Central and Mississippi flyways. Among these flows, the transition rate from Mississippi (R = 0.05 exchange/year) is the most strongly supported with BF = 86. In comparison, the link with the Pacific (R = 0.02 exchange/year) and Central (R = 0.04 exchange/year) were weaker (BF = 4; Figure 5B and Table S7). There was no supported link between the Atlantic flyway and Mexico. The results indicated that the HPAI H7N3 in Mexico probably originated from AIV transmitted by wild birds from three different flyways.

PPT PowerPoint slide

PowerPoint slide PNG larger image

larger image TIFF original image Download: Figure 5. Inferred phylogeographic transmission networks of Mexican outbreak AIV. A: Flyway. AIV transmission among 4 N. American flyways with links to the Mexican outbreak strains. Arrows show the direction of transmission between two flyways; arrow weight and the number above each arrow indicates the per capita transmission rate. Node size reflects the number of AIV for each flyway (Table S6). B: Location. AIV transmission among states/provinces in North America and Jalisco (the Mexican state where the outbreak strains were isolated). Arrows show the direction of transmission between the two states; the arrow weight and the number above each arrow indicates the per capita transmission rate. Node size reflects the number of AIV for each flyway (Table S5). https://doi.org/10.1371/journal.pone.0107330.g005

Given these findings it appeared likely that the precursor strains were generated somewhere near Mexico as it is the place where birds from the different migration routes meet during winter. To test this hypothesis, we further reduced the number of transitions in the locations transition matrix: transitions between two flyway regions were switched off (by forcing the initial indicator of a given transition pair from 1 to 0, so that this transition pair will not be counted), and those for within flyway transitions and transitions linked to Mexico were maintained. There are 98 non-reversible transition pairs in the new reduced matrix (Table S10). AICM tests (see Methods) revealed that the non-reversible BSSVS model with reduced number of transitions was significantly favoured over the other models with the original matrix (Table 2), indicating the number of transitions has an effect on the performance of discrete trait models. This reduced model has better support than a randomly reduced model with the same number of transition pairs and the same non-reversible BSSVS setting (Table 3); we conclude that gene transitions within flyways and Mexico alone better explain the gene flow of North America AIV than a model incorporating the between flyway transitions.

Considering individual locations, we found 11 locations were linked to the HPAI H7N3 in Mexico, among which 7 showed significantly strong links (BF>100). These were: Alaska (R = 0.59), Alberta (R = 1.07), California (R = 0.59), Illinois (R = 1.61), Missouri (R = 2.26), Ohio (R = 0.66) and Wisconsin (R = 1.39), belonging to the Pacific, Central and Mississippi flyways. AIV from other four states/provinces (Minnesota, New Brunswick, Quebec and Washington) are also linked to Mexico H7N3 AIV but with weaker support (BF = 21 to 30) and lower rate (0.7 to 0.84) (Figure 5A and Table S8). This result indicates the donor locations of the Mexico outbreak are spread widely across North America.

In addition, an extremely complex pattern of linkage between the 52 ancestral subtypes (Table S9) was identified, which confirmed the extent of the reassortment events which had occurred between, especially in the internal segments (full tree with annotation can be found on the Dryad Digital Repository: doi:10.5061/dryad.j5bf8). However, as for locations, no significant linkage between H7N3 in Mexico and other subtypes was identified due to the large number of candidate donor subtypes.

Mexico H7N3 is confirmed in the phylogenetic trees of the HA segment (Figure S11) to have mutated from a LPAI (low pathogenic avian influenza) virus to HPAI after the ancestral virus was introduced into poultry from wild birds (Figure S11A), rather than being associated with previous HP outbreaks. The virulence of HP avian influenza viruses is associated with the appearance of an insertion of multiple basic amino acids at the cleavage site of the HA protein [18]. Categorizing the cleavage sites in this dataset into three types: 1) Insertion, 2) Partial insertion,3) No insertion (Figure S11 C), we find that the H7N3 Mexico strain has a unique insertion - DRKSRHRR - compared to other HPAI strains (Figure S11 C). We conclude the Mexico H7N3 strains originated from a lineage composed of LP strains with partial insertions in the cleavage sites. Similarly, the H7N3 strains which caused an outbreak of HPAI in Canada in 2004 were also derived from LP strains with the same partial insertions, but from a completely separate HA clade, indicating parallel evolution with respect to the acquisition of the multibasic cleavage site, starting from different lineages (Figure S11 B).