During 2001–2008, as part of the Vietnam National Influenza Surveillance System, nasopharyngeal or throat swab samples were collected from patients seeking care for ILI at hospitals (6). Specimens were tested for influenza A and B viruses and were further subtyped for H1, H3, and H5 by reverse transcription PCR by using primers, probes, and reagents recommended by the Centers for Disease Control and Prevention and the World Health Organization (WHO). Samples that were positive for influenza A by PCR were selected for virus isolation, and isolates reaching a titer of 1:8 in hemagglutination assays were selected for sequencing analysis. All isolates were subtyped by using hemagglutination assays with reference antigens and antiserum from the WHO reagent kit. A total of 242 samples were shipped to the National Institutes of Health Influenza Genome Sequencing Project (USA) (22) for whole-genome sequencing at the J. Craig Venter Institute. Of the 242 samples, 2 were excluded from this analysis (1 that could not be sequenced and 1 from a patient with a mixed infection). The final dataset of the 240 whole-genome sequences comprised 145 influenza subtypes H3N2 and 95 H1N1 (GenBank accession nos. CY103972–CY105893). Table 1 shows the numbers and locations of the viruses.

For phylogeographic analysis, we compiled influenza virus sequences of subtypes H1N1 and H3N2 into 2 datasets: a regional dataset of whole-genome sequences from Asia and Australia/New Zealand and a global dataset of geographically subsampled sequences (50 replicates) of the hemagglutinin (HA) segment (Technical Appendix). For subsampling, we randomly sampled 12 sequences per geographic region per year.

Sequences were aligned by using the MUSCLE program, version 3.8 (23). Maximum-likelihood trees were inferred by using RAxML version 7.3.0 with 2,000 bootstrap replicates (24,25). For the regional HA datasets, phylogenetic trees with sampling date information were inferred by using BEAST version 1.6.2 (26) and a relaxed molecular clock (uncorrelated lognormal). The nucleotide substitution model was SRD06 + HKY85 + Γ, and the demographic models used were constant population size and Bayesian skyride (Technical Appendix).

Migration analysis was conducted by using a straightforward parsimony method in the PAUP* program (27,28). The 2,000 bootstrap trees and the best maximum-likelihood tree inferred by RAxML were read into PAUP*, and nucleotide sequences were replaced by single-letter location codes assigned to a set of predefined global regions (Technical Appendix Figures 5, 6). Changes in location code were mapped onto the branches of the trees by using an ACCTRAN parsimony criterion (28). Analysis was performed on all segments to determine whether migration histories differed among them, which could have been caused by reassortment among influenza virus RNA segments. A strong reassortment signal was verified when standard mosaic and phylogenetic criteria were used (29) (Technical Appendix Figure 2).

For comparison, we performed the same migration analysis on the global dataset of 50 subsampled replicates for influenza virus subtypes H3N2 and H1N1. Migration matrices were built describing numbers of connections between 27 subtype H3N2 or 29 subtype H1N1 predefined geographic regions. We used Gephi software (30) to visualize the connections in the matrix. For regions with sufficient samples, we computed minimum distances to the trunk of the rooted phylogeny for all 50 subsampled trees (subtype H3N2 only) to determine whether viruses from different regions could be described as ancestral (close to the trunk) or derived (far from the trunk).

The relationship between the subtype H3N2 HA sequences from Vietnam and other viruses sampled in the region is shown in Figure 1. Representative samples from Vietnam are available for 2003–2008 but not for 2006 when subtype H1N1 predominated. Viruses isolated in Vietnam show close relationships to viruses isolated in Hong Kong, Taiwan, Singapore, and Australia/New Zealand. Inferred from this tree were 20 parsimony-unambiguous migration events, 9 showing Vietnam–Hong Kong migration, 7 showing Vietnam–Australia/New Zealand migration, 2 showing Vietnam–Taiwan migration, 1 showing Vietnam–South Korea migration, and 1 showing Vietnam–Qatar migration. Clearly, because Hong Kong (68 sequences) and Australia/New Zealand (>500 sequences) were overrepresented in the regional dataset, most Vietnam migrations were associated with these 2 locations. Because we were initially uncertain how well a geographically regional phylogenetic tree would reflect the true migration patterns of influenza viruses in Vietnam, we performed a validation exercise to determine what proportion of global migration of influenza subtype H3N2 virus from Vietnam could be detected in a regional phylogenetic tree (Technical Appendix). Approximately 70% of Vietnam influenza migrations from a global analysis were also observed in the regional tree.

The relationship between influenza subtype H1N1 HA sequences from Vietnam and other regional viruses is shown in Figure 2. Representative samples from Vietnam are available for 2001–2008 but not for 2004 and 2007. Inferred from this tree are 10 parsimony-unambiguous migration events, 6 showing Vietnam–Taiwan migration, 2 showing Vietnam–Japan migration, and 2 showing Vietnam–Australia/New Zealand migration. As in the analysis for subtype H3N2, these migration links correspond with the viruses that were sequenced from the region during 2001–2008.

Because the regional trees included only sequences for which whole genomes were available, migration patterns were compared systematically across all 8 influenza segments. Because influenza viruses reassort, different event histories should be visible in phylogenies inferred separately for the 8 virus segments. Indeed, for the subtype H3N2 dataset, we observed a median of 14 parsimony-unambiguous migration events for the neuraminidase segment and a median of 41 for the matrix protein segment; the other segments fell somewhere in between (Table 2). Again, most migrations were with Australia/New Zealand and Hong Kong, indicating that the pattern of migration is similar across segments, although different numbers of migrations and different individual migration events are visible when different segments are analyzed. For the matrix protein and nonstructural protein segments of subtype H3N2 viruses, the large number of migration events may result from the larger number of topologically uncertain and polytomic nodes in these trees that had to be randomly resolved to compute the number of migration events; that is, sequences from one country could be mistakenly mixed with sequences from other countries, thus generating some artificial migration events in the parsimony analysis. The low confidence in the Vietnam–Singapore migration link for subtype H3N2 may result from the small number of whole-genome sequences available from Singapore, all of which were collected in 2003.

For subtype H1N1 viruses, we observed 6–16 migration events across the trees inferred for the 8 segments (Table 2). In the bootstrapped data, the Vietnam–Taiwan and Vietnam–Japan migratory connections seem to be approximately equal, despite the fact that the best maximum-likelihood tree showed 6 Vietnam–Taiwan connections and 2 Vietnam–Japan connections. The migratory connection between Vietnam and Australia/New Zealand seems to be somewhat weaker, possibly because of substantially less sampling of Australia/New Zealand viruses in the subtype H1N1 dataset. In the subtype H1N1 dataset, the number of migration events for the matrix protein and nonstructural protein segments did not increase.

Because it is clear that the presence and number of samples from different regions influence migration analysis, migratory patterns were reanalyzed on the global subsampled dataset to reduce the geographic bias present from having higher numbers of samples available from some regions than others. Using global HA trees for subtype H3N2 and H1N1 sequences, we constructed full migration matrices including all parsimony-unambiguous migration events among our predefined regions (Technical Appendix Figures 5, 6). These migration networks are shown in Figure 3, where the United States is a major hub of influenza migration and eastern Asia and Australia/New Zealand play major roles. The subtype H3N2 data show Vietnam connected with most other countries in eastern and Southeast Asia, with the United States, and weakly with southern Asia. The subtype H1N1 data show Vietnam connected with the United States and Europe but weakly with other Asian countries. For both subtypes, the total number of migration events associated with each node in the network is correlated with the number of samples available for that node (all p values were <10⁻⁵; Kendall and Spearman tests). Note that sample numbers are not identical for each node because for some regions <12 sequences per year were available, and these regions did not need to be subsampled for those years. Hence, undersampling and oversampling can generate this correlation. Despite our attempt to reduce geographic bias in the global dataset, inference on migration events is still closely associated with regional availability of samples; this bias appears to affect all phylogeographic studies.

For H3N2 sequences, to determine whether any region has the characteristics of an ecological source, we computed the phylogenetic distance of sequences from 6 well-sampled regions (China, Hong Kong, Japan, Vietnam, Australia/New Zealand, and the United States) to the trunk of the global maximum-likelihood phylogenetic tree (Figure 4). In 2003, for example, across all 50 subsampled trees, sequences isolated in China were typically closest to the trunk of the phylogenetic tree, indicating that these sequences are ancestral to other viral sequences sampled in 2003; this finding is consistent with the global replacement of subtype H3N2 viruses by the A/Fujian/411/2002-lineage that occurred in 2003. In general, for the years 2003–2007, in no region were sequences consistently ancestral, indicating that it is unlikely that there is a single global source of human influenza viruses. The more likely global migration model involves periodic global strain replacements originating in different regions in different years (3,4). There were not sufficient samples from all regions/years to perform this analysis on the subtype H1N1 dataset.

Figure 5 shows a Bayesian subtype H3N2 phylogenetic tree inferred from the time-stamped regional sequence data. The insets of this figure detail the 2007–2008 part of the phylogeny (87 sequences) and the coalescent times for the tips of these branches. It is difficult to draw a complete persistence picture for these viruses because of undersampling during the second quarter of 2007 and the first half of 2008 despite PCR-confirmed evidence of influenza virus activity during these periods (6). The median coalescent time for viruses from Vietnam sampled during this period is 37 days (interquartile range 21–72 days), and the insets in Figure 5 suggest that one of the lineages persisted in Vietnam for the 13 months from January 2007 through January 2008. An absence of samples from February through May 2008 makes it impossible to determine conclusively if this lineage persisted in Vietnam for the entire 2-year period.

To determine whether the lack of samples from other countries created an artifactual picture of lineage persistence in Vietnam during 2007–2008, we assembled a sequence set of all 672 viruses from Asia and Australia/New Zealand from 2006 through 2008. The maximum clade credibility tree of these sequences (Technical Appendix Figure 4) indicates that the Vietnam lineages separate into >10 distinct lineages when viewed in the context of all Asian/ Australia/New Zealand influenza viruses. One of these lineages persisted in Vietnam for 15 months (Technical Appendix Figure 1 panel A), and another persisted for 10–12 months (Technical Appendix Figure 1, panel B), suggesting that lineage persistence of ≥1 year may have occurred in Vietnam during 2007–2008. However, this type of analysis is very sensitive to phylogenetic uncertainty because the individual lineages (or subclades) contain few sequences and may not be robust to small changes in tree topology.