BACKGROUND
evolutionary and phylogenetic analysis has become an important tool in the study of established and emerging viruses, including influenza e.g.
<cit> , hiv e.g.
<cit> , mers e.g.
<cit> , and most recently, ebola  <cit> .
the application of evolutionary concepts to epidemiological surveillance and outbreak control has been transformed in recent years by the increasing availability of viral genome sequences, the growth in computer processing power, and the development of sophisticated analytical methods .
information gleaned from genetic data has the potential to identify factors that influence the spread and evolution of viral diseases.
an evolutionary approach to epidemiology provides several benefits when combined with traditional methods, as it allows reconstruction of epidemic transmission history from a small number of pathogens sampled shortly after the discovery of an outbreak, the estimation of migration histories, and the inference of transmission links among cases that may not be evident using data on spatio-temporal incidence alone.
evolutionary approaches are particularly powerful in the context of rna viruses because these viruses are characterised by very rapid evolution  <cit> , such that their evolutionary and ecological processes occur on the same time scale  <cit> .
as a consequence the pattern of genetic differences among viruses sampled from a population contains information about transmission history.
specifically, phylogeographic and statistical studies of virus genomes sampled through space and time can answer questions relating to the geographic dissemination of epidemics e.g.
<cit>  and numerous methods are available to infer patterns of spatial spread from genetic data, e.g.
.
however, the explanatory power of viral phylogeography would be greater if it could identify and test the effects of specific aspects of the environment on spatial spread and growth.
this question represents a major focus of current research in viral phylogeography.
lemey et al.
<cit>  developed a method based on relaxed random walk models that can reconstruct virus dispersal in continuous space using viral genomes sampled from known locations at known times.
this method co-estimates both the virus phylogeny and the locations of unsampled common ancestors, thereby producing a full spatial dispersal history of the sampled infections.
this history includes an estimate of the location of epidemic origin and of the velocity, direction and heterogeneity of spatial spread.
the method has been used in several studies to reconstruct the spatial spread of pathogenic viruses, including crop diseases, dengue virus, west nile virus and avian influenza viruses .
however, while this method can place phylogenies in a geographical context, it does not explicitly incorporate environmental differences across the geographic landscape within which transmission occurs.
such heterogeneity may have a significant effect on spatial dissemination, especially for outbreaks in natural animal and plant populations e.g.
.
at present there is only one phylogeographic method that explicitly tests potential predictors of spatial spread.
this is the “phylogeographic glm” approach  <cit>  that estimates rates of lineage movement among a fixed number of discrete locations, and in doing so parameterises each among-location rate as a linear function of one or more predictor variables.
the coefficients of the linear model are then co-estimated with the among-location rates, the phylogeny, and other parameters, using bayesian monte carlo markov chain  inference.
however, this framework is only currently applicable to discretised spatial locations and is more suited to hypotheses concerning human mobility among population centres via transportation hubs, such as airports  <cit>  and is less suitable for hypotheses concerning the dissemination of animal and plant pathogens throughout natural landscapes.
further, it is computationally demanding because predictor variables are assessed whilst sampling phylogenies using a bayesian mcmc algorithm.
there is therefore a need for complementary approaches.
here we present and apply a new analytical framework that aims to integrate landscape ecology with a phylogenetic approach to molecular epidemiology.
this framework uses spatial information obtained from phylogeography to study the impact of environmental variables on epidemic dispersal in continuous space.
our method differs from previous approaches in that it explicitly separates the task of testing correlates of lineage dispersal from the task of inferring the history of movement from genetic data.
this separation has two benefits.
first, it allows large numbers of environmental variables to be assessed without the need for phylogenetic mcmc sampling.
second, it increases flexibility, because it can be applied to any dispersal history, irrespective of the specific phylogeographic method or software package that was used to infer that history.
we illustrate our framework by applying it to a well-characterised viral outbreak, the spread of rabies in raccoons across the north-east of the usa over approximately thirty years .
this data set has been analysed using previous phylogeographic methods  <cit> , thus enabling direct comparison with the results obtained here.
methods
overview of methodology
the structure of our framework can be summarised in the following five steps.
a full description of each step is provided in the next section.
in the first step, the history of lineage dispersal is recovered from one or more spatial- and temporally-referenced phylogenies .
such trees are generated by the continuous phylogeography method implemented in beast  <cit> , but it is important to note that our method is applicable phylogenies from any source that are annotated with dates and locations in the same way .
the velocity, distance and duration of spatial movement along each branch in each tree are extracted and represented by a vector.fig.
<dig> an illustration of the node position randomisation procedure used to generate null distributions of the d statistic.
a the original environmental raster  upon which is superimposed the movement events extracted from one spatiotemporally-referenced phylogeny.
b the result of one randomisation of node positions.
this randomisation procedure is performed within a minimum convex hull , which is defined by the node locations of all selected phylogenies
in a second, optional, step, the spatio-temporal data extracted in step  is used to calculate summary statistics of spatial spread, such as the velocity of the epidemic wavefront.
each of the vectors obtained in step  is assigned a “weight” score, which is calculated using a raster that defines the spatial heterogeneity of a specified environmental variable .
we have implemented three different models of spatial movement  to compute the “weight” allocated to each phylogeny branch for a given environmental raster.
the “straight-line path” model, which assumes that movement is in a straight-line between the start and end locations of each branch.
the “least-cost path” model, which uses a least-cost algorithm  <cit>  to determine the route taken between the start and end points.
the “random walk path” model, which uses circuit theory  <cit>  to accommodate uncertainty in the route taken.
note that for models  and , each environmental raster must be considered twice, once as a conductance factor  and once as a resistance factor .
the correlation between the duration of each phylogeny branch and its “weight”  is then estimated.
correlations are repeated for each environmental raster, and for each of the three path models described above.
finally, the statistical significance of these correlations is tested using a null model.
to generate this null model, we implement a simple randomisation procedure: phylogenetic node positions are randomised within the study area, under the constraint that branch lengths, tree topology and root position are unchanged .
description of methodology
step 1: extracting spatio-temporal information from phylogenies
the input data for our method consists of one or more spatiotemporally-referenced phylogenies .
if multiple phylogenies are used then each is assumed to be equally probable.
at present, such trees are mostly likely to be generated using the phylogeographic models implemented in beast  <cit> , but our framework can be applied to phylogenies from any source that have been annotated in the same way.
in order to extract the spatio-temporal information in the input trees, each phylogeny branch in each tree is summarised as a vector defined by its start and end location, and its start and end dates.
each branch therefore represents an independent viral lineage dispersal event  <cit> .
vectors for each input tree are tabulated.
step 2: estimation of dispersal and epidemiological statistics
as an optional step, our analytical framework calculates several spatio-temporal statistics from these dispersal vectors.
so far, we have implemented three such statistics: the mean lineage dispersal velocity, the mean lineage diffusion coefficient, and a time-series of the maximal epidemic wavefront distance.
these statistics are visually summarised as plots of  the kernel density of lineage dispersal velocity parameters ,  the kernel density of lineage diffusion coefficient parameters , and  the change through time in the spatial and patristic maximal wavefront distances, as measured from the location of the tree root.
the spatial distance corresponds to a straight-line distance, whereas the patristic distance equals the sum of the spatial distances along each phylogeny branch between the tree root and its tips.
step 3: computation of environmental weights
in this step, each of the vectors obtained in step  <dig>  is assigned a specific “weight” based on a raster of values that represent the type or magnitude of an environmental variable .
in order to compute these “weights”, we have implemented three different “path models” which represent the path taken by a phylogeny branch as it travels between its start and end locations: straight-line path: each lineage travels in a straight line between its start and end location.
in this case, the branch “weight” is computed as the sum of the values of the raster cells through which the straight line passes.
least-cost path: each lineage travels via the least-cost path between its start and end location  <cit> .
with this model, the branch “weight” is computed as the sum of the transition values between adjacent cells along the least-cost path  <cit> .
random walk path: each lineage travels via a random walk between its start and end location.
in this case, the “weight” is a graph-theoretic metric based on circuit theory, which takes into account multiple possible pathways connecting a given pair of locations, and the values of the raster cells through which they pass  <cit> .
specifically, for each lineage, one node is connected to a one-ampere current source, while the other is connected to ground.
the “weight” between each pair of nodes is then defined by the effective electric resistance or conductance connecting them on the grid, i.e.
the environmental raster in question  <cit> .
more details about the circuit theory model can be found in  <cit>  and  <cit> .
for the latter two models, an environmental variable can be treated either as a conductance or a resistance factor.
for example, if an environmental raster is treated as a resistance factor then raster cells with low values will be more permeable to dispersal and those with high values will correspond to poor dispersal habitat or to movement barriers  <cit> .
when there is no prior information about whether a given variable will facilitate or impede lineage movement, it is prudent to consider it twice, once as a potential conductance factor, and once as a potential resistance factor.
by default, we use the environmental values as resistance or conductance values.
log-transformed and/or standardised values may be necessary if the current univariate approach is in future extended to a multivariate one .
note, when all cells in a raster are identically-valued  then the total “weight” for any path is simply proportional to its geographic distance.
a special case exists when all cells have values equal to  <dig>  which we here define as the “null raster”.
step 4: correlation analyses
we next calculate the regression between the duration of each movement event  and the “weight” computed for that branch.
a separate regression is performed for each environmental factor of interest.
the absolute strength of the regression can be measured using its coefficient of determination .
specifically, we use the statistic d = , where r2env is the r <dig> obtained when branch durations are regressed against weights defined by the environmental raster, and r2null is the r <dig> obtained when branch durations are regressed against weights defined by the null raster .
all environmental raster cell values are increased by  <dig>  to enable a direct comparison of r2env and r2null.
d is therefore a correlation measure relative to a null hypothesis and represents the degree to which the regression is strengthened when spatial variation in the environmental variable is taken into account.
step 5: significance testing using randomisation
in the final step, the statistical significance of d is tested.
in order to calculate a null distribution for the d statistic, we have implemented a randomisation procedure that randomises phylogenetic node positions under the constraint that branch lengths , the tree topology and the root position are unchanged .
furthermore, this randomisation procedure is not applied to the entire raster but to a subset of it defined by the minimum convex hull around a set of locations.
the set of locations comprises the start and end positions of each branch in the spatiotemporally-referenced input phylogenies.
during the randomisation, if a randomised branch position falls outside the area defined by this convex hull, the algorithm randomises its position again until it does not fall outside this area.
by constraining permutation within this convex hull we are, in effect, using the data to inform the relevant study area.
in contrast, the square raster inputted by the user may be of arbitrary size and orientation.
without this constraint, randomised phylogenies will fall in areas with different environmental values, potentially leading to type i errors.
the result of this procedure is one p-value per input phylogeny, which equals the proportion of randomisation replicates that generated d values larger than that generated by the empirical input phylogeny.
if multiple input trees are used then a distribution of p-values is obtained and we report the percentage of p-values >  <dig> .
it is worth noting that the statistical significance of our test is derived from a randomisation distribution, not from a distribution derived from the regression model.
hence the coefficients of determination r2env and r2null are treated only as data statistics, and therefore our test is not dependent on regression model assumptions .
application to rabies in north america
the raccoon rabies virus data set comprises  <dig> sequences with known sampling dates and locations.
the sequences are ~2800 nt long and span the viral nucleoprotein n gene and the 5' end of the phosphoprotein p gene; see biek et al.
<cit>  for further details).
we extracted the spatial-temporal information contained in  <dig> spatiotemporally-referenced trees sampled from the posterior distribution of trees inferred for this data set by lemey et al.
<cit> .
we first calculated three epidemiological statistics from these  <dig> input trees: mean lineage dispersal velocity, mean lineage diffusion coefficient and a plot of the maximal wavefront distance.
for ease of explanation, we initially show how our analytical framework works on a single input phylogeny.
subsequently we show how the method is extended to incorporate phylogenetic uncertainty by analysing multiple input trees.
we investigated six environmental variables to determine if they were associated with the dispersal rate of raccoon rabies lineages.
these included the three most important igbp  land cover variables for the study area , as well as elevation, human population density and “inaccessibility” .
the “forests” and “savannas” layers combine several igbp land cover layers .
the environmental rasters are shown in fig.
<dig>  the sources of the original raster files are given in table  <dig>  the original data presented a resolution of  <dig>  arcmin, corresponding to cells ~1 km square.
we generated distinct land cover rasters from the original data by creating lower resolution rasters  whose cell values equalled the number of occurrences of each land cover category within the  <dig> arcmin cells.
the resolution of the three other original rasters was also decreased to  <dig> arcmin for tractability.
for each selected raster and path model combination, we performed  <dig> randomisations using the randomisation procedure outlined in step  <dig> above.fig.
<dig> the six environmental variables that were tested in the analysis of the raccoon rabies virus data set.
the region shown corresponds to the northeast of the usa, centered approximately on harrisburg, pa. details of the construction and source data for these rasters is provided in the main text
www.igbp.net
webmap.ornl.gov
www.map.ox.ac.uk
bioval.jr.ec.europea.eu
performance on simulated data sets
we used artificial data sets to measure the potential type i error rate of the randomisation procedure implemented to assess the significance of observed d values.
artificial data sets were generated using four different approaches:  the duration of each movement event was randomly permuted among branches,  the duration of each movement event was equal to the spatial distance along the path between its start and end locations,  the duration of each movement event was equal to the environmental “weight” assigned to the path between the locations, and  the duration is equal to a value randomly and uniformly chosen between the values defined in cases  and .
case  thus corresponds to case  but with the addition of random noise.
simulation approaches  and  generate data sets under the null hypothesis , whereas approaches  and  generate data sets under the alternate hypothesis .
we generated  <dig> artificial data sets for each of these three simulation approaches.
the artificial data sets used the empirical phylogenies from the raccoon rabies data set , but replaced the empirical branch durations with new durations as specified above.
for approach , environmental “weights” were calculated using the “elevation” raster treated as a resistance factor.
this raster was chosen because elevation appears to impede the movement of raccoon rabies virus lineages .
for each path model and simulation method, we report the percentage of artificial data sets for which the randomisation null hypothesis test p-value was <  <dig> .
for simulation approaches  and , these percentages correspond to estimates of the type i error rate of the test.
for approaches  and , these percentages indicate statistical power.
RESULTS
we use the methods outlined above to analyse previously published virus gene sequences sampled from an epidemic of rabies virus in raccoons in north america  <cit> .
spatio-temporal phylogenies were estimated from these sequences using the continuous phylogeographic model implemented in beast  <cit> .
the analyses below were performed on  <dig> phylogenies sampled regularly from the post burn-in posterior tree distribution generated by beast.
estimation of spatio-temporal statistics
we used our framework to calculate several statistics that summarise the spatio-temporal information inherent in the rabies virus phylogenies .
these statistics include the epidemic wavefront and patristic distances through time , and parameters that describe dispersal velocities and diffusion coefficients .
the epidemic wavefront plot  shows that the spatial extent of the raccoon rabies epidemic increased at a relatively constant rate until around  <dig>  the patristic distance  at this time is greater  than the spatial distance .
the difference between these two measures is to be expected, because the patristic distance includes movements in all directions, including those back towards the epidemic origin .
hence dispersal was diffusive in nature and did not always follow the shortest or most direct path from epidemic origin to epidemic wavefront.
kernel density graphs of the mean lineage velocity  and phylogenetic diffusion coefficient  indicate significant variation in these parameters among lineages.
the estimated mean lineage velocity was ~37 km/year , similar to that reported for dengue virus in vietnam  but substantially smaller than that of the west nile virus epidemic in north america .
these statistics are helpful in quantifying the spatial dynamics of an epidemic and can also be used to compare different outbreaks.
some, but not all, of them are implemented in the software package spread  <cit> .
we hope that future work will lead to the development and implementation of further summary statistics.fig.
<dig> epidemiological statistics estimated from the raccoon rabies virus data set.
a time-series of the spatial distance between epidemic origin and maximal epidemic wavefront, and b evolution of the patristic distance between epidemic origin and maximal epidemic wavefront, c kernel density estimates of lineage velocity parameters and d kernel density estimates of lineage diffusion coefficient parameters .
in parts a and b the grey area corresponds to the 95 % credible region of the estimated wavefront position.
in parts c and d the three contours show, in shades of decreasing darkness, the 25 %, 50 %, and 75 % highest posterior density regions via kernel density estimation
impact of environmental variables
next, we used our framework to measure and test the correlation of various environmental variables with the movement of rabies virus lineages .
we explored rasters that represent six different environmental variables, specifically  cropland land cover,  forest land cover,  savanna land cover,  elevation,  human population density and  inaccessibility to major cities.
for illustrative purposes, we first present the results based on only one tree .
we later show how the framework can be extended to incorporate phylogenetic uncertainty, by analysing multiple input trees.
the purpose of the analysis is to calculate a test statistic, d, for each environmental variable, e. the statistic d represents the degree to which the correlation between the duration and “weight” of each lineage is strengthened when spatial heterogeneity in e is taken into account .
formally, d = r2env - r2null, where r2env is the coefficient of determination obtained when branch durations are regressed against branch weights calculated using raster e, and r2null is the coefficient of determination obtained when branch durations are regressed against branch weights defined by a null raster with no spatial variation .
if d ≤  <dig>  then the environmental variable e does not explain variation in branch duration better than geographic distance alone.
however, if d is strongly positive, then the values of e are strongly associated with branches that move more rapidly or more slowly than average.
figure  <dig> displays the linear regressions between branch durations and branch weights for a rabies virus phylogeny randomly chosen from the posterior tree distribution.
regressions are shown for “weights” calculated using the null raster  and for “weights” calculated using the elevation raster, when elevation is treated as a factor that impedes movement .
results are shown for each of the three path models .
in this example, there is comparatively little difference among path models in the strength of the relationship between branch durations and branch weights.
for all three path models, r2env is approximately twice as great as r2null, so the corresponding d values are positive .
thus while geographic distance alone explains a small part of among-lineage variation in velocity , this variation is better explained when the geographic elevation of the path taken by each lineage is taken into account .
table  <dig> reports d values for other environmental variables when applied to the same single phylogeny as that used in fig.
<dig> fig.
<dig> linear regressions between branch durations and branch weights for one phylogenetic tree.
plots in  were generated by calculating the branch weights using the “null” raster, whereas plots in  were generated by calculating the branch weights using the “elevation” raster treated as a resistance factor.
plots are shown for each of the three path models .
the d value obtained under each path model is shown at the top.
the plots show that the r <dig> of the regression is approximately doubled when spatial heterogeneity in elevation is taken into account
for each combination of environmental variable and path model, the test is applied to a single phylogeny, and to a set of  <dig> trees.
for the former, the d statistic and p-value of the test are shown.
for the latter, we report the percentage of trees for which p <  <dig>
p-value <  <dig>
in each case, artificial data sets were applied to the “elevation” raster treated as a resistance factor.
values equal the percentage of  <dig> artificial data sets for which the p-value of the null hypothesis test was <  <dig>
the statistical significance of each d value can be assessed by comparing it to a null distribution generated by randomising the empirical node positions .
in each case,  <dig> randomisation replicates were performed.
table  <dig> reports the p-values obtained when this randomisation procedure is applied to the rabies virus phylogeny introduced above.
for this data set only some of the d values are statistically significant.
specifically, significant positive d values were obtained for the “elevation” and “inaccessibility” rasters  under the least-cost path model .
although these results are interesting they are based on only one tree and therefore do not take in account the statistical uncertainty arising from phylogenetic inference.
we will now demonstrate how the null hypothesis test is performed on  <dig> phylogenies sampled from the posterior distribution of trees.
first, we estimated d values for each of the  <dig> trees.
this generates, for each combination of environmental variable and path model, a distribution of  <dig> d values.
the randomisation procedure is then applied to each d value, precisely as outlined above.
this results in  <dig> p-values for each raster/path model combination.
we therefore report the percentage of trees which give rise to p-values <  <dig>  ; this can be interpreted as the posterior probability of observing a significant correlation between lineage movements and the environmental variable.
we illustrate this procedure in fig.
<dig>  which compares the empirical distribution of d  with five replicates of the null distribution of d generated by the randomisation procedure .
in fig.
5a the distributions were calculated using the “elevation” raster  and in fig.
5b they were calculated using the “forests” raster .
in both cases the least-cost path model was used.
in fig.
5a, the empirical distribution is greater than the randomisation replicates, indicating that the environmental variable in question explains viral lineage movement, whereas in fig.
5b there is no noticeable difference between the empirical and null d values.fig.
<dig> empirical distributions of the d statistic , calculated from  <dig> trees sampled using bayesian mcmc inference.
these are compared with five replicates of the null distribution of d generated by the randomisation procedure .
in  the distributions were calculated using the “elevation” raster  and in  they were calculated using the “forests” raster .
in both cases the least-cost path model was used.
for visual clarity, discrete histograms were converted into density curves using a gaussian smoothing kernel
the environmental variables in table  <dig> can be grouped into three categories:  rasters that have little or no association with lineage movement ,  rasters for which there is weak evidence that they affect lineage movement , and  rasters for which there relatively stronger evidence that the environmental variable is associated with slower or faster movements .
the “forests” and “savannas” rasters fall within the first category, whilst “croplands”  and “human population density”  belong to the second category.
only the “elevation” raster  clearly belongs to the third category.
the importance of the “inaccessibility” raster  is less clear and is sensitive to the path model chosen.
for both the “elevation” and “inaccessibility” rasters, the fraction of trees with p <  <dig>  is smaller under the straight-line path model than under the other two path models.
this may reflect the over-simplified nature of the straight-line model, which permits biologically unrealistic scenarios, such as paths that traverse large water bodies.
performance on simulated data sets
table  <dig> shows the results of null hypothesis tests performed on artificial data sets.
these data sets were simulated using three different approaches.
simulation approaches  and  generate artificial data sets under the null hypothesis .
therefore for these two approaches the fraction of significant results should equal the critical value of the test.
this is indeed the case, as only 0–7 % of tests on artificial data produce a p-value <  <dig>  and hence the type i error of the test appears to be appropriate.
simulation approaches  and  generates data sets under the alternate hypothesis, i.e.
branch duration is equal or influenced by the environmental raster.
the null hypothesis was almost always rejected  when the test was applied to data generated using approaches  or , indicating that the test has reasonable statistical power in the context of this data set.
note that very similar results were also obtained on data sets whose simulations were based on alternative environmental rasters .
CONCLUSIONS
in this study we show that it is possible to integrate genomic and environmental data in order to test hypotheses concerning the mode and tempo of virus dispersal during emerging epidemics.
the raccoon rabies data set explored here was chosen for illustrative purposes and our results strongly support the notion that increasing elevation is associated with slower movements of this virus.
this result is biologically plausible and the impact of elevation on the dissemination of raccoon rabies was previously addressed by biek et al.
<cit> .
their analysis used spatial kriging to compute annual contours from the date of the first reported case of each county and the contours were coloured according to temporal periods identified by phylogenetic molecular clock analysis.
they then visually compared the resulting contours overlaid on an elevation map, enabling them to posit that mountain ranges likely formed a barrier to raccoon movement.
we intend to explore the application of our new method to other viral outbreaks in future work.
we have implemented three path models in our analytical framework.
although the straight-line path model is perhaps overly simplistic, the appropriateness of one model over another will depend on the locomotive behaviour of the pathogen’s host species .
in the case of the raccoon rabies epidemic, we had no prior assumption about which path model might be most appropriate and therefore we decided to test them all.
in other instances, there may be good reasons to choose a priori one path model over another.
for instance, one might chose the least-cost path model over random walk path model in some cases, for example viruses carried by migratory birds; i.e.
the least-cost model may be more appropriate when host-species intentionally reach specific locations by avoiding non-suitable landscape areas.
simulations indicate that the randomisation procedure used to test the null hypothesis  has appropriate type i error rates and acceptable statistical power .
however, it is important to note that a significant result represents evidence only for a correlation between virus dispersal and an aspect of the geographic landscape.
in future work we intend to explore whether generalized linear models can be applied to this framework in order to account for possible correlations among the environmental variables in question.
in the meantime, the broader ecology and epidemiology of the pathogen, host, and habitat concerned should be taken into account when interpreting results.
currently, our analytical framework does not allow for the analysis of environmental variables, such as temperature or humidity, that vary notably over the course of the epidemic under investigation.
the assumption of constant values is unlikely to be unrealistic for the environmental variables investigated here, although some variation has undoubtedly occurred in human population density.
the extension of our approach to time-varying environments is theoretically feasible, but poses significant technical and practical problems, specifically  the acquisition of a series of environmental rasters that represents change in the variable concerned at a sufficiently precise temporal resolution, and  generalisation of the least-cost and random-walk path models to three rather than two dimensions, with the third dimension corresponding to rasters that represent different points in time.
a second improvement for future work would be to allow the co-analysis of several environmental factors in a multivariate framework.
such a framework would require correlation statistics, for example those based on a generalised linear model , as in  <cit> .
implementation
the methods introduced here are implemented in r  and the scripts are freely available from http://evolve.zoo.ox.ac.uk/evolve/software.html.
example files and a tutorial are available as additional file  <dig>
additional file
additional file 1:
example files and tutorial related to the present study.
competing interests
the authors declare that they have no competing interests.
authors’ contributions
sd, rr and ogp designed the study and the methodology.
sd developed the analytical framework and performed the study.
sd and ogp wrote the manuscript.
all authors read, edited and approved the final manuscript.
we thank two anonymous reviewers, roman biek, kenneth hoehn, nuno faria, philippe lemey and stéphane dray for their useful comments, as well as dan weiss, peter gething, clément tisseuil and marius gilbert for their advice regarding the processing and management of gis data.
this work was funded by the wiener-anspach foundation.
ogp received funding from the european research council under the european union's seventh framework programme  / erc grant agreement no.
614725-pathphylodyn.

 evolutionary and phylogenetic analysis has become an important tool in the study of established and emerging viruses, including influenza e.g., hiv e.g., mers e.g., and most recently, ebola. the application of evolutionary concepts to epidemiological surveillance and outbreak control has been transformed in recent years by the increasing availability of viral genome sequences, the growth in computer processing power, and the development of sophisticated analytical methods. information gleaned from genetic data has the potential to identify factors that influence the spread and evolution of viral diseases. an evolutionary approach to epidemiology provides several benefits when combined with traditional methods, as it allows reconstruction of epidemic transmission history from a small number of pathogens sampled shortly after the discovery of an outbreak, the estimation of migration histories, and the inference of transmission links among cases that may not be evident using data on spatio-temporal incidence alone. evolutionary approaches are particularly powerful in the context of rna viruses because these viruses are characterised by very rapid evolution, such that their evolutionary and ecological processes occur on the same time scale. as a consequence the pattern of genetic differences among viruses sampled from a population contains information about transmission history. specifically, phylogeographic and statistical studies of virus genomes sampled through space and time can answer questions relating to the geographic dissemination of epidemics e.g. and numerous methods are available to infer patterns of spatial spread from genetic data, e.g.. however, the explanatory power of viral phylogeography would be greater if it could identify and test the effects of specific aspects of the environment on spatial spread and growth. this question represents a major focus of current research in viral phylogeography. lemey et al. developed a method based on relaxed random walk models that can reconstruct virus dispersal in continuous space using viral genomes sampled from known locations at known times. this method co-estimates both the virus phylogeny and the locations of unsampled common ancestors, thereby producing a full spatial dispersal history of the sampled infections. this history includes an estimate of the location of epidemic origin and of the velocity, direction and heterogeneity of spatial spread. the method has been used in several studies to reconstruct the spatial spread of pathogenic viruses, including crop diseases, dengue virus, west nile virus and avian influenza viruses. however, while this method can place phylogenies in a geographical context, it does not explicitly incorporate environmental differences across the geographic landscape within which transmission occurs. such heterogeneity may have a significant effect on spatial dissemination, especially for outbreaks in natural animal and plant populations e.g.. at present there is only one phylogeographic method that explicitly tests potential predictors of spatial spread. this is the “phylogeographic glm” approach that estimates rates of lineage movement among a fixed number of discrete locations, and in doing so parameterises each among-location rate as a linear function of one or more predictor variables. the coefficients of the linear model are then co-estimated with the among-location rates, the phylogeny, and other parameters, using bayesian monte carlo markov chain inference. however, this framework is only currently applicable to discretised spatial locations and is more suited to hypotheses concerning human mobility among population centres via transportation hubs, such as airports and is less suitable for hypotheses concerning the dissemination of animal and plant pathogens throughout natural landscapes. further, it is computationally demanding because predictor variables are assessed whilst sampling phylogenies using a bayesian mcmc algorithm. there is therefore a need for complementary approaches. here we present and apply a new analytical framework that aims to integrate landscape ecology with a phylogenetic approach to molecular epidemiology. this framework uses spatial information obtained from phylogeography to study the impact of environmental variables on epidemic dispersal in continuous space. our method differs from previous approaches in that it explicitly separates the task of testing correlates of lineage dispersal from the task of inferring the history of movement from genetic data. this separation has two benefits. first, it allows large numbers of environmental variables to be assessed without the need for phylogenetic mcmc sampling. second, it increases flexibility, because it can be applied to any dispersal history, irrespective of the specific phylogeographic method or software package that was used to infer that history. we illustrate our framework by applying it to a well-characterised viral outbreak, the spread of rabies in raccoons across the north-east of the usa over approximately thirty years. this data set has been analysed using previous phylogeographic methods, thus enabling direct comparison with the results obtained here. methods overview of methodology the structure of our framework can be summarised in the following five steps. a full description of each step is provided in the next section. in the first step, the history of lineage dispersal is recovered from one or more spatial- and temporally-referenced phylogenies. such trees are generated by the continuous phylogeography method implemented in beast, but it is important to note that our method is applicable phylogenies from any source that are annotated with dates and locations in the same way. the velocity, distance and duration of spatial movement along each branch in each tree are extracted and represented by a vector.fig. an illustration of the node position randomisation procedure used to generate null distributions of the d statistic. a the original environmental raster upon which is superimposed the movement events extracted from one spatiotemporally-referenced phylogeny. b the result of one randomisation of node positions. this randomisation procedure is performed within a minimum convex hull, which is defined by the node locations of all selected phylogenies in a second, optional, step, the spatio-temporal data extracted in step is used to calculate summary statistics of spatial spread, such as the velocity of the epidemic wavefront. each of the vectors obtained in step is assigned a “weight” score, which is calculated using a raster that defines the spatial heterogeneity of a specified environmental variable. we have implemented three different models of spatial movement to compute the “weight” allocated to each phylogeny branch for a given environmental raster. the “straight-line path” model, which assumes that movement is in a straight-line between the start and end locations of each branch. the “least-cost path” model, which uses a least-cost algorithm to determine the route taken between the start and end points. the “random walk path” model, which uses circuit theory to accommodate uncertainty in the route taken. note that for models and, each environmental raster must be considered twice, once as a conductance factor and once as a resistance factor. the correlation between the duration of each phylogeny branch and its “weight” is then estimated. correlations are repeated for each environmental raster, and for each of the three path models described above. finally, the statistical significance of these correlations is tested using a null model. to generate this null model, we implement a simple randomisation procedure: phylogenetic node positions are randomised within the study area, under the constraint that branch lengths, tree topology and root position are unchanged. description of methodology step 1: extracting spatio-temporal information from phylogenies the input data for our method consists of one or more spatiotemporally-referenced phylogenies. if multiple phylogenies are used then each is assumed to be equally probable. at present, such trees are mostly likely to be generated using the phylogeographic models implemented in beast, but our framework can be applied to phylogenies from any source that have been annotated in the same way. in order to extract the spatio-temporal information in the input trees, each phylogeny branch in each tree is summarised as a vector defined by its start and end location, and its start and end dates. each branch therefore represents an independent viral lineage dispersal event. vectors for each input tree are tabulated. step 2: estimation of dispersal and epidemiological statistics as an optional step, our analytical framework calculates several spatio-temporal statistics from these dispersal vectors. so far, we have implemented three such statistics: the mean lineage dispersal velocity, the mean lineage diffusion coefficient, and a time-series of the maximal epidemic wavefront distance. these statistics are visually summarised as plots of the kernel density of lineage dispersal velocity parameters, the kernel density of lineage diffusion coefficient parameters, and the change through time in the spatial and patristic maximal wavefront distances, as measured from the location of the tree root. the spatial distance corresponds to a straight-line distance, whereas the patristic distance equals the sum of the spatial distances along each phylogeny branch between the tree root and its tips. step 3: computation of environmental weights in this step, each of the vectors obtained in step is assigned a specific “weight” based on a raster of values that represent the type or magnitude of an environmental variable. in order to compute these “weights”, we have implemented three different “path models” which represent the path taken by a phylogeny branch as it travels between its start and end locations: straight-line path: each lineage travels in a straight line between its start and end location. in this case, the branch “weight” is computed as the sum of the values of the raster cells through which the straight line passes. least-cost path: each lineage travels via the least-cost path between its start and end location. with this model, the branch “weight” is computed as the sum of the transition values between adjacent cells along the least-cost path. random walk path: each lineage travels via a random walk between its start and end location. in this case, the “weight” is a graph-theoretic metric based on circuit theory, which takes into account multiple possible pathways connecting a given pair of locations, and the values of the raster cells through which they pass. specifically, for each lineage, one node is connected to a one-ampere current source, while the other is connected to ground. the “weight” between each pair of nodes is then defined by the effective electric resistance or conductance connecting them on the grid, i.e. the environmental raster in question. more details about the circuit theory model can be found in and. for the latter two models, an environmental variable can be treated either as a conductance or a resistance factor. for example, if an environmental raster is treated as a resistance factor then raster cells with low values will be more permeable to dispersal and those with high values will correspond to poor dispersal habitat or to movement barriers. when there is no prior information about whether a given variable will facilitate or impede lineage movement, it is prudent to consider it twice, once as a potential conductance factor, and once as a potential resistance factor. by default, we use the environmental values as resistance or conductance values. log-transformed and/or standardised values may be necessary if the current univariate approach is in future extended to a multivariate one. note, when all cells in a raster are identically-valued then the total “weight” for any path is simply proportional to its geographic distance. a special case exists when all cells have values equal to which we here define as the “null raster”. step 4: correlation analyses we next calculate the regression between the duration of each movement event and the “weight” computed for that branch. a separate regression is performed for each environmental factor of interest. the absolute strength of the regression can be measured using its coefficient of determination. specifically, we use the statistic d =, where r2env is the r obtained when branch durations are regressed against weights defined by the environmental raster, and r2null is the r obtained when branch durations are regressed against weights defined by the null raster. all environmental raster cell values are increased by to enable a direct comparison of r2env and r2null. d is therefore a correlation measure relative to a null hypothesis and represents the degree to which the regression is strengthened when spatial variation in the environmental variable is taken into account. step 5: significance testing using randomisation in the final step, the statistical significance of d is tested. in order to calculate a null distribution for the d statistic, we have implemented a randomisation procedure that randomises phylogenetic node positions under the constraint that branch lengths, the tree topology and the root position are unchanged. furthermore, this randomisation procedure is not applied to the entire raster but to a subset of it defined by the minimum convex hull around a set of locations. the set of locations comprises the start and end positions of each branch in the spatiotemporally-referenced input phylogenies. during the randomisation, if a randomised branch position falls outside the area defined by this convex hull, the algorithm randomises its position again until it does not fall outside this area. by constraining permutation within this convex hull we are, in effect, using the data to inform the relevant study area. in contrast, the square raster inputted by the user may be of arbitrary size and orientation. without this constraint, randomised phylogenies will fall in areas with different environmental values, potentially leading to type i errors. the result of this procedure is one p-value per input phylogeny, which equals the proportion of randomisation replicates that generated d values larger than that generated by the empirical input phylogeny. if multiple input trees are used then a distribution of p-values is obtained and we report the percentage of p-values > . it is worth noting that the statistical significance of our test is derived from a randomisation distribution, not from a distribution derived from the regression model. hence the coefficients of determination r2env and r2null are treated only as data statistics, and therefore our test is not dependent on regression model assumptions. application to rabies in north america the raccoon rabies virus data set comprises sequences with known sampling dates and locations. the sequences are ~2800 nt long and span the viral nucleoprotein n gene and the 5' end of the phosphoprotein p gene; see biek et al. for further details). we extracted the spatial-temporal information contained in spatiotemporally-referenced trees sampled from the posterior distribution of trees inferred for this data set by lemey et al.. we first calculated three epidemiological statistics from these input trees: mean lineage dispersal velocity, mean lineage diffusion coefficient and a plot of the maximal wavefront distance. for ease of explanation, we initially show how our analytical framework works on a single input phylogeny. subsequently we show how the method is extended to incorporate phylogenetic uncertainty by analysing multiple input trees. we investigated six environmental variables to determine if they were associated with the dispersal rate of raccoon rabies lineages. these included the three most important igbp land cover variables for the study area, as well as elevation, human population density and “inaccessibility”. the “forests” and “savannas” layers combine several igbp land cover layers. the environmental rasters are shown in fig.  the sources of the original raster files are given in table  the original data presented a resolution of arcmin, corresponding to cells ~1 km square. we generated distinct land cover rasters from the original data by creating lower resolution rasters whose cell values equalled the number of occurrences of each land cover category within the arcmin cells. the resolution of the three other original rasters was also decreased to arcmin for tractability. for each selected raster and path model combination, we performed randomisations using the randomisation procedure outlined in step above.fig. the six environmental variables that were tested in the analysis of the raccoon rabies virus data set. the region shown corresponds to the northeast of the usa, centered approximately on harrisburg, pa. details of the construction and source data for these rasters is provided in the main text www.igbp.net webmap.ornl.gov www.map.ox.ac.uk bioval.jr.ec.europea.eu performance on simulated data sets we used artificial data sets to measure the potential type i error rate of the randomisation procedure implemented to assess the significance of observed d values. artificial data sets were generated using four different approaches: the duration of each movement event was randomly permuted among branches, the duration of each movement event was equal to the spatial distance along the path between its start and end locations, the duration of each movement event was equal to the environmental “weight” assigned to the path between the locations, and the duration is equal to a value randomly and uniformly chosen between the values defined in cases and. case thus corresponds to case but with the addition of random noise. simulation approaches and generate data sets under the null hypothesis, whereas approaches and generate data sets under the alternate hypothesis. we generated artificial data sets for each of these three simulation approaches. the artificial data sets used the empirical phylogenies from the raccoon rabies data set, but replaced the empirical branch durations with new durations as specified above. for approach, environmental “weights” were calculated using the “elevation” raster treated as a resistance factor. this raster was chosen because elevation appears to impede the movement of raccoon rabies virus lineages. for each path model and simulation method, we report the percentage of artificial data sets for which the randomisation null hypothesis test p-value was < . for simulation approaches and, these percentages correspond to estimates of the type i error rate of the test. for approaches and, these percentages indicate statistical power.  we use the methods outlined above to analyse previously published virus gene sequences sampled from an epidemic of rabies virus in raccoons in north america. spatio-temporal phylogenies were estimated from these sequences using the continuous phylogeographic model implemented in beast. the analyses below were performed on phylogenies sampled regularly from the post burn-in posterior tree distribution generated by beast. estimation of spatio-temporal statistics we used our framework to calculate several statistics that summarise the spatio-temporal information inherent in the rabies virus phylogenies. these statistics include the epidemic wavefront and patristic distances through time, and parameters that describe dispersal velocities and diffusion coefficients. the epidemic wavefront plot shows that the spatial extent of the raccoon rabies epidemic increased at a relatively constant rate until around the patristic distance at this time is greater than the spatial distance. the difference between these two measures is to be expected, because the patristic distance includes movements in all directions, including those back towards the epidemic origin. hence dispersal was diffusive in nature and did not always follow the shortest or most direct path from epidemic origin to epidemic wavefront. kernel density graphs of the mean lineage velocity and phylogenetic diffusion coefficient indicate significant variation in these parameters among lineages. the estimated mean lineage velocity was ~37 km/year, similar to that reported for dengue virus in vietnam but substantially smaller than that of the west nile virus epidemic in north america. these statistics are helpful in quantifying the spatial dynamics of an epidemic and can also be used to compare different outbreaks. some, but not all, of them are implemented in the software package spread. we hope that future work will lead to the development and implementation of further summary statistics.fig. epidemiological statistics estimated from the raccoon rabies virus data set. a time-series of the spatial distance between epidemic origin and maximal epidemic wavefront, and b evolution of the patristic distance between epidemic origin and maximal epidemic wavefront, c kernel density estimates of lineage velocity parameters and d kernel density estimates of lineage diffusion coefficient parameters. in parts a and b the grey area corresponds to the 95 % credible region of the estimated wavefront position. in parts c and d the three contours show, in shades of decreasing darkness, the 25 %, 50 %, and 75 % highest posterior density regions via kernel density estimation impact of environmental variables next, we used our framework to measure and test the correlation of various environmental variables with the movement of rabies virus lineages. we explored rasters that represent six different environmental variables, specifically cropland land cover, forest land cover, savanna land cover, elevation, human population density and inaccessibility to major cities. for illustrative purposes, we first present the results based on only one tree. we later show how the framework can be extended to incorporate phylogenetic uncertainty, by analysing multiple input trees. the purpose of the analysis is to calculate a test statistic, d, for each environmental variable, e. the statistic d represents the degree to which the correlation between the duration and “weight” of each lineage is strengthened when spatial heterogeneity in e is taken into account. formally, d = r2env - r2null, where r2env is the coefficient of determination obtained when branch durations are regressed against branch weights calculated using raster e, and r2null is the coefficient of determination obtained when branch durations are regressed against branch weights defined by a null raster with no spatial variation. if d ≤  then the environmental variable e does not explain variation in branch duration better than geographic distance alone. however, if d is strongly positive, then the values of e are strongly associated with branches that move more rapidly or more slowly than average. figure  displays the linear regressions between branch durations and branch weights for a rabies virus phylogeny randomly chosen from the posterior tree distribution. regressions are shown for “weights” calculated using the null raster and for “weights” calculated using the elevation raster, when elevation is treated as a factor that impedes movement. results are shown for each of the three path models. in this example, there is comparatively little difference among path models in the strength of the relationship between branch durations and branch weights. for all three path models, r2env is approximately twice as great as r2null, so the corresponding d values are positive. thus while geographic distance alone explains a small part of among-lineage variation in velocity, this variation is better explained when the geographic elevation of the path taken by each lineage is taken into account. table  reports d values for other environmental variables when applied to the same single phylogeny as that used in fig.  fig. linear regressions between branch durations and branch weights for one phylogenetic tree. plots in were generated by calculating the branch weights using the “null” raster, whereas plots in were generated by calculating the branch weights using the “elevation” raster treated as a resistance factor. plots are shown for each of the three path models. the d value obtained under each path model is shown at the top. the plots show that the r of the regression is approximately doubled when spatial heterogeneity in elevation is taken into account for each combination of environmental variable and path model, the test is applied to a single phylogeny, and to a set of trees. for the former, the d statistic and p-value of the test are shown. for the latter, we report the percentage of trees for which p <  p-value <  in each case, artificial data sets were applied to the “elevation” raster treated as a resistance factor. values equal the percentage of artificial data sets for which the p-value of the null hypothesis test was <  the statistical significance of each d value can be assessed by comparing it to a null distribution generated by randomising the empirical node positions. in each case, randomisation replicates were performed. table  reports the p-values obtained when this randomisation procedure is applied to the rabies virus phylogeny introduced above. for this data set only some of the d values are statistically significant. specifically, significant positive d values were obtained for the “elevation” and “inaccessibility” rasters under the least-cost path model. although these results are interesting they are based on only one tree and therefore do not take in account the statistical uncertainty arising from phylogenetic inference. we will now demonstrate how the null hypothesis test is performed on phylogenies sampled from the posterior distribution of trees. first, we estimated d values for each of the trees. this generates, for each combination of environmental variable and path model, a distribution of d values. the randomisation procedure is then applied to each d value, precisely as outlined above. this results in p-values for each raster/path model combination. we therefore report the percentage of trees which give rise to p-values < ; this can be interpreted as the posterior probability of observing a significant correlation between lineage movements and the environmental variable. we illustrate this procedure in fig.  which compares the empirical distribution of d with five replicates of the null distribution of d generated by the randomisation procedure. in fig. 5a the distributions were calculated using the “elevation” raster and in fig. 5b they were calculated using the “forests” raster. in both cases the least-cost path model was used. in fig. 5a, the empirical distribution is greater than the randomisation replicates, indicating that the environmental variable in question explains viral lineage movement, whereas in fig. 5b there is no noticeable difference between the empirical and null d values.fig. empirical distributions of the d statistic, calculated from trees sampled using bayesian mcmc inference. these are compared with five replicates of the null distribution of d generated by the randomisation procedure. in the distributions were calculated using the “elevation” raster and in they were calculated using the “forests” raster. in both cases the least-cost path model was used. for visual clarity, discrete histograms were converted into density curves using a gaussian smoothing kernel the environmental variables in table  can be grouped into three categories: rasters that have little or no association with lineage movement, rasters for which there is weak evidence that they affect lineage movement, and rasters for which there relatively stronger evidence that the environmental variable is associated with slower or faster movements. the “forests” and “savannas” rasters fall within the first category, whilst “croplands” and “human population density” belong to the second category. only the “elevation” raster clearly belongs to the third category. the importance of the “inaccessibility” raster is less clear and is sensitive to the path model chosen. for both the “elevation” and “inaccessibility” rasters, the fraction of trees with p <  is smaller under the straight-line path model than under the other two path models. this may reflect the over-simplified nature of the straight-line model, which permits biologically unrealistic scenarios, such as paths that traverse large water bodies. performance on simulated data sets table  shows the results of null hypothesis tests performed on artificial data sets. these data sets were simulated using three different approaches. simulation approaches and generate artificial data sets under the null hypothesis. therefore for these two approaches the fraction of significant results should equal the critical value of the test. this is indeed the case, as only 0–7 % of tests on artificial data produce a p-value <  and hence the type i error of the test appears to be appropriate. simulation approaches and generates data sets under the alternate hypothesis, i.e. branch duration is equal or influenced by the environmental raster. the null hypothesis was almost always rejected when the test was applied to data generated using approaches or, indicating that the test has reasonable statistical power in the context of this data set. note that very similar results were also obtained on data sets whose simulations were based on alternative environmental rasters.  in this study we show that it is possible to integrate genomic and environmental data in order to test hypotheses concerning the mode and tempo of virus dispersal during emerging epidemics. the raccoon rabies data set explored here was chosen for illustrative purposes and our results strongly support the notion that increasing elevation is associated with slower movements of this virus. this result is biologically plausible and the impact of elevation on the dissemination of raccoon rabies was previously addressed by biek et al.. their analysis used spatial kriging to compute annual contours from the date of the first reported case of each county and the contours were coloured according to temporal periods identified by phylogenetic molecular clock analysis. they then visually compared the resulting contours overlaid on an elevation map, enabling them to posit that mountain ranges likely formed a barrier to raccoon movement. we intend to explore the application of our new method to other viral outbreaks in future work. we have implemented three path models in our analytical framework. although the straight-line path model is perhaps overly simplistic, the appropriateness of one model over another will depend on the locomotive behaviour of the pathogen’s host species. in the case of the raccoon rabies epidemic, we had no prior assumption about which path model might be most appropriate and therefore we decided to test them all. in other instances, there may be good reasons to choose a priori one path model over another. for instance, one might chose the least-cost path model over random walk path model in some cases, for example viruses carried by migratory birds; i.e. the least-cost model may be more appropriate when host-species intentionally reach specific locations by avoiding non-suitable landscape areas. simulations indicate that the randomisation procedure used to test the null hypothesis has appropriate type i error rates and acceptable statistical power. however, it is important to note that a significant result represents evidence only for a correlation between virus dispersal and an aspect of the geographic landscape. in future work we intend to explore whether generalized linear models can be applied to this framework in order to account for possible correlations among the environmental variables in question. in the meantime, the broader ecology and epidemiology of the pathogen, host, and habitat concerned should be taken into account when interpreting results. currently, our analytical framework does not allow for the analysis of environmental variables, such as temperature or humidity, that vary notably over the course of the epidemic under investigation. the assumption of constant values is unlikely to be unrealistic for the environmental variables investigated here, although some variation has undoubtedly occurred in human population density. the extension of our approach to time-varying environments is theoretically feasible, but poses significant technical and practical problems, specifically the acquisition of a series of environmental rasters that represents change in the variable concerned at a sufficiently precise temporal resolution, and generalisation of the least-cost and random-walk path models to three rather than two dimensions, with the third dimension corresponding to rasters that represent different points in time. a second improvement for future work would be to allow the co-analysis of several environmental factors in a multivariate framework. such a framework would require correlation statistics, for example those based on a generalised linear model, as in. implementation the methods introduced here are implemented in r and the scripts are freely available from http://evolve.zoo.ox.ac.uk/evolve/software.html. example files and a tutorial are available as additional file additional file additional file 1: example files and tutorial related to the present study. competing interests the authors declare that they have no competing interests. authors’ contributions sd, rr and ogp designed the study and the methodology. sd developed the analytical framework and performed the study. sd and ogp wrote the manuscript. all authors read, edited and approved the final manuscript. we thank two anonymous reviewers, roman biek, kenneth hoehn, nuno faria, philippe lemey and stéphane dray for their useful comments, as well as dan weiss, peter gething, clément tisseuil and marius gilbert for their advice regarding the processing and management of gis data. this work was funded by the wiener-anspach foundation. ogp received funding from the european research council under the european union's seventh framework programme / erc grant agreement no. 614725-pathphylodyn.
however, a remaining challenge is to develop an analytical framework that can test hypotheses about the effect of environmental variables on pathogen spatial spread. we present and apply a new statistical framework that combines genomic and geographic data to test the impact of environmental variables on the mode and tempo of pathogen dispersal during emerging epidemics. we test the impact of six different environmental variables on this epidemic and demonstrate that elevation is associated with a slower rabies spread in a natural population. recent phylogeographic approaches can reconstruct the history of virus dispersal from sampled viral genomes and infer the locations of ancestral infections. such methods provide a unique source of spatio-temporal information, and are exploited here. phylogenetic analysis is now an important tool in the study of viral outbreaks. this study shows that it is possible to integrate genomic and environmental data in order to test hypotheses concerning the mode and tempo of virus dispersal during emerging epidemics. a randomisation procedure is used to assess statistical confidence and we validate this approach using simulated data. next, tests measure the association between each environmental variable and lineage movement. we apply our new framework to a set of gene sequences from an epidemic of rabies virus in north american raccoons. phylogeographyphylodynamicsmolecular epidemiologyrabies virusissue-copyright-statement© the author 2016 the inferred dispersal path for each phylogenetic lineage is then assigned a “weight” using environmental data.
 phylogenetic analysis is now an important tool in the study of viral outbreaks. it can reconstruct epidemic history when surveillance epidemiology data are sparse, and can indicate transmission linkages among infections that may not otherwise be evident. however, a remaining challenge is to develop an analytical framework that can test hypotheses about the effect of environmental variables on pathogen spatial spread. recent phylogeographic approaches can reconstruct the history of virus dispersal from sampled viral genomes and infer the locations of ancestral infections. such methods provide a unique source of spatio-temporal information, and are exploited here.  we present and apply a new statistical framework that combines genomic and geographic data to test the impact of environmental variables on the mode and tempo of pathogen dispersal during emerging epidemics. first, the spatial history of an emerging pathogen is estimated using standard phylogeographic methods. the inferred dispersal path for each phylogenetic lineage is then assigned a “weight” using environmental data. next, tests measure the association between each environmental variable and lineage movement. a randomisation procedure is used to assess statistical confidence and we validate this approach using simulated data. we apply our new framework to a set of gene sequences from an epidemic of rabies virus in north american raccoons. we test the impact of six different environmental variables on this epidemic and demonstrate that elevation is associated with a slower rabies spread in a natural population.  this study shows that it is possible to integrate genomic and environmental data in order to test hypotheses concerning the mode and tempo of virus dispersal during emerging epidemics. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords phylogeographyphylodynamicsmolecular epidemiologyrabies virusissue-copyright-statement© the author 2016 