Bayesian Phylogeography Finds Its Roots

Introduction

Phylogenetic inference from molecular sequences is becoming an increasingly popular tool to trace the patterns of pathogen dispersal. The time-scale of epidemic spread usually provides ample time for rapidly evolving viruses to accumulate informative mutations in their genomes [1]. As a consequence, spatial diffusion—among other processes—can leave a measurable footprint in sampled gene sequences from these viruses [1]. Reconstructing both the evolutionary history and spatial process from these sequences provides fundamental understanding of the evolutionary dynamics underlying epidemics, e.g. [2],[3]. It is also hoped that these insights can be translated to effective intervention and prevention strategies [4] and elucidating the key factors in viral transmission and gene flow over larger distances is central in formulating such strategies, e.g. [5].

Phylogeographic analyses are a common approach in molecular ecology, connecting historical processes in evolution with spatial distributions that traditionally scale over millions of years [6]. Many popular phylogeographic approaches [7],[8] can be remiss in ignoring the interaction between evolutionary processes and spatial-temporal domains. One first reconstructs a phylogeny omitting spatial information and then conditions the phylogeographic inferences on this reconstruction [1],[9], exploiting non-parametric tests to evaluate the significance of this conditional structure, e.g. [7],[10],[11]. To draw conclusions about the epidemic origin or epidemiological linkage between locations, however, we require a reconstruction of the dispersal patterns and process throughout the evolutionary history. Considering locations as discrete states, this boils down to the well-known problem of ancestral state inference [7]. Parsimony is a popular heuristic approach to map characters onto a single phylogenetic tree [12]. Unfortunately, parsimony reconstructions ignore important sources of model uncertainty, including both uncertainty in the dispersal process as well as in the unknown phylogeny [13]. In addition, minimizing the number of state exchanges over a phylogeny is misleading when rates of evolution are rapid and when the state exchange probabilities are unequal [14].

Probabilistic methods draw on an explicit model of state evolution, permitting the ability to glimpse the complete state history over the entire phylogeny and conveniently draw statistical inferences [15]–[17]. These analyses typically employ continuous-time Markov chain models for discrete state evolution analogous to common nucleotide, codon or amino acid substitution models [18]. In contrast to parsimony, maximum likelihood-based reconstructions incorporate branch length differences in calculating the conditional probability of each ancestral state given the observed states at the phylogeny tips [14]. Bayesian reconstruction methods enable further generalization of this conditional probability analysis by removing the necessity to fix the Markov model parameters to obtain ancestral states and the necessity to specify a fixed tree topology with known branch lengths. Bayesian inference integrates conclusions over all possible parameter values but to achieve this, however, requires prior probability distributions for all aspects of the model.

While probabilistic methods have been previously presented in a bio- or phylogeographic context, in particular Bayesian methods that integrate over phylogenetic uncertainty and Markov model parameter uncertainty [19], viral phylogeography studies have rarely made use of these developments. This may be a consequence of low awareness of existing software implementations for arbitrary continuous-time Markov chain models [20],[21] or a lack of appreciation for the uncertainty intrinsic in these reconstructions and the ease with which one can formally access epidemiological linkage through probabilistic approaches. A recent phylogeographic study of influenza A H5N1 introduces a heuristic non-parametric test to evaluate whether parsimony-inferred migration events between two particular locations occur at significantly high frequency [22]. Null distributions for these frequencies arise from randomizing tip localities after false discovery rate correction to control for simultaneous testing issues. Although this procedure addresses concerns about statistical inference on sparse frequency matrices, the multiple comparison correction still results in a conservative estimate of significant migration events. Fully probabilistic approaches may further ease statistical inference, yet similar tests remain lacking for likelihood-based phylogeographic models.

Advances in evolutionary inference methodology have frequently demonstrated how novel approaches can be appended to a sequence of analyses, in many cases starting from alignment to parameter estimation conditional on tree reconstructions. For example, demographic inference has involved genealogy reconstruction, estimating a time scale for the evolutionary history, and coalescent theory to quantify the demographic impact on this tree shape [23]. It is well acknowledged that such sequential procedures ignore important sources of uncertainty because they generally purge error associated with each intermediate estimate. With the advent of novel computational techniques like Markov chain Monte Carlo (MCMC) sampling, it has become feasible to integrate many of the models involved and simultaneously estimate parameters of interest. Demographic inference is a well-known example of genealogy-based population genetics that benefited from these advances [24],[25]. Bayesian MCMC methods also enable ancestral state reconstruction while simultaneously accounting for both phylogenetic and mapping uncertainty. Although this adds much needed credibility to ancestral reconstruction [13], phylogeographic analysis would benefit even more from fully integrating spatial, temporal and demographic inference.

Here, we implement ancestral reconstruction of discrete states in a Bayesian statistical framework for evolutionary hypothesis testing that is geared towards rooted, time-measured phylogenies. This allows character mapping in natural time scales, calibrated under a strict or relaxed molecular clock, in combination with several models of population size change. We use this full probabilistic approach to study viral phylogeography and extend the Bayesian implementation to a mixture model in which exchange rates in the Markov model are allowed to be zero with some probability. This Bayesian stochastic search variable selection (BSSVS) enables us to construct a Bayes factor test that identifies the most parsimonious description of the phylogeographic diffusion process. We also demonstrate how the geographical distribution of the sampling locations can be incorporated as prior specifications. Through feature-rich visual summaries of the space-time process, we demonstrate how this approach can offer insights into the spatial epidemic history of Avian influenza A-H5N1 and rabies viruses in Africa.

The highly pathogenic avian influenza A-H5N1 viruses have been present for over a decade in Southern China and spread in multiple waves to different types of poultry in countries across Asia, Africa and Europe [26]. As a result, highly pathogenic A-H5N1 is now a panzootic disease and represents a continuous threat for human spill-over. Strong surveillance has been in place since these viruses caused extensive outbreaks, but the source and early dissemination pathways have remained uncertain. Because parsimony analysis has attempted to shed light on the latter [22], A-H5N1 provides an ideal example for comparison with Bayesian phylogeographic inference. Rabies is endemic in Asia and Africa, where the primary reservoir and vector for rabies virus (RABV) is the domestic dog. Phylogenetic analysis has revealed several genotypes of lyssaviruses (family Rhabdoviridae); genotype 1 has been found responsible for classical rabies, a fatal disease in terrestrial mammals throughout the world [27],[28]. Here, we explore the phylogeographic history of RABV in domestic dogs in West Central Africa, using recently obtained sequence data, and evaluate the role of viral dispersal in maintaining RABV epidemic cycles.

Results

We examine the evolution and spatial dispersion of two viral pathogens, Avian influenza A-H5N1 and rabies, to demonstrate the strengths and limitations of our discretized stochastic model for phylogeography.

Discussion

The Bayesian phylogeographic inference framework we present here incorporates the spatial and temporal dynamics of gene flow. In this study, we focus on pathogen diffusion because viral sequence sampling on a time-scale commensurate with the rate of substitution permits the inference of spatial patterns in real-time units. In addition, elucidating the phylodynamics of viral epidemics has important implications for public health management. We selected the Avian influenza A-H5N1 example to allow a convenient comparison of Bayesian ancestral state inference with the previous parsimony analysis; on the other hand, statistical analysis of the rabies migration in Africa up to this point has been largely unexplored. Both zoonoses represent a clear threat to human health. The frequent transmission of A-H5N1 from poultry or wild birds to humans suggest that the virus could emerge as or contribute genetically to the next human flu pandemic. Although the lack of a human-to-human transmission mechanism means that rabies will not emerge as a purely human disease, rabies infection causes a fatal neurological disease and at least 55,000 people die from this disease every year, mainly in the developing world [40].

A Bayesian statistical approach presents many advantages over parsimony inference of ancestral states. First, MCMC offers a computational technique to integrate over an unknown phylogeny and unknown migration process as the former is not directly observable in nature and the latter is poorly understood. Accommodating this lack of knowledge protects against potentially severe bias, but can reduce the power to make inferences; phylogeographic analyses are no exception to this. One can regard this uncertainty itself as a ‘mixed blessing’ because whilst it can hamper drawing definitive conclusions [13], it protect us from making overstated conclusions. For example, parsimony analysis of the influenza data establishes an epidemiological link between Guangdong and Indonesia [22]. Bayesian inference does not confirm this conclusion and phylogenetic analysis of more recently obtained sequence data now identifies the progenitors of Indonesian strains in the Chinese province of Hunan [26], a site which our shared analysis also identified as a possible source. Further, unlike parsimony, likelihood-based probabilistic methods consider branch lengths in ancestral reconstructions. The impact of the tree depth on root state reconstructions for the A-H5N1 genome segments clearly illustrates the importance of branch lengths. Moreover, probabilistic methods allow for estimating the relative posterior probability of each location state at any position along the phylogeny; this ability is indispensable in a hypothesis testing framework. As introduced by [19] in a different setting, we also demonstrate how phylogeographic parameters can be estimated from different genomic segments without assuming the same evolutionary history. H5N1 reassortment, however, will have not have fully unlinked HA and NA evolutionary histories and partially shared ancestry may lead to overstated credibility in some aspects of the phylogeographic inference.

Bayesian inference also proffers particular benefits within the class of likelihood-based methods, for example, by allowing for straightforward approaches to control model complexity. BSSVS naturally provides a BF test to identify significant non-zero migration rates. Further prior specification easily incorporates geographical detail of the sequence data. Although distance-informed priors appear to have little impact on the phylogeographic analyses presented here, both BSSVS and informed priors furnish new opportunities for hypothesis testing when comparing competing prior scenarios of the diffusion process. Examples include “gravity models” [41] in which infinitesimal rates become functions of the host population-sizes at the end-point locations and a priori structurally-fixed graphs [19]. Finally, it has been recognized that an MCMC-based Bayesian framework is well-suited to bring together information of different kinds [42]. The BEAST software, which has a strong focus on calibrated phylogenies and genealogies, elegantly illustrates this by offering a large number of complementary evolutionary models including substitution models, demographic and relaxed clock models that can be combined into a full probabilistic model [43]. By adding spatial reconstruction to this arsenal of evolutionary models, the full probabilistic inference now brings us much closer to biogeographical history reconstruction from genetic data.

Our primary motivation for exploiting BSSVS to select among all possible migration graphs is to elucidate the limited number of epidemiological links that appropriately explain the viral diffusion process. This parsimonious set both informs major modes of migration and reduces the high statistical variance that burdens estimation of all pairwise transition rates. Following this argument, less uncertainty on node state reconstructions would be expected when focusing on a parsimonious parameterization of the instantaneous rate matrix. The A-H5N1 analysis indeed indicates lower uncertainty of root state reconstructions. However, for some other internal nodes, we note the opposite behavior. We attribute this to the reversibility assumption in the rate matrix. Selection of reversible rates by BSSVS imposes more balanced transitions in the phylogeny among locations that could have unidirectional links in reality. Therefore, work is in progress to develop non-reversible models that may better fit a spatially expanding epidemic like A-H5N1 or recurring epidemic influenza emergence through source-sink dynamics [2]. Considerable technical hurdles remain to incorporate BSSVS procedures under such models. Because BSSVS places non-negligible probability on structural zeros in the rate matrix, we can not guarantee that all resulting rate matrices are diagonalizable, challenging stable computation. Bearing in mind the reversibility assumption, we pass no judgement on the origin of the A-H5N1 epidemic based on the frequency by which a location is present in well-supported rates, as was previously done in the parsimony analysis [22]. Instead, we focus on node location state reconstructions throughout the phylogeny and their posterior probabilities. Figure 1 suggests that, although Hong Kong and Guangdong both receive posterior support as root location states, the dominant location throughout the phylogeny and hence the more likely hub of diffusion is Guangdong. An inherent assumption of the discrete model of location change is that ancestral viruses necessarily reside at only the sampled locations of the extant viruses. In this respect, it is important to realize that the CTMC process should describe the underlying spatial dynamics more accurately as the sampling density increases. E.g., for A-H5N1 [26], provide more recently obtained sequence data across a larger set of geographic locations; the data could inform further pathways seeding remote localities that remain elusive in our present analysis. In addition to tackling the reversibility assumption, it may also prove necessary to relax constant diffusion rates through time to realistically model phylogeographic processes in many situations. Covarion-like models [44] and allowing different diffusion matrices across different time intervals in the phylogeny may help achieve these aims.

Our rabies phylogeographic analysis confirms a longstanding presence of this viral lineage in West Africa [28]. The virus appears to have a constant population size for about 150 years during which, extrapolating from the more recent spatial dynamics, diffusion occurs continuously with no particular directionality (Figure 6A). These continuous dynamics explain why we can not achieve precise root state location inference based only on samples from the last 20 years. In the light of the constant population dynamics, however, the location of the MRCA may be epidemiologically irrelevant as the location probably does not necessarily represent the ultimate source of the rabies endemic. We note that our analysis does not include all currently available strains originating from Chad, which may add to the weak East-to-West dispersal signal revealed by a recent parsimony analysis [39]. Our analysis confirms the model proposed for dog RABV in general; that is, of a series of spatially distinct clusters that experience relatively little contact among them [28],[39]. By providing a time scale for the seeding of these spatial clusters, we again demonstrate a clear advantage of the Bayesian inference over parsimony analysis. The ability to draw migrations over time also promotes a more precise dissection of local and temporal RABV movement on smaller geographical scales. After migration, the virus appears to establish local populations maintaining the viral lineage for at least a limited amount of time. These dynamics are reminiscent of a metapopulation model with continuous turnover of locally established viral populations. A long-standing rabies presence in West Africa is not surprising; already recognized since the late 60s, the territory plays a major role in the rabies-canid ecological balance [45]. It remains a remarkable feat that an acute and mainly fatal disease achieves prolonged endemicity. Because disease-induced mortality can rapidly deplete the number of susceptibles in a population, one expects epidemic cycles with oscillatory dynamics to occur. Rabies cycles and traveling waves have been well documented in wildlife across Europe and North America, e.g. [46], and more recently [47], demonstrate such cycles in African dogs. Because their periodicity is notably shorter than expected from epidemiological models, the authors argue that intervention responses also impact the epidemic cycles [47]. Importantly, there is also a remarkable phase synchrony in rabies outbreaks across southern and eastern Africa, most pronounced for distances up to 1,000 km [47]. For oscillating systems in particular, it is well known that dispersal can generate population synchrony [48]. Because previous studies illustrate that even limited amounts of relatively local dispersal can generate synchrony in cyclical dynamics over large spatial scales [48], and that the resulting synchrony tends to decline as distance increases and varies through time [49], [47] argue that dispersal could enforce synchrony in dog rabies epidemics across different countries. Our analysis clearly reveals rabies dispersal as a continuous dynamic process that could indeed be essential in maintaining epidemic cycles. As argued by [39], however, the rate of dispersal is probably not sufficiently high to explain the short epidemic cycles as suggested by [47]. Nevertheless, we underscore that sustained and coordinated responses across political boundaries are necessary to control domestic dog rabies in Africa.

Many questions in evolutionary biology require a biogeographical perspective on the population under investigation. We hope to have demonstrated that Bayesian phylogeographic framework can contribute significantly to evolutionary hypothesis testing, and, although we have focused on viral phylodynamics, this approach is generally applicable in molecular evolution. Employing geographically-informed priors delivers a first step in incorporating GIS information. Future developments like irreversible CTMC processes may offer even more biological realism.

Methods

For many spatial scales and problems, geography can naturally be partitioned into a finite number of discrete sites for . Examples of these situations include individual cities, islands or countries. Starting from the observed data, at the tips of the phylogeny we record discretized locations , where pin-points the sampling site of taxon . Unobserved in the spatial process are the locations of the most common recent ancestor drawing from root distribution , the times at which the descendent taxa move and amongst which discrete sites these moves occur, a process which ultimately gives rise to . Conditioning on and the unobserved locations realized at each internal node [13],[16], suggest modeling the instantaneous locations for taxa along each branch in as independent continuous-time Markov chains (CTMCs). CTMCs are the same processes one commonly exploits to model sequence character evolution [15],[50]. Although many readers are familiar with CTMCs, we here highlight several chain properties to which we turn later when discussing CTMC modeling limitations. CTMCs are the simplest stochastic processes that emit discrete outcomes as a continuous function of time. The processes are memoryless, in that the probability of transitioning to a new location only depends on the current location and not the past history. A infinitesimal rate matrix completely characterizes the CTMC process. Rate matrix contains non-negative off-diagonal entries and all rows sum to 0, yielding a stochastic matrix upon exponentiation. Solving the Chapman-Kolmogorov equation that specifies the behavior of the chain yields the finite-time transition probabilities . In matrix form,

(1)

Determining the finite-time transition probabilities involves matrix exponentiation, generally accomplished through an eigen-decomposition of . Here, we restrict ourselves to infinitesimal rate matrices that yield only real eigen-values and eigen-vectors. Any matrix similar to a symmetric matrix ensures a real eigen-decomposition; consequentially, we formulate

(2)

where is an overall rate scalar, is a symmetric matrix and . Infinitesimal rate matrices of this form generate reversible Markov chains, such that

(3)

placing many restrictions on the underlying geographic process. We discuss these limitations and modeling extensions that allow for irreversible chains in the Discussion. In its most general time-reversibile (GTR) form, contains free parameters, with donating together with 's off-diagonal entries. Following standard practice, we normalize entries in such that measures the expected (with respect to ) number of transitions per unit time .

One illuminating perspective from which to view a CTMC is that of a random walk on a graph . From this perspective, the possible realizations of the chain correspond to the vertex set of . Between the vertices lie edges that record the infinitesimal transition rates. For example, between and sits . As a continuous-time random walk, a particle, starting at vertex at time , first waits an Exponential amount of time with rate and then randomly decides to which neighboring vertix to move with probability . Now on , the process repeats. Neighboring vertices are those for which a single edge connects them. For character evolution, “complete” graphs find almost exclusive use, such that edges exist between all pairs of vertices. At a minimum, however, the graph must remain “connected”, such that it remains possible to walk between any two vertices on .

Bayesian stochastic search variable selection

When GTR models find use modeling nucleotide substitution, most of the possible transitions have non-neglible probability of occurring and are observed over the evolutionary history. Such is unlikely to be the case for geographical locations; given that there may be many sites and each taxon only has one location (the equivalent of just one single alignment site), we expect most transitions to rarely, if ever, occur. Consequentially, we suspect a priori that many infinitesimal rates are zero. From a statistical perspective, so many degrees of freedom fit to the limited data lead to extremely high variance estimates. These poor estimates arise not only for , but, more critically, for inference of the unobserved ancestral locations and . We circumvent this sparse data problem by invoking BSSVS to select a parsimonious parameterization of . BSSVS enables us to simultaneously determine which infinitesimal rates are zero depending on the evidence in the data and efficiently infer the ancestral locations. As a beneficial by-product of BSSVS, directly quantifying the evidence about which rates are non-zero furnishes both the most likely migration patterns and the ability to test between competing migratory hypotheses.

BSSVS is traditionally applied to model selection problems in a linear regression framework, in which statisticians start with a large number of potential predictors and ask which among these associate linearly with an -dimensional outcome variable . For example, the full model becomes , where is a -dimensional vector of regression coefficients and is an -dimensional vector of normally distributed errors with mean . When for differs significantly from 0, helps predict , otherwise contributes little additional information and warrants removal from the model via forcing . Given potentially high correlation between the predictors, deterministic model search strategies tend not to find the optimal set of predictors unless one explores all possible subsets. This exploration is generally computationally impractical as there exist such subsets and completely fails for .

Recent work in BSSVS [51],[52] efficiently performs the exploration in two steps. In the first step, the approach augments the model state-space with a set of binary indicator variables and imposes a prior on the regression coefficients that has expectation and variance proportional to a diagonal matrix with its entries equal to . If , then the prior variance on shrinks to 0 and enforces in the posterior. In the second step, MCMC explores the joint space of simultaneously.

To map BSSVS into the phylogeography setting, we consider selection among the random graphs in which each of the edges either exists or does not exist in . Let be the binary indicator that an edge exists connecting and . An equivalent parameterization specifies that when and otherwise. So, rate matrix plays an analogous role to the regression coefficients in BSSVS. An important difference is that while , mandating alternative prior formulations.

Prior specification

To specify a prior distribution over , we assume that each indicator acts a priori as an independent Bernoulli random variable (RV) with small success probability . The sum of independent Bernoulli RVs yields a Binomial distribution over their sum . In the limit that , this prior conveniently collapses to

(4)

where is the prior expected number of edges in graph .

We entertain two prior choices for . Diagonal vector is the stationary distribution for the CTMC when all edges are included in the graph . For this complete graph, as the length of the random walk , one expects that with probability . One natural choice says that there exists no spatial preference over the long-run and fixes for all . However, sites may expound spatial preference over the long-run; for example, such preference can relate to known or unknown quantities such as population-size or geographic size of the site. In these situations, we estimate simultaneously with the rest of the model by imposing the flat prior . Also non-informatively for small values, we take .

To complete the CTMC specification, we assume that all unnormalized rates in are a priori independent and Gamma-distributed with prior expectation and variance , following in the vein of Bayesian SSVS. However, little previous work on prior formulation helps inform our choices of and . This represents a critically important area of research. A common, yet arbitrary choice in the Bayesian phylogenetic literature assumes that rates draw from Exponential distributions, forcing . We follow this practice in light of there being no obvious way to elicit information on the variance of these rates. Finally, we explore two choices for setting the means. The first assumes no preference over rates, setting all , where is an arbitrary constant; as, after normalization, only ratios of infinitesimal rates participate in the data likelihood, the actual value of has no influence on the likelihood. The second is informed by the geographical distance between sites.

Distance informed prior

Considerable additional information exists about the sites and remains unused. Most notably, the geographic distances between (the centroids) of sites is readily measurable. A priori we may believe that more distantly separated sites have the smallest infinitesimal migration rates, yielding

(5)

Other information is also surely helpful and application-specific. One example involving human hosts quantifies the availability of motorized transportation, such as air flights, between sites. We explore the utility BSSVS and distance-informed priors in our phylogeographic models.

Bayes factor test of significant diffusion rates

The Bayes factor (BF) for a particular rate contributing to the migration graph is the posterior odds that rate is non-zero divided by the equivalent prior odds,

(6)

where is the posterior probability that rate is non-zero, in this case the posterior expectation of indicator . Since we employ a truncated Poisson prior with mean , that assigns 50% prior probability on the minimal rate configuration (), the prior probability reduces to

(7)

We consider rates yielding a BF3 as well supported diffusion rates constituting the migration graph.

Sequence data, joint modeling and posterior sampling

A strength of the Bayesian approach we exploit in this paper is the ability to integrate together into a joint model of spatial locations and aligned molecular sequence data collected from the taxa. The joint model affords a natural way to incorporate uncertainty about the unobserved phylogeny and the character substitution process giving rise to . We take a standard statistical phylogenetic approach and assume that a separate CTMC characterized by generates . While we discuss specific choices about this process in the Results sections, we do assume that the sequence and location CTMCs are independent given , enabling us to write the joint model posterior distribution as

(8)

Likelihoods and follow directly from Felsenstein's pruning algorithm [15], efficiently integrating over all possible locations and sequences at the root and internal nodes in .

We approximate the joint posterior (8) and its marginalizations using MCMC implemented in the software package BEAST [43]. We employ standard transition kernels over the parameter spaces of and . To sample realizations of , we consider random-walk operators on the continuous portions and a specialized “bit-flip” operator on the Bernoulli rate indicators . [53] discuss this transition kernel further. Finally, in many situations, inference on the posterior distribution of the root and internal node states is of paramount interest. We implement a pre-order, tree-traversal algorithm in BEAST that allows researchers to generate realizations of the root and internal node states following [20] and produce posterior summaries. Importantly, this procedure is not limited to phylogeographic models, making general ancestral state reconstruction available in BEAST for the first time.

Summarizing posterior location uncertainty

An important statistical question asks to what extent the data inform our inference when fitting different phylogeographic models. A model of low statistical power makes poor use of the information in the data, while a successful model exploits this information to generate posterior distributions that are maximally different from prior beliefs. One primary outcome of a Bayesian phylogeographic study is the marginal posterior distribution of the root location . We calculate the Kullback-Leibler (KL) divergence [54] from the root location prior to summarize this information gain,

(9)

where . When the posterior and prior distributions are equal, . In the examples in this paper, we fix and achieves its maximum when the posterior places all estimable mass on a single location. From this perspective, plays the role of a measure of dispersion [55] or uncertainty. As a simple numerical summary, we also use to explore the utility of BSSVS and distance-informed priors on drawing inference from phylogeographic models. Larger values signify that the model extracts more information from the data. To calculate KL divergence, we employ a uniform discrete distribution as reference distribution.

Association index

Following existing phylogeographic approaches, we finally score the degree of spatial admixture using a modified association index (AI) [35]. For a given phylogeny and tip locations , we obtain the association value by summing over each internal node ,

(10)

where counts the number of sampled locations descendent to and is the frequency of the highest frequency location amongst these descendents. Similar to [10], we report the posterior distributions and the AI compares these distributions to those obtained by random permutation of the tip locations . Deviation from these permuted distributions reflected in low AI values suggests phylogeographic structure whereas AI values close to 1 suggest spatial admixture.

Visualizing phylogeographic diffusion

To summarize the posterior distribution of ancestral location states, we annotate nodes in the MCC tree with the modal location state for each node using TreeAnnotator, and visualize this tree using FigTree (available at http://tree.bio.ed.ac.uk/software). To provide a spatial projection, we convert the tree into a keyhole markup language (KML) file suitable for viewing with Google Earth (http://earth.google.com). We introduce the temporal information on the marked-up tree using the KML-function to animate viral dispersal over the time. Example KML files for the Avian Influenza A HA and NA genes are included as supplementary files and software to convert annotated trees to KML is available from the authors on request.

Supporting Information