Comparative performance of viral landscape phylogeography approaches

Abstract

The rapid evolution of RNA viruses implies that their evolutionary and ecological processes occur on the same time scale. Genome sequences of these pathogens therefore can contain information about the processes that govern their transmission and dispersal. Landscape phylogeographic approaches use phylogeographic reconstructions to investigate the impact of environmental factors and variables on the spatial spread of viruses. Here, we extend and improve existing approaches and develop three novel landscape phylogeographic methods that can test the impact of continuous environmental factors on the diffusion velocity of viral lineages. In order to evaluate the different methods, we also implemented two simulation frameworks to test and compare their statistical performance. The results enable us to formulate clear guidelines for the use of three complementary landscape phylogeographic approaches that have sufficient statistical power and low rates of false positives. Our open-source methods are available to the cientific community and can be used to investigate the drivers of viral spread, with potential benefits for understanding virus epidemiology and designing tailored intervention strategies.

Keywords: continuous phylogeography; landscape phylogeography; molecular epidemiology; viruses.

Fig. 1.

Illustration of a cartogram and MDS transformations based on an environmental layer tested as a potential resistance factor slowing the diffusion velocity of viral lineages. We first display an original elevation raster (A; here rescaled from 1 to 11) along with the sampling coordinates (black dots) of the raccoon rabies virus (RABV) dataset sampled in northeastern USA (38), followed by a cartogram transformation (B) of this raster. We then also display an MDS transformation of those coordinates based on environmental distances computed with the Circuitscape algorithm among all pairs of sampling points while treating the elevation raster as a resistance factor (C). In this third panel, sampling points are colored according to the altitude of their position prior to the transformation of the space. As shown in panel B, the cartogram transformation of this elevation raster — here also tested as a potential resistance factor — tends to increase the relative distance between sampling points located in higher altitude areas while decreasing the relative distance between sampling points located in lower altitude areas. As shown in panel C, a similar observation can be made when visually inspecting the relative distance among sampling points before and after the MDS transformation: sampling points in relatively low altitude areas tend to be closer to each other after the MDS transformation. The rationale behind these two prior-informed landscape phylogeographic approaches based on a transformation of the study space (with a cartogram transformation or with a MDS analysis) is to assess whether such a transformation would lead to a more Brownian diffusion process than in the untransformed study space, which would then indicate an association between the tested environmental factor and the diffusion velocity of viral lineages.

Fig. 2.

Example of continuous phylogeographic simulations based on a birth–death process and a RRW diffusion process during which an environmental layer impacts the diffusion velocity of viral lineages. The two first panels display the phylogenetic tree obtained from a first simulation: (A) the time-scaled visualization of the simulated topology and (B) the mapped visualization on the environmental layer impacting the movement velocity of evolving lineages — here an elevation raster rescaled from 1 to 11 (see Materials and Methods for further detail) and acting as a resistance factor (i.e., slowing down movement velocity in relatively higher altitude areas). The third panel (C) displays the mapped visualization of a second example of simulation. In all panels, tree nodes are colored according to time, with internal and tip nodes colored according to their age and collection time, respectively. The RRW simulation process illustrated here is inspired by the phylogeographic analysis of the RABV spread in North America (13, 30, 38) and has been used in the present study to assess the statistical performance of four landscape phylogeographic approaches aiming at testing the impact of environmental factors on the diffusion velocity of viral lineages.

Fig. 3.

Simulation and analytical workflow implemented to assess the statistical performance of four landscape phylogeographic approaches to test the impact of environmental factors on the diffusion velocity of viral lineages or on the deviation from an IBD pattern. We here assessed the performance of two post hoc (IBR and lineage diffusion velocity analyses) and two prior-informed landscape phylogeographic approaches (analyses of continuous phylogeographic reconstructions based on sampling coordinates obtained either after a MDS or cartogram transformation). Continuous phylogeographic simulations illustrated for step 1 either consisted in (1.1) birth–death simulations based on a RRW diffusion process with a velocity impacted by an underlying environmental layer (referred to in the text as RRW simulations) or (1.2) tree topologies obtained by modifying the branch lengths of an actual MCC tree according to three distinct scenarios (MCC simulations). In the subsequent analytical steps illustrated in the graphical workflow, the path taken by those RRW and MCC simulations can then be followed by the reddish and bluish arrows, respectively. For the MCC simulations, three distinct scenarios were considered: no impact of the environmental factor (scenario 1; light blue arrows), as well as a moderate (scenario 2; intermediate blue arrows) and strong (scenario 3; dark blue arrows) impact of the environmental layer on lineage diffusion velocity. In step 2, three distinct versions of each simulated dataset were obtained by transforming or not the tip node coordinates: as illustrated by the different arrow paths, tip node coordinates were left untransformed, transformed using an MDS transformation (2.1), and transformed according to a cartogram transformation (2.2). (*) MDS transformations were conducted based on environmental distances among pairs of tip nodes, which were computed using either the least-cost (5) or Circuitscape (3, 32) path model while either considering the environmental raster (here an elevation raster with values rescaled between 1 and 11) or a corresponding null raster with accessible raster cell values uniformly equal to 1 (see the text for further detail). In step 3, a continuous phylogeographic analysis was conducted for each resulting simulated dataset. In the case of the RRW simulations, these analyses were directly based on the tree topologies and associated tip node coordinates obtained after step 2. In the case of the MCC simulations however, the continuous phylogeographic analyses were based on genomic sequences simulated along the tree topologies. Finally, step 4 consisted of the different landscape phylogeographic analyses developed and evaluated in the present study: two post hoc (IBR and lineage diffusion velocity analyses; 4.1) and two prior-informed landscape phylogeographic approaches (analyses of continuous phylogeographic reconstructions based on sampling coordinates obtained either after a MDS or cartogram transformation; 4.2). As detailed in the text, each of those four landscape phylogeographic approaches is based on the estimation of a distinct Q statistic: Q₁ for the lineage diffusion analyses, Q₂ for the IBR analyses, as well as Q₃ and Q₄ for the prior-informed analyses based on an MDS or cartogram transformation, respectively. All four Q statistics are defined as the difference between two coefficients of determination, R²_env and R²_null. For Q₁, R²_env is the coefficient of determination obtained from the linear regression between four times the branch durations t and the squared environmental distances D_env computed on the environmental raster (4t ~ D²_env), and R²_null the coefficient of determination obtained from the linear regression between four times the branch durations and the environmental distances D_null computed on the null raster (4t ~ D²_null). For Q₂, R²_env is obtained from the linear regression between the patristic distances and the log-transformed environmental distances computed on the environmental raster for each pair of tip nodes, and R²_null from the linear regression between patristic distances and log-transformed environmental distances computed on the corresponding null raster. For Q₃, R²_env is obtained from the linear regression between four times the branch durations t and the squared Euclidean distances d_Eucl,env traveled by each phylogenetic branch in the space transformed according to environmental distances computed on the environmental raster (4t ~ d²_Eucl,env), and R²_null from the linear regression between four times the branch durations t and the squared Euclidean distances d_Eucl,null traveled by each phylogenetic branch in the space transformed according to environmental distances computed on the corresponding null raster (4t ~ d²_Eucl,null). Finally, for Q₄, R²_env is obtained from the linear regression between four times the branch durations t and the squared Euclidean distances d_Eucl,env traveled by each phylogenetic branch in the space transformed according to a cartogram transformation, and R²_null from the linear regression between four times the branch durations t and the squared great-circle geographic distances d_gc traveled by each phylogeny branch in the nontransformed space (4t ~ d²_gc), i.e., following a standard phylogeographic reconstruction based on the nontransformed sampling coordinates. In step 4, for Q₁ and Q₂, a BF support was eventually computed by comparing the resulting posterior distribution obtained for the Q statistic (Q_posterior) with the distribution obtained according to a null dispersal model (Q_{randomizations}); and for Q₃ and Q₄, a BF support was computed by assessing the probability of getting a positive distribution of Q_posterior values (see the text for further detail). Continuous phylogeographic simulations (step 1), transformations of sampling coordinates (step 2), and landscape phylogeographic analyses (step 4) were conducted with the R package “seraphim” (41), the sequence simulations with the program πBUSS (42), and the continuous phylogeographic reconstructions (step 3) with the software package BEAST X version 1.10.5 (19).