Host ecology determines the dispersal patterns of a plant virus.

Since its isolation in 1966 in Kenya, rice yellow mottle virus (RYMV) has been reported throughout Africa resulting in one of the economically most important tropical plant emerging diseases. A thorough understanding of RYMV evolution and dispersal is critical to manage viral spread in tropical areas that heavily rely on agriculture for subsistence. Phylogenetic analyses have suggested a relatively recent expansion, perhaps driven by the intensification of agricultural practices, but this has not yet been examined in a coherent statistical framework. To gain insight into the historical spread of RYMV within Africa rice cultivations, we analyse a dataset of 300 coat protein gene sequences, sampled from East to West Africa over a 46-year period, using Bayesian evolutionary inference. Spatiotemporal reconstructions date the origin of RMYV back to 1852 (1791-1903) and confirm Tanzania as the most likely geographic origin. Following a single long-distance transmission event from East to West Africa, separate viral populations have been maintained for about a century. To identify the factors that shaped the RYMV distribution, we apply a generalised linear model (GLM) extension of discrete phylogenetic diffusion and provide strong support for distances measured on a rice connectivity landscape as the major determinant of RYMV spread. Phylogeographic estimates in continuous space further complement this by demonstrating more pronounced expansion dynamics in West Africa that are consistent with agricultural intensification and extensification. Taken together, our principled phylogeographic inference approach shows for the first time that host ecology dynamics have shaped the historical spread of a plant virus.

1. Introduction

Although phylodynamics have become a burgeoning area of research focused on many human and animal viruses, comparatively fewer analyses have targeted the interaction between evolutionary and ecological dynamics in plant viruses. On the one hand, this may be explained by a biased interest in viruses that directly impact human health or that may emerge as zoonotic pathogens. On the other hand, it is unclear to what extent phylodynamic concepts apply to plant viruses because their evolutionary and ecological dynamics may not necessarily occur on the same time scale. Lower rates of plant virus evolution have been inferred based on co-divergence assumptions, but also from sequence analysis of old samples (Rodríguez-Cerezo ; Fraile ; Gibbs ). In recent years however, evidence has accumulated for a rapid evolutionary rate in specific plant viruses, as first demonstrated for rice yellow mottle virus (RYMV) (Fargette ) and zucchini yellow mosaic virus (Simmons, Holmes, and Stephenson, 2008). Nucleotide substitution rates falling within the range of animal RNA virus rates have also been reported for particular Geminiviridae (Duffy and Holmes, 2008, 2009; Monjane ) and Luteoviridae (Pagán and Holmes, 2010).

In addition to clarifying the tempo and time scale of plant virus evolution, molecular sequence analyses may also probe spatial population structure and shed light on the transmission dynamics that gave rise to the current spatial distribution of plant viral lineages. It is therefore not surprising that the field of plant virus epidemiology has started to adopt recent statistical inference methodology that integrates temporal and spatial dynamics in a phylogenetic context (Lemey , 2010; Drummond ). As an example of this, the ongoing global spread of tomato yellow leaf curl virus (TYLCV) has attracted significant interest as a potential threat to tomato production in all temperate parts of the world. Motivated by the need to unravel the ecological and economic risks associated with such viral invasions (Lefeuvre ) applied Bayesian phylogeographic methods to reconstruct the spatiotemporal history of TYLCV spread and diversification. This revealed that, while the virus likely originated in the Middle East during the first half of the 20th century, this area remained epidemiologically relatively isolated. Instead, many global movements of TYLCV appear to have been seeded from the Mediterranean basin. As another example of a tropical plant virus that poses a threat to African food security, maize streak virus (MSV) has caused severe epidemics throughout the maize growing regions of Africa. Recent insights gained from Bayesian spatiotemporal reconstructions point at southern Africa as the most probable location from which MSV emerged at the beginning of the 20th century, and subsequently spread transcontinentally at an average rate of 32.5 km/year (Monjane ). As the etiological agent of the most damaging plant virus disease in the world, cassava mosaic-like virus (CMV) has caused devastating crop losses across sub-Saharan Africa. This epidemic was estimated to have originated in the late 1930s in mainland Africa with subsequent introductions to the southwest Indian ocean islands between 1988 and 2009 (De Bruyn ).

Among the fast evolving plant viruses, RYMV is also of particular interest because it circulates in most rice growing countries on the African continent (Bakker ; Abubakar ), impacting the lives of millions of impoverished Africans that rely on rice agriculture for subsistence (Abo, Sy, and Alegbejo, 1998). Symptoms of RYMV infection range from discolouration, stunting and ultimately sterilisation of the plant, resulting in devastating epidemics with yield losses that vary from 10 to 100% depending on how early the infection sets in, the type of rice cultivation and the rice cultivars used (Allarangaye ). RYMV is a member of the Sobemovirus genus with a genome composed of a single-stranded positive RNA molecule encompassing about 4450 nucleotides, organised into five open reading frames (ORFs) that overlap (except for ORF1) (Ling ; Sõmera, Sarmiento, and Truve, 2015). The virus is transmitted by chrysomelid beetles (Bakker ), by mammals (Sarra and Peters, 2003), and by contact during cultural practices (Traoré ), but no evidence of seed transmission has been found (Konate ). The known natural host range of RYMV is limited to the two species of cultivated rice Oryza sativa L. and Oryza glaberrima Steud, and a few related wild grasses (Bakker ).

Although the history of rice agriculture in Africa dates back many centuries, RYMV was only first reported in 1966 in Kenya (Bakker ). With nucleotide substitution rates ranging from 4 × 10−4 to 1.2 × 10–3 nucleotides/site/year, evolutionary studies have characterised the virus as a measurable evolving population with a most recent common ancestor (MRCA) dating back to around 1811 (Fargette ; Pinel-Galzi ). In addition to rapid evolutionary rates, specific RYMV gene sequences also show evidence for recombination between particular ORFs, but not within individual ORFs (Pinel-Galzi ). Early spatial genetic analyses have suggested a fairly regular pattern of spread with a correlation between genetic and geographic distances and no evidence of long-range dispersal. Based on comparisons of genetic diversity, these analyses have also implicated East Africa as the area of early diversification (Abubakar ). Specifically, more recent surveys confirm a large concentration of RYMV diversity in eastern Tanzania (Pinel-Galzi ), a region that is isolated by the Indian Ocean to the east and by the Eastern Arc Mountains to the west. A relatively long history of co-existence of RYMV strains in conditions that support habitat fragmentation indeed point at this region as a putative origin for the virus (Fargette ).

RYMV diversity shows a pronounced and characteristic geographic structure, and has been classified into S1–S6 strains based on serological typing and phylogenetics. Five serological profiles have been identified: three in West and Central Africa (Ser1, Ser2, and Ser3) and two in East Africa (Ser4 and Ser5). Apart from Ser5, which is divided into the S5 and S6 strains, these serotypes also correspond to the S1–S6 strain types (Pinel ; Traore ). The spatial structure of the epidemic, with different strains circulating in different countries, suggests a relatively recent expansion perhaps driven by the intensification of agricultural practices (Konaté and Fargette, 2001; Abubakar ). Evolutionary and ecological hypotheses about the origin and spread of RYMV have however not been examined in a coherent statistical framework. Recent extensions of Bayesian phylogenetic diffusion models for discrete traits now offer the opportunity to formally evaluate predictors of spatial spread. In particular, the recently developed GLM approach parameterises rates of diffusion as a function of potential predictors (Lemey ). This approach has for example identified human and animal transportation measures as the drivers of spatial spread for different influenza viruses (Lemey ; Nelson ), and it may also be useful for identifying the factors responsible for plant virus spread.

Here, we demonstrate the value of Bayesian phylodynamic inference methodologies in plant molecular epidemiology by focusing on the patterns of RYMV spread across Africa and reconstructing its phylogeographic history. We test spatiotemporal hypotheses about the origins of RYMV using state-of-the-art Bayesian statistical inference, quantify the dynamics of spatial spread in both East and West Africa, and formally assess the relationship between RYMV spread and the history of rice cultivation in Africa.

2. Methodology

2.1 Dataset compilation

We have assembled a RYMV sequence dataset by retrieving all publicly available ORF4 CP gene sequences from GenBank (on 4 September 2012) and combining these with additional samples made available by collaborators, which have now been published in recent studies (Hubert ; Longué ). The sequences were aligned using MAFFT version 6.864 b (Katoh and Toh, 2008) and manually edited in Se-Al (tree.bio.ed. ac.uk/software/seal). The final dataset consists of 300 sequences that were sampled between 1966 and 2012 in 20 countries across East and West Africa (Supplementary Fig. S1) and covers all countries in which RYMV has been reported. Although we used countries as locations in the discrete analyses, more specific location coordinates for all the 300 isolates were available and these were used in the continuous phylogeographic reconstructions.

To evaluate the impact of sampling bias on the root location estimate, we applied two different subsampling procedures to the sequence data. Specifically, we subsampled Tanzanian samples down to the next highest sampled country (Côte d’Ivoire) by: (1) randomly selecting 51 Tanzanian isolates and (2) selecting the same number of sequences that best represent the Tanzanian RYMV diversity using the Phylogenetic Diversity Analyzer tool (www.cibiv.at/software/pda) Minh, Klaere, and von Haeseler (2006, 2009). Both down-sampling procedures resulted in datasets of 266 sequences.

To contextualise particular continuous phylogeographic estimates, we also assembled two datasets for MSV and East African CMV (EACMV), both imposing an enormous burden to crops worldwide and in Africa specifically. For MSV, we resorted to the 333 full-genome recombinant-free dataset of Monjane . For EACMV, we obtained a dataset comprising 65 full-genomes from De Bruyn by focusing on non-recombinant sequences originating from mainland Africa by and excluding outliers in a linear regression analysis of root-to-tip divergence (see ‘Temporal signal’ section below).

2.2 Temporal signal

In order to visually examine the degree of temporal signal—or signal for divergence accumulation over the sampling time interval—in the RYMV CP sequence data, we employed an exploratory linear regression approach. We first followed a standard approach of estimating a maximum likelihood (ML) tree under a non-clock (unconstrained) generalised time-reversible (GTR) substitution model with discrete Γ-distributed rate variation among sites using PhyML (Guindon ) and plotted the root-to-tip divergences as a function of sampling time according to a rooting that maximises the Pearson product-moment correlation coefficient using Path-O-Gen (tree.bio.ed. ac.uk/software/pathogen). For comparison, we only plot the root-to-tip divergences for the same subset of taxa that is used in the procedure discussed below.

We also explored an alternative approach that attempts to avoid rate heterogeneity imposed on the deep branches connecting different RYMV clusters, and only considers the overall divergence accumulation within these specific clusters. To identify rooted clusters, we used the maximum clade credibility (MCC) tree from the Bayesian analyses (see below) and re-estimated branch lengths under a non-clock GTR + Γ substitution model using PAUP* v4.0b10. We then selected distinct phylogenetic clusters that contain taxa with a minimum sampling time interval of 15 years and that were associated with a posterior probability support >0.75 (Supplementary Fig. S2). For each of these clusters, we obtained cluster-specific MRCA-to-tip divergences as a function of sampling time based on the branch lengths estimated under the non-clock model. To explicitly model a cluster effect in the MRCA-to-tip divergence data dij for taxon i assigned to cluster j, we fit the following regression model:

where β is the intercept for cluster j, δ is the phylogenetically unadjusted rate of substitution, xij is the sampling time and ϵ an independent 0-mean error term. To visually plot a single regression line through the MRCA-to-tip divergence data with δ as slope, we subtract the estimated cluster effect from the divergence measurements and plot the resulting values, , as a function of sampling time. We acknowledge that linear regression techniques are not appropriate estimators of divergence through time as sequences do not represent independent data, but we merely employ these approaches to tentatively examine temporal signal in our data.

In addition to the visual exploration, we also conducted a date-randomisation test to evaluate to what extent Bayesian evolutionary rate estimates (using the Bayesian Evolutionary Analysis Sampling Trees (BEAST) package, see below) from the time-stamped data deviate significantly from estimates based on randomised tip dates, for which no particular relationship between sampling time and root-to-tip divergence is expected (Firth ). For this purpose, we here propose a novel implementation of this test that avoids having to analyse multiple date-randomised datasets. Because of the computationally expensive nature of Bayesian phylogenetic analyses, the number of these randomisations is generally limited (e.g. 20) in the standard test procedure (Duchêne ). Our BEAST implementation makes use of novel transition kernels that effectively randomise dates during the Markov chain Monte Carlo (MCMC) sampling procedure. Therefore, we do not have to rely on a number of specific randomisations, but conveniently, we average over all possible randomisations in a single analysis. We follow Duchêne and use as our criterion for a significant temporal signal that the 95% credible interval (CI) for the rate estimate obtained from correct sampling times should not overlap with the CI for the estimate obtained while randomising sampling times.

2.3 Bayesian evolutionary inference

We reconstructed time-calibrated phylogenetic and phylogeographic histories using a Bayesian statistical framework implemented in the software package BEAST v1.8 (Drummond ). BEAST uses MCMC integration to average over tree space, so that each tree is weighted proportional to its posterior probability. All analyses were performed using the Broad-platform Evolutionary Analysis General Likelihood Evaluator (BEAGLE) library to enhance computation speed (Suchard and Rambaut, 2009; Ayres ).

2.4 Sequence evolution

To model the nucleotide substitution process, we partitioned the codon positions into first+second and third positions (Shapiro, Rambaut, and Drummond, 2006) and applied a separate Hasegawa–Kishino–Yano 85 (HKY85) substitution model (Hasegawa, Kishino, and Yano, 1985) to the two partitions, each with a discretised Γ distribution () to model rate heterogeneity across sites. To accommodate among-lineage rate variation we applied an uncorrelated relaxed molecular clock that models branch rate variation according to a lognormal distribution (Drummond ). To investigate the sensitivity of the time to MRCA (TMRCA) estimates with respect to the coalescent prior, we tested all currently available flexible, non-parametric demographic priors: the skyride (Minin ) (using uniform smoothing over all inter-coalescent intervals), skyline (Drummond ), and skygrid (Gill ) (using a cut-off to 200 years with 100 grid points) model. Whereas it is generally recommended to employ ‘time-aware’ smoothing for the skyride model, which weighs the smoothing such that the effective population size changes between small, consecutive inter-coalescent intervals are penalised more than changes between intervals of larger size (Minin ), this appeared problematic for our dataset without strong temporal signal and resulted in MRCA estimates that were close to the oldest sample. We ran three independent runs for 100 million generations, sampling every 10 000th and discarded 10% as the chain burn-in. Stationarity and mixing (e.g. based on effective sample sizes ≥ 200 for the continuous parameters) were examined using Tracer version 1.6 (tree.bio.ed.ac.uk/software/tracer), and MCC trees were summarised using TreeAnnotator.

2.6 Discrete phylogeography

We modelled discrete location transitioning of RYMV between the 20 African countries throughout the phylogenetic history using both a reversible and non-reversible continuous-time Markov chain (CTMC) process (Lemey ; Edwards ) and performed the analysis with and without a Bayesian stochastic search variable selection (BSSVS) procedure to identify a sparse migration graph, which includes a restricted number of non-zero rates in the CTMC matrix. These analyses were performed both on the full dataset and the two subsampled datasets. We evaluated model fit using (log) marginal likelihood estimates obtained through path sampling (Lartillot and Philippe, 2006) and stepping-stone sampling (Xie ) procedures as implemented in BEAST (Baele , 2013; Baele and Lemey, 2013). We ran various computational settings to assess convergence of the (log) marginal likelihood estimates. The number of location transitions (‘Markov jumps’) and the time spent in each location state (‘Markov rewards’) were estimated using stochastic mapping techniques (Minin and Suchard, 2008a,b).

In order to quantify the spatial structure, we measured the phylogenetic association in the location trait data by applying the association index (AI) to our posterior set of trees (Wang ; Lemey ). This metric quantifies the degree to which the same traits tend to cluster together relative to the expectation for randomised trait assignments. AI values close to 0 reflect strong phylogeny-location correlation whereas AI values close to 1 reflect the absence of phylogenetic structure for the trait (Wang ; Lemey ).

For both the discrete as well as the continuous phylogeographic analysis (cfr. below), we use TreeAnnotator (Drummond ) to summarise the location estimates on a MCC tree and visualise the tree with annotations using FigTree (tree. bio. ed.ac.uk/software/figtree). We converted the location-annotated trees to keyhole markup language format using the Spatial Phylogenetic Reconstruction of Evolutionary Dynamics software package (Bielejec ) and visualise the spatial projections using Cartographica (www.macgis.com). We also used GenGIS (Parks ) to visualise the MCC tree as a tanglegram in a map adapted from Natural Earth (www.naturalearthdata.com).

In order to test the contribution of various predictors to the patterns of RYMV spread, we adopted a recent GLM extension of discrete phylogeographic diffusion (Lemey ). This approach models diffusion rates as a log linear function of a number of explanatory variables, and performs Bayesian model averaging to identify the combination of variables that is predictive of spatial spread while simultaneously reconstructing the phylogeographic history. The support and effect size for each predictor is estimated using inclusion probabilities and GLM coefficients, respectively (Lemey ). We considered the following predictors in our GLM-diffusion model: (1) great-circle distances between the centroids of each pair of countries; (2) intensities of rice cultivation by country (area of cultivated rice divided by the total country area (hectares per year)) at two different time points (1960 and 1990, obtained from faostat3.fao.org); (3) spatially disaggregated rice production statistics (area harvested in hectares) around the year 2000, obtained using the Spatial Production Allocation Model (HarvestedChoice, 2011) (Supplementary Fig. S3; since this is expressed in hectares of harvested rice for a 5-arc minute grid cell, we consider this as a measure of host connectivity); (4) precipitation by country in millimetres per year (www.climatemps.com); and (5) sample sizes (number of sequences included per country).

Because we model predictors of diffusion rates between pairs of locations, we include both an ‘origin’ and ‘destination’ predictor for location-specific measures such as intensity of rice cultivation, precipitation and sample size. In order to derive pairwise predictor values from the spatially mapped rice production statistics, we employ circuit theory to measure distances on a heterogeneous landscape, with rice area harvested as the heterogeneity factor. Specifically, we use Circuitscape version 3.5 (Shah and McRae, 2008) to compute the distances among pairs of locations in the rice production landscape based on a map of about 320 000 cells that encompasses all 20 sampling regions from East to West Africa. Cells with lower area of rice harvested provide higher resistance than cells with higher rice production. The landscape therefore represents a resistance surface that models small distances between nearby locations that are separated by high rice production and large distances between distant locations that are separated by low rice production. We estimated distances between all pairs of sampling regions and chose to connect cells to their eight neighbouring cells, not only to connect cells to their four cardinal neighbours but also to connect diagonally adjacent cells (Shah and McRae, 2008). All predictors were log transformed and standardised prior to their inclusion in the GLM analyses. We follow Lemey and specify prior inclusion probabilities that put 50% prior probability on no predictor being included, and a normal prior with a mean of 0 and a standard deviation of 2 on the coefficients in log space. Bayes factor (BF) support for predictors was calculated based on the ratio of posterior to prior odds for predictor inclusion.

2.6 Continuous phylogeography

To study the geographic spread of RYMV in continuous space and quantify its tempo of dispersal, we used a phylogenetic Brownian diffusion approach that models the change in coordinates (latitude and longitude) along each branch in the evolutionary history as a bivariate normal random deviate (Lemey ). As an alternative to homogeneous Brownian motion, we adopt a relaxed random walk (RRW) extension that models branch-specific variation in dispersal rates similar to uncorrelated relaxed clock approaches (Drummond ; Lemey ). Specifically, we independently draw branch-specific scalers of the RRW precision from a log-normal distribution to relax the assumption of a constant precision (=1/variance) among branches (Lemey ). The original implementation of multivariate diffusion models in BEAST (Lemey ) resorted to data augmentation of the unobserved locations of ancestral nodes in the phylogeny to compute multivariate trait likelihoods. Here, we employ a more recent dynamic-programming approach that integrates over all possible realisations of the unobserved traits (Pybus ), and provides a more tractable, efficient, and stable inference for large datasets with considerable diffusion rate heterogeneity.

Bayesian estimates under continuous diffusion models yield a posterior distribution of phylogenetic trees, each having ancestral nodes annotated with location estimates. To quantify the spatial epidemic dynamics, we summarise several statistics from the posterior estimates of the continuous phylogenetic diffusion process, as previously introduced by Pybus . Specifically, we provide mean posterior estimates and 95% highest posterior density (HPD) intervals for: (1) dispersal rate (km/year), summarised as the total great-circle distance traveled across the phylogenetic branches divided by the total time elapsed on the branches; (2) wavefront rate (km/year), summarised as the largest great-circle distance traveled from the root location estimate divided by the time since the MRCA; and (3) diffusion coefficient (km2/year), which reflects the diffusivity or the area that an infected host explores per time unit. Here, we use a ‘weighted average’ alternative of the diffusion coefficient () introduced by Pybus because this has recently been shown to provide estimates with considerably lower variances (Trovão ). This statistic is defined as follows (Trovão ): where g and t represent the great-circle distance and time, respectively, along branch of the random phylogeny.

3. Results

3.1 Evolutionary rate and divergence time estimation

As a standard check prior to fitting dated-tip molecular clock models, we first explored to what extent our dataset contained visually-detectable signal for sequence divergence throughout the sampling time interval. Despite the fact that previous evolutionary rate estimates for RYMV are fairly consistent (Fargette ), our standard linear regression exploration of root-to-tip distances as a function of sampling time did not reveal clear evidence for temporal signal in the complete dataset (Fig. 1A).

Figure 1.

Root-to-tip divergence as a function of sampling time for ML tree clusters (A) and for MCC tree clusters after removing the deep branch effects (B). Colour-coding identifies the 5 different clusters included.

We therefore hypothesised that clusters of more closely related variants may still contain temporal information, but extensive rate heterogeneity along the deeper branches connecting these clusters may confound visual detection of such signal. Indeed, particular clusters in the rooted tree with branch lengths estimated using an unconstrained (non-clock) model (Fig. 1A and Supplementary Fig. S2), have tips that are systematically more divergent from the root than other clusters. To examine the impact of this on root-to-tip divergence as a function of sampling time, we perform a similar analysis based on MRCA-to-tip divergences as a function of sampling time for specific clusters and level-out differences in cluster heights prior to plotting all the divergence data (cfr. ‘Methods’ section and Supplementary Fig. S2). This effectively ignores the rate heterogeneity on the deeper branching and removes the cluster effects on the root-tip-regression (Fig. 1B), resulting in a somewhat more discernible divergence accumulation through time. Together with an improved fit (adjusted R2 increase from 0 to 0.12), this suggests that the rate variation among the deeper branches can indeed affect the temporal signal estimate for the complete dataset.

The presence of temporal signal remains however questionable, which urged us to complement this exploration with a date-randomisation test implemented in BEAST (cfr. ‘Methods’ section). Although the rate was estimated to be () nucleotide substitutions per year per site by BEAST using the correct sampling dates, averaging over all possible date randomisations resulted in a far lower estimate of (). Given that the 95% HPDs do not overlap for these estimates, we follow Duchêne in considering this as evidence for significant temporal signal in the time-stamped data. We estimated these rates using a relaxed molecular clock model (Drummond ), which is better supported by the data than a strict clock model (see Supplementary Table S1). This is perhaps not surprising given the lack of a clear divergence accumulation over the sampling time interval in the exploratory linear regression analyses (Fig. 1A), and the substantial variation of the rate about its mean (coefficient of variation = 0.75; Table 1).

Table 1.

Impact of coalescent model on the TMRCA estimate

	Date for the MRCAa (year)	Evolutionary rate (μ) (substitutions/site/year)	Coefficient of variation for μ
Skyline	1818 [1727–1889]	9.71 × 10−4 [7.18 × 10−4–1.25 × 10−3]	0.76 [0.58–0.96]
Skyride	1818 [1736–1891]	9.39 × 10−4 [7.01 × 10−4–1.21 × 10−3]	0.75 [0.58–0.95]
Skygrid	1852 [1791–1903]	1.01 × 10−3 [7.51 × 10−4–1.32 × 10−3]	0.78 [0.59–0.99]

Values in between brackets represent 95% HPD intervals.

aMRCA, most recent common ancestor.

Because the absence of strong temporal signal may lead to a more pronounced impact of tree priors on divergence time estimates, we estimated TMRCAs using three different flexible non-parametric approaches (Table 1). Whereas the rate and TMRCA estimates under the skyline (Drummond ) and skyride (Minin ) models are very similar, the skygrid model results in a somewhat higher rate and younger TMRCA estimate (Table 1), but HPDs remain widely overlapping for estimates under the different models. We note that the estimates under the skyride model were sensitive to the way population sizes estimates are smoothed across the evolutionary time scale (cfr. ‘Methods’ section), which is likely due to the lack of strong temporal signal. As previously shown through simulations (Gill ), the skygrid model performs the best for divergent time estimates and the substitution rate under this model is also more consistent with previous studies (Fargette ). We therefore use this coalescent tree prior in all further analyses.

3.2 Discrete geography

By mapping the tip locations of a cladogram representation of the MCC tree in geographic space (Fig. 2), we highlight a clear separation of East and West African RYMV diversity. The modal location state estimates obtained by discrete phylogeographic reconstruction (represented by branch colours in the tanglegram and in the equivalent time-measured tree in Fig. 3), also reveal a strongly spatially structured viral population. We quantified this through the degree of phylogenetic clustering by location as summarised using the AI (Wang ), and found a low AI of 0.109 (0.08–0.14) indicating that the degree of spatial structuring is not so far from absolute (AI = 0).

Figure 2.

Tanglegram representation of the RYMV history by mapping the tip locations of the MCC cladogram to the geographic location of sampling. Branches are coloured according to the modal location state estimates obtained by a discrete phylogeographic reconstruction under a reversible CTMC model with BSSVS.

Figure 3.

Time-calibrated MCC tree inferred for 300 CP sequences of RYMV. Branches are coloured according to the most probable location state, indicated in the coloured legend (A). Posterior probability densities of the root location state for discrete reversible model with BSSVS; 66.87% of the posterior mass for the root location supports Tanzania as the origin location of RYMV epidemics (B). Early separation of the East-West epidemics in RYMV history, with 87 and 43% of posterior mass for the West (green) and East (violet) lineage clades, respectively.

To infer the discrete ancestral location states, we applied both reversible and non-reversible discrete diffusion models with and without a BSSVS procedure (Lemey ; Edwards ) and compared model fit for the four combinations using (log) marginal likelihood estimation (Baele , 2013; Baele and Lemey, 2013). Although we report the results for the best fitting model (reversible with BSSVS, see Supplementary Table S2), we note that the ancestral reconstructions are robust with respect to diffusion model specification. For example, all four model combinations find support for Tanzania as the geographical origin of RYMV (Fig. 3B), and this remained the best supported root location when the Tanzanian strains were downsampled to the same number as for Côte d’Ivoire (see Supplementary Table S3). The support for Tanzania emerges from a relatively high RYMV diversity in this country, encompassing most of the diversity in the East Africa clade (lineage S4, S5, and S6 in Fig. 3A), and hence a strong support for this location state at ancestral nodes up to the root node. West Africa (lineage S1–S3) was seeded relatively early in the RYMV history as its MRCA dates back to 1887 (1840–1919) and was estimated to have originated from Côte d’Ivoire.

Using BSSVS (Lemey ), we quantify the support for different diffusion pathways under the form of BF support for non-zero rates. Not surprising, we find support for a separate East and West African diffusion network (Fig. 4). We complement the support for the rates by estimating the number of transitions that occurred between the states involved using Markov jump counting (Minin and Suchard, 2008a) (Fig. 4). In East Africa, we find support for diffusion out of Tanzania to Kenya, Uganda, and Rwanda, and from the latter country also to Burundi and the Democratic Republic of Congo (DRC). The western diffusion network involves more countries and is characterised by a high degree of seeding from Côte d’Ivoire with diffusion pathways that extend eventually to Central Africa (the Central African Republic). Taken together, the diffusion pathways in East and West Africa we display in Figure 4 account for 84% of the location state transitions recovered in the RYMV evolutionary history.

Figure 4.

BF test support for discrete diffusion rates. Rates supported by a BF > 4 are indicated. The line colour represents the relative strength by which the rates are supported: green lines and red lines suggest relatively weak and strong support, respectively. The thickness of the arrows indicates increasing number of Markov jumps between locations.

To investigate what process of RYMV spread has led to the spatial genetic patterns we describe here, we apply a recent GLM extension of the discrete phylogeographic diffusion model that allows to test different potential predictors of viral dispersal (Faria ; Lemey ). We consider geographic distances, distances measured on a resistance landscape of harvested rice area (cfr. ‘Methods’ section and Supplementary Fig. S3), location-specific rice intensities at two different time points and precipitation as potential explanatory variables of the patterns of spread (Fig. 5). To examine whether the predictor support is robust to sample size heterogeneity we also include sample sizes as an explanatory variable. The GLM procedure, which attempts to identify the linear combination of predictors of spatial diffusion while reconstructing the phylogeographic history, finds maximal support for distances measured on a landscape of harvested rice area as a predictor of RYMV dispersal (Fig. 5). This predictor has a negative log effect size implying an inverse relationship with transition rates. That is, a high distance in the resistance landscape, as reflected by a large geographic distance and/or low harvested rice area between these locations, correlates with less intense viral dispersal. The importance of harvested rice intensity is also reinforced by a relative modest additional support for origin rice intensity in 1990 (BF =14.2), which is accompanied by a positive log effect size, suggesting higher viral dispersal out of locations with higher rice intensity. In addition to a host connectivity component, the distances in the harvested rice area landscape also incorporate geography. The fact that both are important is demonstrated by a GLM analysis that excludes the landscape distance predictor, which results in clear support for distance as well as origin rice intensity in 1990 (Supplementary Table S4). No other predictor yielded noticeable support in our analyses, and remarkably, also sample sizes did not help to explain viral diffusion intensities (Fig. 5). By repeating the analysis separately on the East and West african clade, we demonstrate that the signal for predictor support can be entirely attributed to the more pronounced and dynamic West African spread. Whereas highly similar predictor support and effect sizes are obtained for West Africa, none of the predictors yield noticeable support in East Africa (Supplementary Fig. S4).

Figure 5.

Predictors of RYMV dispersal across Africa. For each potential predictor, the BF support and the conditional effect size obtained using the GLM diffusion approach implemented in BEAST are shown (posterior mean and 95% Bayesian CI). Note that the credibility intervals for the cES of the predictors with BF > 14 exclude zero, which can be considered as additional evidence for its importance.

3.3 Continuous phylogeography

In order to quantify the dynamics of RYMV spread in continuous space, we also applied multivariate phylogenetic diffusion models to the CP sequences and their geographic coordinates (Lemey ). Because of the clear separation between the spread of East and West African RYMV lineages (Figs. 2 and 4), we perform separate analyses on the data for both regions. We tested a model of strict Brownian diffusion against several versions of RRW models using marginal likelihood estimation (Baele , 2013), and found that a lognormal-RRW provided the best fit to the dispersal dynamics (Supplementary Table S5).

The spatiotemporal patterns of spread under this model are summarised in Figure 6. In agreement with the discrete phylogeographic results, the eastern root location is estimated in Tanzania, and from the early 1930s, the virus spreads in the direction of Kenya and Uganda. The viral expansion continues within these countries and finally, by 2012, the spread of RYMV also includes Burundi, Rwanda and the DRC (Ndikumana ; Hubert ). In West Africa, the credible contour for the origin location overlaps mostly with Côte d’Ivoire, Mali and Senegal from where the virus spreads to the south, west and east. By 1932, the eastward expansion includes Nigeria and continues towards Chad. Extensive diffusion dynamics further develop within and between Côte d’Ivoire and Mali, and also western locations including Sierra Leone and Guinea appear to be seeded from these countries, respectively. Recent years mark the arrival in the most eastern location (Central African Republic) (Longué ).

Figure 6.

Reconstruction of the continuous spatiotemporal dispersal of RYMV in West and East Africa, shown from 1852 to 2012 at intervals that capture the major dispersal events. Black lines show a spatial projection of the representative phylogeny. Coloured clouds represent statistical uncertainty in the estimated locations of RYMV internal nodes (95% HPD intervals).

Based on quantitative summaries of the eastern and western dynamics (listed in Table 2), both datasets are characterised by similar dispersal rates and diffusion coefficients with overlapping HPDs. In terms of wavefront dynamics; however, we observe a slower invasion rate in the East as compared with its western counterpart. In line with different invasion rates on similar time-scales, our phylogeographic reconstruction estimated East and West wavefront distances of ∼1025 (95% HPD: 743–1,239) km and 2,869 (2,370–3,356) km, respectively. By plotting how these wavefront distances evolved over time (Fig. 7), we show that from the mid-1800s viral expansion begins in both geographical regions, but whereas it levels off at around 1960 in the east, RYMV continued to expand its spread in West Africa. Similar continuous diffusion statistics for the East and West Africa dynamics (Table 2) indicate that these account for a considerable degree of heterogeneity in the spatial spread dynamics. For instance, the dispersal rate and diffusivity are similar in both East and West Africa whereas the wavefront rate is three times lower in East than in West Africa, possibly because of the numerous barriers to spread in East Africa, whereas the Niger-Bénoué river axis in West Africa may have been an efficient means of viral propagation. We note that RYMV, MSV, and EACMV are characterised by dispersal statistics that are generally within the same order of magnitude, although some statistics suggest more pronounced dynamics for the latter two. This might be explained by transmission through leafhoppers for MSV and whiteflies for EACMV, but also human-mediated dispersal through infected cuttings.

Table 2.

Dispersal and ecological parameters

	Dispersal rate (km/year)	Diffusivity (km2/year)	Wavefront rate (km/year)
West	13.08 [10.25–16.03]	1559.25 [1186.13–1976.88]	23.13 [14.33–31.85]
East	16.16 [11.55–21.19]	1595.53 [1106.99–2151.71]	7.51 [3.80–12.32]
MSV	33.17 [28.70–7.82]	11665.18 [9133.97–14979.72]	74.79 [45.25–109.36]
EACMV	13.30 [8.07–0.65]	3748.74 [2196.90–5890.32]	32.45 [15.67–56.12]

Values in between brackets represent 95% HPD intervals.

Figure 7.

Mean wavefront distances for the West (blue) and East (pink) epidemics. Mean values indicated by darker lines and 95% HPD intervals indicated by coloured shadows.

4 Discussion and conclusion

As the main viral disease of rice, RYMV has been reported in all major rice producing countries in sub-Saharan Africa. Yielding losses up to 100%, it represents one of the tropical plant emergent diseases with the highest socio-economical impact (Fargette ). Evolutionary studies have only relatively recently characterised RYMV as a rapidly evolving plant virus based on heterochronous sequence data. In this study, we expand on the work of Fargette et al. (2008a,b) and Abubakar by reconstructing the RYMV phylogeographic history in both discrete (Lemey ) and continuous (Lemey ) space using Bayesian inference, and specifically test and quantify a range of potential predictors of spatial spread (Lemey ).

Although our RYMV evolutionary rate estimate is consistent with previous studies (Fargette ), it remains difficult to clearly detect accumulation of sequence divergence over the sampling time interval in the currently available data. We note that such temporal signal depends on both the evolutionary process and how we are able to sample from this process. A high overall tempo of evolution and a constant pattern of substitution accumulation will both increase the probability of measurable evolution over a particular time interval. A large temporal spread in sampling dates and more homogenous sampling throughout this interval will further contribute to the temporal signal (Seo ). An average RYMV substitution rate of about 0.001 substitutions per site per year and a sampling time interval of 44 years—even if sampling is more dense towards the present as is generally the case—may provide a relatively good opportunity to detect temporal signal in the CP gene. Substitution rate variability, however, appears to be a major confounding factor for RYMV, in particular because this may have acted on a relatively long evolutionary time scale of about 160 years. This is apparent through the consistently higher or lower tip divergences from particular clusters in the RYMV phylogeny (Fig. 1). Several factors may be responsible for rate heterogeneity in the phylogenetic history, including variation in mutation rate and replication rate as well as variation in selective pressure and host population sizes, but it remains difficult to disentangle these factors for RYMV. In general, intrinsic differences in mutation rate and replication dynamics are more likely to act between more distantly related viruses, such as different viral families. For more closely related viruses, host factors have been shown to impact viral evolutionary rates (Streicker ; Worobey, Han, and Rambaut, 2014). This together with the dynamics that impact the fixation of substitutions, represent interesting subjects for further RYMV research.

Even if exploratory linear regression plots do not suggest clear temporal signal, tip calibration may still prove useful if rate variation among branches is satisfactorily accommodated (Firth ). Relaxed molecular clocks may indeed perform reasonably well in modelling rate along the relatively long branches that separate distinct RYMV clusters. However, a test is needed to assure that tip calibrations will lead to meaningful estimates. For this purpose, a date-randomisation procedure has been proposed that tests whether the real rate estimate deviates from rates obtained in the absence of temporal structure in the tip-calibrations (Ramsden ). Here, we provide a convenient BEAST implementation of this test that does not require multiple independent randomisations, but averages over all plausible randomisations during the rate estimation process. Based on a relatively stringent test criterion (Duchêne ), we still find a significant association between substitutions and time in our RYMV sequence data.

Although significant, the relative weakness of the temporal signal explains the sensitivity to the coalescent prior we observed in our analyses. Our TMRCA estimate using the skygrid model was somewhat more recent than that obtained under the skyline and skyride model, but more importantly, different settings in the skyride model strongly impacted the evolutionary time-scale. This provides a general warning that weak temporal signal may not simply be reflected in uncertainty of date and/or rate estimation in a Bayesian coalescent framework, but coalescent priors may also affect mean TMRCA estimates. In our study, we highlight the skygrid TMRCA estimate of 1852 [1791-1903] because this model has been shown to outperform the other flexible coalescent priors for divergence time estimation (Gill ). We acknowledge that there is a limit to TMRCA estimation for rapidly evolving viruses in general because saturation and strong purifying selection can lead to an underestimation of old viral origins (Wertheim and Kosakovsky Pond, 2011). However, a TMRCA of about 160 years may be well below the time-scale on which saturation becomes truly important (Bielejec ).

Our study represents a natural extension of earlier descriptions of RYMV spatial genetics (e.g. Abubakar ). Using statistical reconstructions of discrete phylogeographic diffusion, we confirm a clear east-west separation, a strong spatial structure in general, and a likely origin in Tanzania. In this region, and maybe in the Eastern Arc Mountains biodiversity hotspot in particular, RYMV may have emerged in cultivated rice from an ancestor infecting wild graminaceous (e.g. perennial wild rice species such as O. longistaminata) before spreading to other parts of Africa. However, currently identified RYMV isolates in perennial wild hosts appear to be spill-over events from cultivated rice (Traore ), and no related viruses have been identified remote from rice crop areas. In line with an estimated origin in Tanzania, recent analyses of complete genome data have shown that the West and Central African RYMV diversity is nested as a monophyletic clade within the Tanzania diversity (Ochola ). Although we did not consider this in our analysis, we note that a Tanzanian origin also appears to be supported by three coding insertion-deletion polymorphisms at two positions of the CP gene (amino acid position 18 and 60). The three forms are distributed over different clades, but they are only found together in Eastern Tanzania. In the remaining parts of East Africa and in West Africa, only one of the three forms has been identified, with S5 strain (Eastern Arc Mountain in Eastern Tanzania) with both K19 and R60, S6 strain (Eastern Arc Mountain in Eastern Tanzania) with R60 and a deletion of codon K19, and strains S1, S2, S3, Sa, Sg (West Africa), and S4 (Eastern Arc Mountain in Eastern Tanzania) have K19 and a deletion of codon R60. Following its emergence in East Africa, we can only speculate on how the virus was introduced relatively early in West Africa. The long-distance movement event appears to be unique to the natural history of RYMV and could have occurred through human trading practices (Carpenter, 1978).

Phylogeographic analyses allow the description of the spatiotemporal patterns of viral spread, which in turn may lead to the formulation of hypotheses about the underlying processes that shape the dynamics of spread. Agricultural intensification and extensification are considered to strongly facilitate the establishment and epidemic spread of emerging viruses (Thresh, 1982; Elena, 2011), and this has also been invoked as a potential driver of RYMV expansion by molecular epidemiology (Konaté and Fargette, 2001; Abubakar ). Specifically, the increasing adoption of new production modes such as water-fed rice farming, annual double cropping and high-yielding Asian varieties highly susceptible to RYMV are likely to have contributed to its spread (Fargette ; Konaté and Fargette, 2001). By identifying the well supported rates of diffusion, we delineated the major RYMV pathways of spread, but until very recently it has remained challenging to formally test the drivers of spatial spread. We address this here using a GLM extension of the discrete phylogeographic model (Faria ; Lemey ), which aims at determining which subset of explanatory variables helps to explain the relative intensities of viral dispersal among pairs of locations. Agricultural intensification leads to higher rice densities and harvest, which we incorporated as a predictor in our analyses. To this purpose, we used circuit theory to build a resistance landscape with the harvested area of rice in 2000 as resistance factor. This strongly predicted the patterns of RYMV spread, and both the geographic and host ecology component of the distances measured on the resistance landscape appeared to be important. The inclusion of origin rice intensity in 1990 as an additional predictor points at a degree of asymmetry in spread facilitated by rice connectivity, with stronger effect on the dispersal out of areas with high rice intensity, which seems in line with spread facilitated by agricultural extensification. The fact that a more recent measure of rice intensity (1990 vs. 1960) provides better explanatory power may be related to the higher branch density in the recent evolutionary history covering more dispersal events around that time. To our knowledge, our study is the first to formally demonstrate the role of host ecology in plant virus spread using genetic data. It is interesting to note that a historical approach based on a historical map of rice distribution Portères (1950, 1957, 1962) (Pinel-Galzi ), and our statistical approach using a spatiotemporal reconstruction incorporating present day rice statistics, converged towards the same conclusion that harvested rice intensity or connectivity is the main determinant of RYMV emergence and spread.

Earlier RYMV phylogeographic analyses have established an isolation-by-distance pattern for RYMV (Abubakar ), and spread as a function of geographic distance was also evident from our phylogeographic test approach. This motivated the complementary application of phylogeographic reconstruction in continuous space (Lemey ) allowing us to quantify the tempo and mode of spread using several spatial summary statistics (Pybus ). Given the clear distinction of East and West African RYMV lineages, we compared separate estimates from both regions. Considerable differences exist in terms of climate, ecology and host range between these regions, with East Africa growing only the Asiatic rice O. sativa whereas both the African rice O. glaberrima, which is genetically quite different from O. sativa, and the Asiatic rice are cultivated in West Africa. Despite such differences, we found that the overall rate of RYMV spread and diffusivity was highly similar. These measures do however not take into account the directionality of spread, and when the directionality is considered to be the distance from the estimated origin in both regions (the wavefront distance), we find higher wavefront velocities in the West. So, it seems that rice densities, which are more pronounced in the West, do not necessarily increase the overall rate of spread, but they facilitate more extensive expansion dynamics. By summarising these expansion dynamics over time (Fig. 7), we revealed that agricultural intensification and extensification had a more prominent impact in the West, which is in accordance with the fact that our GLM-diffusion estimates were essentially informed by the data from West Africa. The major expansion dynamics in the west are characterised by spread from Côte d’Ivoire or Mali in an eastern direction, towards countries that have relatively lower rice production like Niger, Chad, or Central African Republic, in line with the asymmetry suggested by the origin rice intensity predictor in the GLM analysis. In Mali, various lineages have been identified in the Inner Niger Delta, suggesting that this specific area may have been the West African centre of diversification (Traore ; Fargette ). The propagation towards Central Africa is likely to have followed the more accessible routes of transmission, in particular along the Niger-Bénoué rivers. In East Africa, the comparatively sparser and less intense rice production, along with physical (mountains) and ecological (tropical forests) boundaries have restricted viral expansion. The continuous phylogeographic reconstruction shows only recent spread out of its likely area of origin. However, further spread may continue in the future and molecular surveillance will be needed to track these dynamics.

Several aspects of our phylogeographic analyses may be further improved or fine-tuned in the future. Longer genome regions offer more phylogenetic resolution and are likely to increase the temporal signal, but they may also require taking into account recombination (Pinel-Galzi ). Sampling biases represent an important challenge for ancestral reconstructions, and we acknowledge that such biases may also burden the sample we analysed, even though the GLM analysis did not associate sampling numbers with diffusion intensities. Structured coalescent approaches are expected to be less sensitive to sampling biases and represent interesting alternatives for discrete phylogeographic reconstructions (De Maio ). Furthermore, it may prove interesting to expand on the predictors of RYMV dispersal if systematic data would be available, including for example, on vector demographics, rice cultivar resistance to RYMV, mode of watering and other agricultural practices. Despite these areas of potential improvement, our current analysis takes an important step towards hypothesis testing in plant virus epidemiology and ecology. Over the last decade, RYMV has become the main threat to rice cultivation in Africa and Madagascar (Konaté and Fargette, 2001). The finding that host ecology has shaped RYMV spread suggests predictable patterns of spread that may help to inform predictive models for RYMV control and public health policies. More generally, our results reinforce the concept that host population ecology is crucial for the onward transmission and epidemic potential of any emerging virus (Woolhouse and Gowtage-Sequeria, 2005).

Supplementary data

Supplementary data are available at Virus Evolution online.