{'fraction_non_alphanumeric': 0.03138444567015996, 'fraction_numerical': 0.023496966354109212, 'mean_word_length': 5.106601549343213, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Identifying the patterns and drivers of Puumala hantavirus enzootic dynamics using reservoir sampling

Abstract Hantaviruses are zoonotic hemorrhagic fever viruses for which prevention of human spillover remains the first priority in disease management. Tailored intervention measures require an understanding of the drivers of enzootic dynamics, commonly inferred from distorted human incidence data. Here, we use longitudinal sampling of approximately three decades of Puumala orthohantavirus (PUUV) evolution in isolated reservoir populations to estimate PUUV evolutionary rates, and apply these to study the impact of environmental factors on viral spread. We find that PUUV accumulates genetic changes at a rate of ∼10−4 substitutions per site per year and that land cover type defines the dispersal dynamics of PUUV, with forests facilitating and croplands impeding virus spread. By providing reliable short-term PUUV evolutionary rate estimates, this work facilitates the evaluation of spatial risk heterogeneity starting from timed phylogeographic reconstructions based on virus sampling in its animal reservoir, thereby side-stepping the need for difficult-to-collect human disease incidence data.


Introduction
Hantaviruses are important zoonotic pathogens that can cause hemorrhagic fever with renal syndrome (HFRS) and hantavirus pulmonary syndrome, potentially life-threatening diseases in humans (Maes et al. 2004). Concomitant with the increasing public health relevance of hantaviruses (Kruger et al. 2015), there is a growing interest in uncovering the present-day patterns and drivers of the spatiotemporal dynamics of hantavirus circulation in its host reservoirs (Guivier et al. 2011;Barrios et al. 2013a,b;Haredasht et al. 2013a,b;Korva et al. 2013;Weber de Melo et al. 2015;Laenen et al. 2016;Tian et al. 2017) due to its potential in assisting risk determination and targeting surveillance efforts. Human infections result from inhaling infectious aerosols from excreta of infected animals. HFRS incidence is strongly correlated with the reservoir host population dynamics, resulting in seasonal and multi-annual fluctuations in incidence rates and sporadic local outbreaks (Kallio et al. 2009;Haredasht et al. 2014;Luis et al. 2015). Urbanization and landscape factors such as habitat composition, size and quality affect small mammal biodiversity and shape host distribution and metapopulation dynamics, consequently influencing hantavirus transmission within the host population and spillover to humans (Guivier et al. 2011;Ecke et al. 2017;Tian et al. 2018). Because of the difficulties associated with collecting reliable incidence data covering longer time periods and large geographic areas (e.g. due to subclinical infections (Clement et al. 1997) and differences in notification status of disease among regions (Heyman, Vaheri, and Members 2008)), the imprint of epidemic events that is contained in virus genetic data constitutes an attractive alternative source of information (see Biek and Real 2010;Bloomquist, Lemey, and Suchard 2010;Pybus, Tatem, and Lemey 2015;Holmes et al. 2016; for recent reviews on this topic). To use such imprints in genomic sequences for identifying the drivers of virus emergence and spread; however, it is crucial that the epidemiological history of the virus of interest that underlies all inferences, is accurately dated. This can be achieved by exploiting the temporal dimension in samples collected over the last few years or decades. Such temporally spaced virus genomic data can reliably inform evolutionary rate estimates when they foster a statistically significant number of genetic differences (Drummond, Pybus, and Rambaut 2003), and it is now well-recognized that this assumption should be thoroughly tested in order to avoid spurious results (Firth et al. 2010;Duchê ne et al. 2015;Murray et al. 2016;Rambaut et al. 2016).
In Western and Northern Europe, Puumala orthohantavirus (PUUV) is the principal causative agent of hantavirus disease in humans (Maes et al. 2004;Kruger et al. 2015). The natural host of PUUV is the bank vole (Myodes glareolus), which has broadleaf forests with dense undergrowth as a preferred habitat in Western Europe (Mazurkiewicz 1994;Linard et al. 2007). Bank voles remain persistently infected with PUUV, but there can be remarkable dissimilarity in the persistence of PUUV variants over time within local vole populations, with both high rates of lineage turnover and preservation of identical variants for several years (Johansson et al. 2008;Weber de Melo et al. 2015). This extensive variation in the substitution rate among lineages likely results from a complex interplay between virus and host ecology. For example, Voutilainen et al. (2015) reported a strong variability in individual PUUV excretion dynamics of bank voles in natural settings, possibly reflecting differences in virus replication due to fluctuating immune responses as a result of changes in resource availability. Bank voles can get infected with PUUV either through direct contact or indirectly via contaminated food or environment (Sauvage et al. 2003). The relative importance of each transmission path depends on seasonal fluctuations in ecological factors and bank vole behavior (Tersago, Verhagen, and Leirs 2011) and on multiannual population density cycles (Voutilainen et al. 2016). Additionally, hantavirus ex vivo stability varies substantially depending on the environmental conditions (Kallio et al. 2006;Linard et al. 2007), thereby influencing the efficiency of the indirect transmission path. Thus, population-specific differences in transmission dynamics can affect the PUUV evolutionary rate as well as its constancy. This is important because substantial deviations of the linear accumulation of divergence through time can confound evolutionary rate estimations (Trovao et al. 2015;Raghwani et al. 2016). Furthermore, within-host life-history traits that counter the effects of short-sighted evolution may add to the high heterogeneity in evolutionary rates along different lineages, as is often observed in chronic viral infections (Vrancken et al. 2014;Lythgoe et al. 2017;Vrancken, Suchard, and Lemey 2017).
We here capitalize on the longest available sampling through time of well-defined PUUV populations infecting bank voles to estimate the rate of PUUV evolution at a timescale relevant for the present-day epidemic emergence and spread. We first exhaustively investigate the veracity of the temporal signal before estimating substitution rates and, in the timestructured clades, evaluate the impact of several landscape factors on the rate of PUUV spread.

Dataset compilation
A total of 664 bank voles were trapped at eight locations in Belgium between 2006 and 2016, and ninety-eight (14.8%) were PUUV IgG positive. Complete S segment sequences were generated from eighty-two (83.7%) of these samples. In addition, 2 S segments from early passage PUUV isolates originating from bank voles captured between 1984 and 1986, were also sequenced. The new sequence data were complemented with the available PUUV S-segment data from NCBI GenBank.

Identifying isolated populations
We first focused on the fragment of the S segment that is most frequently sequenced (nucleotide (nt) positions 400-900 in PUUV strain Sotkamo, NC_005224) to frame the Belgian PUUV diversity in a global perspective. An unconstrained analysis (with no molecular clock enforced) using MrBayes showed that all Belgian PUUV lineages group in a well-supported clade that also includes strains from nearby locations in France (Charleville (Castel et al. 2015)), Germany (Cologne (Essbauer et al. 2007), Erft (Heiske et al. 1999), Koblenz (Schilling et al. 2007)), and The Netherlands (Middelbeers, Oirschot, and Kaatsheuvel (de Vries et al. 2016)) Supplementary Fig. S1. Next, we focused on this 'Belgian' cluster and extended the scope to the entire coding region of the S segment (1,299 nt) to more reliably identify PUUVs that have been evolving in isolated bank vole populations. This analysis showed that the Belgian lineages group in three well-defined clades (each having 100% posterior support), that can roughly be assigned to the north, middle, and south of Belgium (Fig. 1). These lineages will in this manuscript be further referred to as the Campine, Sonian Forest, and Ardennes clades, respectively.

Validation of temporal structure and estimating evolutionary rates
To explore the presence of temporal structure and the data quality for the Campine, Sonian Forest, and Ardennes clades ( Supplementary Fig. S2) the genetic divergence was regressed against time. This revealed a complete lack of temporal structure in the Sonian Forest clade ( Supplementary Fig. S2C), that only has 4 months between the earliest and the most recent available sequence. In contrast, there appeared to be a timerelated accumulation of mutations in the Campine and Ardennes clades, and for this reason, we focused on these PUUV populations in subsequent explorations (Supplementary Fig. S2A and B). For neither clade, there was evidence of confounding between the genetic and temporal structure (Table 1), and the regression-based analyses indicate that the observed correlation between the rootto-tip distances and time only differs significantly from what is expected in the absence of a temporal structure for the Ardennes lineage (Murray et al. 2016). Additionally, the proportion of explained variation (R 2 ) is low for both the Ardennes and Campine datasets (Table 1), and large residuals indicate that the data may better fit a model where the substitution rate can vary along the branches (Supplementary Fig. S2A and B).
The suboptimal R 2 values of the regression analysis and the possible confounding impact of the non-independence of the data (Drummond, Pybus, and Rambaut 2003;Vrancken, Suchard, and Lemey 2017), motivated us to further evaluate the presence of temporal signal using a tip-date randomization test in conjunction with BEAST (Firth et al. 2010). In addition to appropriately accounting for the shared ancestry, this also allows to capture differences in the rate of evolution along lineages with a relaxed molecular clock (Drummond et al. 2006). Relaxed clock models provide an estimate of the variation in evolutionary rate among lineages, the coefficient of variation (CoV) statistic that can be used to evaluate the clocklikeness of evolution. Specifically, the CoV statistic depicts the scaled variance in evolutionary rate among lineages, meaning that a strict clock model may fit the data better than a relaxed clock model when   clades was estimated using a relaxed clock model with branch rates drawn from an underlying discretized log-normal distribution (Drummond et al. 2006). For the Campine clade, the CoV was also assessed after imposing a monophyly constraint on the Gierle-Zevendonk-Turnhout subclade. Correspondence between the posterior densities and the datasets is indicated in this figure.
the 95 per cent highest posterior density (HPD) interval of the CoV statistic borders zero (Drummond et al. 2006;Gray et al. 2011). Reconstructions of the demographic histories of the Ardennes and Campine clade using a non-parametric coalescent model were in agreement with a constant population size scenario for both clades (Fig. 2). This was implemented in subsequent phylogenetic inferences. The credible interval of the CoV for the Ardennes dataset ( Fig. 3) buts to zero, but the rather liberal range of plausible CoV values indicates that a relaxed and strict clock model may have a comparable goodness of fit of the data. The Campine CoV posterior density (Fig. 3) also shows that the CoV cannot be used to unambiguously set apart different clock model hypotheses, albeit the support for rate heterogeneity is more outspoken. Because of these ambiguities, we formally compared the fit of strict and relaxed clock models to both datasets using a Generalized Stepping Stone marginal likelihood estimator in BEAST (Table 2) (Baele, Lemey, and Suchard 2016). In agreement with the patterns that emerge from Fig. 3, we find that both types of clock models fit the Ardennes data comparably well and the relaxed clock model fits the Campine data slightly better.
Additionally, the Campine CoV posterior density indicates that levels of rate variation comparable to those observed for the Ardennes clade are slightly better supported than higher levels (Fig. 3). Furthermore, a comparison of the maximum clade credibility (MCC) trees, corresponding to the sampled states of the converged Markov Chain Monte Carlo (MCMC) chain with a CoV larger respectively smaller than 1, shows that extreme levels of among branch rate variation (CoV > 1) involve unrealistically high rates ($10 À2 substitutions/site/year) on particular internal branches linking the subclades ( Supplementary  Fig. S3). This is likely consequential to the combination of the parameter-rich relaxed clock model and data with noisy and/or weak temporal signal. Because very non-clock like data impede reliably estimating substitution rates (Trovao et al. 2015;Raghwani et al. 2016), we attempted to shrink the level of among-lineage rate heterogeneity in the Campine clade by imposing a monophyly constraint that involves the Gierle-Zevendonk-Turnhout subclade (see Supplementary Fig. S3, red box). The rationale is that this subclade has 100 per cent posterior support for being monophyletic in the unconstrained analysis and requiring these taxa to co-cluster does not allow for the topology of the MCC tree that is seen for high levels of among branch rate variation (CoV > 1, Supplementary Fig. S3).
Because of these uncertainties, the tip-date randomization test was done with both types of clock models (Fig. 4). This reveals that the Ardennes PUUVs clearly represent a 'measurably evolving population' (Drummond, Pybus, and Rambaut 2003). The evolutionary rate estimates for the Ardennes clade are 1.8 Â 10 À4 s/s/y (95% HPD: 1.1 Â 10 À4 À 2.5 Â 10 À4 ) under a strict clock and 2.0 Â 10 À4 s/s/y (95% HPD: 1.2 Â 10 À4 À 3.0 Â 10 À4 ) under a relaxed clock model and are robust to the clock model specification (Figs 4 and 5). Although the monophyly constraint on the Gierle-Zevendonk-Turnhout subclade improved the model fit (Table 2) and indeed reduces the level of rate heterogeneity (Fig. 3), it is clear from the large overlap of the 95 per cent Bayesian credible interval (BCI) of the tip-date informed evolutionary rate estimate with the 95 per cent BCIs of the estimates obtained from tip-date randomized datasets, that the substitution rate estimate from the Campine dataset is not as clearly informed by the data as for the Ardennes clade (Fig. 4). Yet, that the posterior density of the tip-date informed substitution rate is shifted towards higher values when compared with the null distribution (i.e. the expectation when randomizing the sampling dates), shows that there is signal for a build-up of genetic divergences over the sampling time. Indeed, analogously to the decreasing slope (i.e. substitution rate) with decreasing clock signal in a regression of root-to-tip divergence versus sampling time, it is expected that with decreasing temporal signal the substitution rate will increasingly tend towards the null distribution, which is opposite to what is observed here. Furthermore, the rate estimates of the Campine dataset largely overlap with the Ardennes estimates, irrespective of the type of clock model (Fig. 5). Overall, these results show that the temporally annotated PUUV genetic data can be used to reliably frame this virus' recent evolution in absolute time.

The impact of landscape factors on the Puumala virus dispersal velocity
Combining geographic data with dated sequence evolution enables the reconstruction of virus dispersal dynamics. In turn, this allows relating molecular epidemiology with landscape ecology and testing of hypotheses about the effect of landscape configuration on PUUV dispersal. Here, we pursue this by first reconstructing the dispersal history of the Ardennes and Campine clades using a Bayesian continuous phylogeographic method . In a next step, the SERAPHIM R package (Dellicour et al. 2016;Dellicour, Rose and Pybus 2016) is used to investigate the impact of three landscape features (forest areas, croplands, and primary roads) on the PUUV dispersal velocity.
From the land cover maps with superimposed geo-referenced phylogenies (Fig. 6), it is clear that PUUV migration involves both forested and non-forested (mostly croplands in these regions) areas, and that primary roads form no absolute barriers. When relating landscape features to PUUV virus dispersal velocity, we find modest Bayes factor support (BF > 3) for forest areas facilitating PUUV dispersal in the Ardennes clade, while croplands act as a resistance factor, impeding PUUV dispersal (Fig. 7). For the Campine region, however, no significant support was obtained for any of the tested landscape features. Furthermore, primary roads do not act as a significant resistance factor for PUUV dispersal velocities in either clade (Supplementary Tables S2 and S3). However, the level of statistical support does not quantify the contribution of landscape factors in explaining differences in migration rates. For this, we consider the ratio of explained heterogeneity of the environmental versus null model (R 2 env /R 2 null , Supplementary Tables S2 and S3), that shows that forest and non-forest cover have a similar added value, with an average increase of R 2 of 42.2 and 43.8 per cent above what is captured by geographic distances alone. Importantly, the high absolute values of R 2 env when treating cropland as a resistance (median 0.59, 95% CI: 0.19-0.80) and forest as a conductance factor (median 0.55, 95% CI: 0.16-0.80) testify to the importance of both types of land cover in shaping PUUV dispersal.

Discussion
Collecting genomic data is becoming a routine part of many hantavirus surveillance programs. Through phylogenetic analyses, this enables linking virus evolution to migration processes as well as statistically identifying correlates of epidemic spread, without having to rely on time series of incidence rates. These incidence data are often absent or scarce and can be unpredictable owing to the difficulties in their collection. To reliably bridge evolutionary information with data on epidemiological processes, however, it is crucial that the statistical models used to connect genetic information to time-stamped data-referred to as 'molecular clock models'-are accurately calibrated.
Based on a unique sampling effort covering approximately three decades, we obtained the best-informed short-term  evolutionary rate estimates of the PUUV nucleocapsid gene (S segment) to date, and exhaustively tested their reliability. Our work provides a substantial improvement on previously published PUUV substitution rate estimates in three related aspects. First, by considering PUUV evolution in isolated populations, as opposed to mixing data from populations that evolved in geographic isolation over extended periods of time (Ramsden et al. 2008), we avoid the biasing effect of the time-dependent rate decline that is expected to manifest on the long, deep branches connecting taxa from separated populations (Bennett et al. 2014 Tables S2 and S3 for the complete results). BF supports were estimated with a randomization procedure detailed in Supplementary Appendix S1 cover a long range of the variable under investigation. Third, the authenticity of the temporal signal was comprehensively analyzed, which is a prerequisite for reliance on the reported estimates (Zhang and Holmes 2014). We demonstrate that the rate estimate for the Ardennes population is well-informed by the data. Although the temporal structure in the Campine clade is less pronounced, the estimates for both populations are in close agreement (Fig. 5). Strikingly, our PUUV evolutionary rate estimates are also highly similar to substitution rates reported for S segment in other hantaviruses estimated from datasets with comparable genomic coverage and sampling time frame (Supplementary Fig. S4). This indicates that, contrary to what is often observed for other viruses (Streicker et al. 2012;Worobey, Han, and Rambaut 2014), differences in life-history traits appear to not substantially impact short-term evolutionary rates of these, and perhaps all, hantaviruses. The high transition/transversion ratios observed for the Ardennes and Campine datasets (with respective mean ratios of twelve (95% BCI: 8-16) and seventeen (95% BCI: 9-26)) testify that synonymous mutations are favored over amino acid changing mutations. Presumably, this is due to a history of virus-host co-evolution of these hantaviruses being sufficiently long to have reached a fitness plateau, resulting in strong purifying selection, as was previously reported for PUUV and other hantaviruses (Razzauti et al. 2013;Laenen et al. 2016).
Bank vole population sizes and PUUV prevalence can be highly variable (Tersago et al. 2008;Weber de Melo et al. 2015;Voutilainen et al. 2016). This, however, is not reflected in changes in the effective number of infections through time (Fig. 2), what indicates that expanding lineages are usually transient and have remained unsampled. As a matter of fact, it has been shown that the temporal distribution of sampling determines to what extent fluctuations in cyclic growth patterns are reflected in coalescent-based reconstructions of the demographic history (Stack et al. 2010). This illustrates the need for surveillance programs with an improved spacing of sampling times, preferentially generating long sequences from all three hantavirus segments to capture the multiannual population dynamics in detail. Improved sampling will also allow relating the periodic virus population growth and decline to the HFRS incidence and its potential drivers, such as host population density and climatologic factors, while avoiding the problems associated with post hoc methods (Gill et al. 2016;Volz and Didelot 2018). In addition, enhanced sampling is needed to assess the importance of segment reassortment for PUUV emergence dynamics (Razzauti et al. 2013;Szabo et al. 2017;Laenen et al. 2018).
Persistence of a limited set of lineages can also account for the absence of a clear imprint of variation in the rate of evolution between lineages in the virus genome. Specifically, many evolutionary dead-ends can eliminate the diversifying impact of lineage-specific differences in transmission dynamics on the substitution rate when, of the many variants that evolve at different rates at a single moment, only those with similar transmission dynamics survive in the longer-term. Importantly, both this and the constancy in population size (Fig. 2) align with the recent observation of long-term shedding of PUUV by bank voles in natural settings as a means to overcome host population bottlenecks (Voutilainen et al. 2015) and support the view that this hantavirus' long-term life strategy includes prolonged periods of replication in chronically infected hosts.
The PUUV and, more generally, the hantavirus emergence potential is sensitive to changes in environmental factors, and many modeling efforts have attempted to identify climatic, anthropogenic, and other ecological determinants of hantavirus-associated disease incidence (Linard et al. 2007;Palo 2009;Barrios et al. 2010Barrios et al. , 2013aTersago et al. 2011;Haredasht et al. 2013a,b;Khalil et al. 2014;Voutilainen et al. 2016;Tian et al. 2017Tian et al. , 2018. To date, conditioning inferences on the spatial and temporal dynamics of gene flow has been unattainable for these studies. Thanks to the presence of sufficient temporal signal in our data, however, we could exploit the concordant timescales of epidemiological and evolutionary dynamics to evaluate the influence of environmental variables on PUUV dispersal velocity. As forests are the preferred, yet not unique, habitat of bank voles (Gurnell 1985;Canova and Fasola 1991;Van Apeldoorn et al. 1992;Mazurkiewicz 1994;Wijnhoven et al. 2005;Torre and Arrizabalaga 2008), and roads can act as barriers to bank vole migration (Gerlach and Musolf 2000), we opted to investigate the relevance of forest and non-forest cover, and the presence of main roads in explaining migration rate heterogeneity. There was positive support for forest cover to facilitate PUUV dispersal and, vice versa, for non-forest cover as a resistance factor in the Ardennes clade. On the other hand, the type of land cover could not explain the variation in migration rates better than geographic distance alone in the Campine population. This is not due to a low amount of dispersal velocity heterogeneity in the Campine clade ( Supplementary Fig. S5), but most likely relates to the substantial fragmentation of land cover in the Campine area, leading to less opportunity for a particular land cover type to leave an imprint on the migration rate. However, it should be noted that the environmental factors are assumed constant in time, and our approach currently cannot consider time-variable predictors. This may be particularly relevant for the Campine region, which has been subjected to substantial deforestation in the 20th century (Van der Veken, Verheyen, and Hermy 2004).
Given the impact of landscape structures, in particular forest fragmentation, on host abundance and spatial behavior (Linard et al. 2007;Guivier et al. 2011), it is likely that the extent of land cover fragmentation in the Campine area accounts for the noticeable lower effective number of infections through time (mean = 43, 95% BCI: 18-72) when compared with the Ardennes (mean = 100, 95% BCI: 58-153), where PUUV can more rapidly be reseeded following local extinction (Guivier et al. 2011). The relevance of land cover is further highlighted by historical HFRS incidence data: in the north of Belgium (cf. Campine clade), where forests are more fragmented, a lower incidence of HFRS is seen, while for the south of Belgium (cf. Ardennes clade) where forests are more connected, an endemic region with high HFRS incidence rate is present. These results further strengthen the notion that the type of land cover can be a valuable variable in predictive models of PUUV spread and infection risk (Linard et al. 2007;Guivier et al. 2011;Barrios et al. 2013b). Finally, in line with previous findings that only highways strongly impede bank vole migration (Gerlach and Musolf 2000), the presence of primary roads does not seem to affect the migration rate of PUUV in the Ardennes and Campine regions.
In summary, we demonstrated that PUUV genomic surveillance data are amenable to quantifying the rate of the PUUV dispersal dynamics and identifying its drivers. Our results will assist in identifying intervention points and improve the efficiency of public health resource allocations.

Sample collection and sequencing
Rodent trappings were conducted on eight different locations in Belgium from 2006 to 2016. On each location, one hundred Longworth or Sherman live traps were placed in a 10 Â 10 grid at 10 m intervals at each location. Captured bank voles were sedated with isoflurane followed by a 10 ml blood collection through the retro-orbital sinus. The collected blood was tested in the field for the presence of PUUV IgG antibodies using the ReaScan V R Ab-Dect Puumala IgG kit (Reagena). PUUV antibodypositive bank voles were sacrificed and immediately stored at À20 C. Bank vole kidney and lung tissues were dissected and preserved in RNAlater stabilization solution (Life Technologies). Kidney and/or lung tissue was homogenized with the Minilys homogenizer (Bertin Technologies) followed by total RNA extraction with the RNeasy Mini kit (Qiagen), according to the manufacturer's instructions. Tissues were screened for the presence of hantavirus RNA as described previously (Laenen et al. 2016). In addition, PUUV isolates, obtained from low passaging on VeroE6 cells, from bank voles captured between 1984 and 1986 were included in sequencing. From the positive samples and cultures, the complete S segment was amplified with the One-Step RT-PCR kit (Qiagen) using the following conditions: 50 C for 30 min, 95 C for 15 min followed by 40 cycles of amplification (30 s at 94 C, 30 s at 53 C, 2 min at 72 C) with a final extension of 10 min at 72 C. PCR amplicons were purified using ExoSAP-IT V R PCR Product Cleanup (Affymetrix) and sequenced according to the ddNTP chain termination method with the BigDye Terminator v3.1 cycle sequencing kit (Life Technologies) on an Applied Biosystems 3130xl Genetic Analyzer. Primers used to amplify and sequence the S segment are available in Supplementary Table S1.

Querying the GenBank sequence database
All available PUUV sequence data were downloaded from NCBI GenBank. S segment sequences were filtered by aligning all data to a reference S segment sequence (GenBank accession number NC_005224) with a pairwise codon-aware alignment tool (http:// regatools.med.kuleuven.be/sequencetool/sequencetool.wt), allowing for at most ten frameshifts.

Bayesian inference of sequence evolution
For the analyses where no molecular clock was enforced, MrBayes v.3.2.6 (Ronquist et al. 2012) was used to estimate the posterior probabilities of trees using Bayesian inference through MCMC simulation, employing a General Time Reversible (GTR) substitution model and a discretized C distribution using default priors (Yang 1994).
For the specific clades where preliminary tests (see below) indicated the presence of a measurable build-up of genetic change over the sampling time span (the Ardennes and Campine clades), timed evolutionary histories were estimated using BEAST v1.8.4 (Drummond et al. 2012). Here, a strict or relaxed molecular clock model (Drummond et al. 2006) was combined with the HKY substitution model (Hasegawa, Kishino, and Yano 1985) and substitution rate heterogeneity across positions was again modeled with a discretized C distribution (Yang 1994). After an initial exploration of the demographic history with a flexible non-parametric coalescent model (Gill et al. 2013; see Fig. 2), we decided to assume a constant population size for both datasets when estimating clade-specific substitution rates. The fit of different models to the data was estimated using the Generalized Stepping Stone marginal likelihood estimator as implemented in BEAST (Baele et al. 2016).

Evaluating the temporal signal
The initial evaluation of the temporal signal and data quality with TempEst (Rambaut et al. 2016) was based on a maximum likelihood tree estimated with PhyML in Seaview (Gouy, Guindon, and Gascuel 2010), with GTR correction for multiple hits, C-distributed rate variation among sites (Yang 1994) and not assuming a molecular clock. A Mantel test (Murray et al. 2016) served to identify whether genetically similar taxa were more likely to have been sampled around the same time, and the significance of the correlation between the root-to-tip divergences and sampling times was estimated against a null distribution obtained by 1,000 times randomly reassigning the sampling dates to the taxa. The significance of the tip-date informed evolutionary rate signal was also determined, while accounting for the genetic relationships. This involved creating twenty datasets with randomly permuted tip-dates (Murray et al. 2016) and comparing the mean substitution rate estimate from the observed data with the 95 per cent BCIs estimated from the randomized datasets (Firth et al. 2010), assuming the same evolutionary models as in the tip-date informed analyses.

Testing the impact of landscape factors on the rate of geographic spread
The workflow for statistically analyzing the impact of landscape features (forest areas, croplands, and primary roads) on the rate of PUUV geographic spread is similar to those previously described by Dellicour, Rose and Pybus (2016) and Dellicour et al. (2017) and is summarized below. All required scripts are available in the SERAPHIM R package (Dellicour et al. 2016;R Core Team 2017). A more detailed exposition of this method is available in Supplementary Appendix S1.
Spatially and temporally referenced phylogenies for the Ardennes and Campine clades were inferred with the continuous phylogeographic method implemented in BEAST 1.8.4 ), using the same demographic and substitution models as those used for estimating the substitution rates were specified (see above). Because the evolutionary rate estimates largely overlap irrespective of the clock model (Fig. 5) and a strict clock model is not rejected by the data (Table 2), we chose the strict clock model for an optimal dating precision (Ho et al. 2005;Drummond et al. 2006). Information on the timed history of spread was extracted from a subset of 100 trees sampled at regular intervals from the post burn-in posterior distribution of trees. Specifically, all branches are treated as movement vectors that represent conditionally independent viral lineage dispersal events (Pybus et al. 2012) with known start and end locations (latitude and longitude) and start and end dates. These vectors are subsequently used to investigate the impact of landscape features on the PUUV dispersal velocity by computing 'environmental distances' that is spatial distances that are weighted according to the values of an environmental raster mirroring the spatial heterogeneity of the tested landscape features (Dellicour, Rose and Pybus 2016). Environmental distance computation was performed using two different models: 1, the least-cost path model, which uses a least-cost algorithm to obtain the migration route between the start and end points (Dijkstra 1959), and 2, the Circuitscape path model that can accommodate uncertainty in the migration route (McRae 2006;Mcrae et al. 2008). By correlating environmental distances to the time that was needed by the virus to migrate from start to end location, one can test the impact of landscape features on the dispersal velocity. Practically, this is done by calculating the statistic R = R 2 env /R 2 null , where R 2 env is the coefficient of determination obtained by regressing branch durations against environmental distances computed on the environmental raster, and R 2 null is the coefficient of determination obtained by regressing branch durations against environmental distances computed on the 'null' raster. The R statistic, therefore, represents how much variation in migration rates, beyond that explained by geographic distance alone, is explained by considering the spatial heterogeneity of the environmental variable. Note that this statistic R is slightly different from the statistic Q (= R 2 env À R 2 null ), previously used in Dellicour et al. (2017).
Because it is a priori unknown to what extent the geographic variation of the environmental factor affects the level at which a cell in the raster acts as a conductance/resistance factor relative to the null raster, the environmental distances were calculated from several rasters with rescaled cell values (Dellicour et al. 2017). For this, the original raster cell values were transformed according to v t = 1 þ k*(v o /v max ), where v t and v o are the transformed and original raster cell values, v max is the maximum raster cell value recorded in the raster and k is 10, 100, or 1,000. Also, each landscape feature was tested as a resistance factor (i.e. impeding virus movement) and as a conductance factor (i.e. facilitating virus movement). In summary, for each combination of environmental factor, k parameter value and path model, an R distribution is derived from the 100 sampled phylogenies.
A landscape factor was considered as potentially explanatory, and BF support for the significance of an R distribution was calculated only when the distribution of regression coefficients was positive and when p (R > 1) was at least 90 per cent.

Data availability
Sequences generated for this study have been submitted to NCBI GenBank under accession numbers MG812385-MG812468.

Supplementary data
Supplementary data are available at Virus Evolution online.