# Emerging strains of watermelon mosaic virus in Southeastern France: model-based estimation of the dates and places of introduction

Lionel Roques, Cecile Desbiez, Karine Berthier, Samuel Soubeyrand, Emily Walker, Etienne K. Klein, Jimmy Garnier, Benoît Moury, Julien Papaïx

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Lionel Roques, Cecile Desbiez, Karine Berthier, Samuel Soubeyrand, Emily Walker, et al.. Emerging strains of watermelon mosaic virus in Southeastern France: model-based estimation of the dates and places of introduction. Scientific Reports, 2021, 11, pp.7058. 10.1038/s41598-021-86314-y. hal-03197633

HAL Id: hal-03197633

https://hal.inrae.fr/hal-03197633

Submitted on 14 Apr 2021

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

Check for updates

OPEN

# Emerging strains of watermelon mosaic virus in Southeastern France: model-based estimation of the dates and places of introduction

L. Roques<fnref n="1" />✉, C. Desbiez², K. Berthier², S. Soubeyrand<fnref n="1" />, E. Walker<fnref n="1" />, E. K. Klein<fnref n="1" />, J. Garnier³, B. Moury² & J. Papaïx<fnref n="1" />

Where and when alien organisms are successfully introduced are central questions to elucidate biotic and abiotic conditions favorable to the introduction, establishment and spread of invasive species. We propose a modelling framework to analyze multiple introductions by several invasive genotypes or genetic variants, in competition with a resident population, when observations provide knowledge on the relative proportions of each variant at some dates and places. This framework is based on a mechanistic-statistical model coupling a reaction-diffusion model with a probabilistic observation model. We apply it to a spatio-temporal dataset reporting the relative proportions of five genetic variants of watermelon mosaic virus (WMV, genus *Potyvirus*, family *Potyviridae*) in infections of commercial cucurbit fields. Despite the parsimonious nature of the model, it succeeds in fitting the data well and provides an estimation of the dates and places of successful introduction of each emerging variant as well as a reconstruction of the dynamics of each variant since its introduction.

Plant and animal species—and as a consequence, their pathogens and pests—are translocated across the globe at an ever-increasing rate since the nineteenth century<fnref n="1" />. Most introductions of microorganisms and insects (notably plant pathogens and/or their arthropod vectors) are accidental, while a few are deliberate for biological control purposes². Introductions of plant pathogens can have a dramatic impact on agricultural production and food security³,⁴ and being able to predict the fate of a newly introduced disease based on surveillance data is pivotal to anticipate control actions³.

Where and when alien organisms are successfully introduced are central questions for the study of biologi-cal invasions, to elucidate biotic and abiotic conditions favorable to the introduction and establishment and spread of invasive species⁶,⁷. Indeed, unraveling such conditions is a prerequisite to map the risk of invasion and to design an efficient surveillance strategy⁷ as initial phases are critical for the establishment success of intro-duced pathogens⁸,⁹. It is also during this phase that invaders are the easiest to control⁶,¹⁰. In addition, a precise knowledge of the dates and places of introduction is critical to accurately determine invader reproduction and dispersal parameters, as their estimated values are directly dependent on the time and distance between the actual introduction and the observations¹¹,¹². However, biological invasions by alien organisms are often reported several years after the initial successful introduction event¹³. Thus, monitoring data generally do not provide a clear picture of the date and place of introduction.

Many successful emergences of plant viruses have taken place worldwide in the last decades¹⁴,¹⁵, whereas some have only been observed punctually with no long term establishment of the pathogens¹⁶. As plant viruses are rapidly-evolving pathogens with small genomes and high mutation rates¹⁷, they may present measurably evolving populations¹⁸ over the time scale of decades. Thus, the spatio-temporal histories of invading species can sometimes be reconstructed from georeferenced and dated genomic data, with phylogeographic methodologies¹⁹⁻²¹ or population dynamic models embedding evolutionary processes. Over shorter time-scales, or when only abundance data are available, recent mechanistic modelling approaches have been proposed to infer the date

<footnote n="1">¹INRAE, BioSP, 84914 Avignon, France. ²INRAE, Pathologie Végétale, 84140 Montfavet, France. ³Laboratoire de Mathématiques (LAMA), CNRS and Université de Savoie-Mont Blanc, Chambéry, France. ✉email: lionel.roques@inrae.fr</footnote>
and place of introduction of a single species along with other demographic parameters²². In many modeling approaches, the interactions of the new viruses or new strains with preexisting virus populations are not taken into consideration, even if it is known that synergism or competition between virus species or strains can affect their maintenance and spread²³⁻²⁵.

Mechanistic models are increasingly used in statistical ecology because, compared to purely correlative approaches, their parameters can directly inform on biological processes and life history traits. Among them, the reaction-diffusion framework is widely used in spatial ecology to model dynamic species distributions²⁶,²⁷. In such models, the diffusion coefficient is related to dispersal ability and the growth or competition coefficients may help in understanding the respective interactions between different species, variants or genotypes. Additionally, reaction-diffusion models can easily account for spatio-temporal heterogeneities²⁸,²⁹. A drawback of this type of approach is that the compartments that are modelled typically correspond to continuous population densities, which rarely match with observation data. There is therefore a challenge in connecting the solution of the model with complex data, such as noisy data, binary data, temporally and spatially censored data. Recent approaches have been proposed to bridge the gap between reaction-diffusion models and data (e.g.,³⁰, in a framework known as mechanistic-statistical modelling³¹⁻³³. An advantage of this method is that it allows to estimate simultaneously both the biological parameters and the date and place of introduction in a single framework²². Recently, this type of approach was applied to localize and date the invasion of South Corsica by *Xylella fastidiosa*, based on a single-species reaction-diffusion model and binary data¹¹.

Our objective here is to propose a modelling framework to deal with multiple introductions by several invasive variants, in competition with a resident population, when observations provide knowledge on the relative proportions of each variant at some dates and places rather than absolute abundances of the different variants. We develop a reaction-diffusion mechanistic-statistical model applied to a genetic spatio-temporal dataset reporting the relative proportions of five genetic variants of watermelon mosaic virus (WMV, genus *Potyvirus*, family *Potyviridae*) in infections of commercial cucurbit fields. Our framework allows us to (1) estimate the dates and places of successful introduction of each emerging variant along with other ecological parameters, (2) reconstruct the invasion history of the emerging variants from their introduction sites, (3) detect competitive advantages of the emerging variants as compared to the resident population, and (4) predict the fate of the different genetic groups, in particular the takeover of the emerging variants over the resident population.

## Material and methods

**Data.** *Pathosystem.* WMV is widespread in cucurbit crops, mostly in temperate and Mediterranean climatic regions of the world¹⁶. WMV has a wide host range including some legumes, orchids and many weeds that can be alternative hosts¹⁶. Like other potyviruses, it is non-persistently transmitted by at least 30 aphid species¹⁶. In temperate regions, WMV causes summer epidemics on cucurbit crops, and it can overwinter in several common non-cucurbit weeds when no crops are present¹⁶,³⁴. WMV has been common in France for more than 40 years, causing mosaics on leaves and fruits in melon, but mostly mild symptoms on zucchini squash. Since 2000, new symptoms were observed in southeastern France on zucchini squash: leaf deformations and mosaics, as well as fruit discoloration and deformations that made them unmarketable. This new agronomic problem was correlated to the introduction of new molecular groups of WMV strains. At least four new groups have emerged since 2000 and they have rapidly replaced the native “classical” strains, causing important problems for the producers³⁵. These new groups, hereafter “emerging strains” (ES) are significantly more related molecularly to worldwide strains than to any other isolates from the French populations³⁶. As emphasised in³⁵, this supports that the new group of emerging strains has arisen through introductions, mostly from Southeastern Asia, rather than through local differentiation.

In this study, we focus on the pathosystem corresponding to a classical strain (CS) and four emerging strains (ES<sub>k</sub>, k = 1, ..., 4) of WMV and their cucurbit hosts.

**Study area and sampling.** The study area, located in Southeastern France, is included in a rectangle of about 25,000 km² and is bounded on the South by the Mediterranean Sea. Between 2004 and 2008, the presence of WMV had been monitored in collaboration with farmers, farm advisers and seed companies. Each year, cultivated host plants were collected in different fields and at different dates between May 1st and September 30th. In total, more than two thousand plant samples were collected over the entire study area. All plant samples were analyzed in the INRAE Plant Pathology Unit to confirm the presence of WMV and determine the molecular type of the virus strain causing the infection (see³⁵ for detail on field and laboratory protocols). All infected host plants were cucurbits, mostly melon and different squashes (e.g., zucchini, pumpkins).

**Observations.** In the absence of individual geographic coordinates, all infected host plants were attributed to the centroid of the municipality (French administrative unit, median size about 10 km²) where they have been collected. Then for one date, one observation corresponded to a municipality in which at least one infected host plant was sampled. Table 1 summarizes for each year, the number of observations (i.e. number of municipalities), the number of infected plants sampled and the proportion of each WMV strain (CS, and ES₁ to ES₄) found in the infected host plants. Errors in assignment of virus samples to the CS or ES strains was negligible because of the large genetic distance separating them: 5 to 10% nucleotide divergence both in the fragment used in the study and in complete genomes³⁵, also precluding the possibility of local jumps between groups by accumulation of mutations.
<table><caption>**Table 1.** Number of observations and corresponding proportions of classical and emerging strains.</caption><thead><tr><th></th><th>2004</th><th>2005</th><th>2006</th><th>2007</th><th>2008</th></tr></thead><tbody><tr><td># observations</td><td>67</td><td>64</td><td>68</td><td>50</td><td>40</td></tr><tr><td># infected samples</td><td>408</td><td>371</td><td>422</td><td>280</td><td>212</td></tr><tr><td>Classical strain (%)</td><td>55</td><td>45</td><td>28</td><td>17</td><td>14</td></tr><tr><td>Emerging strain 1 (%)</td><td>21</td><td>23</td><td>22</td><td>37</td><td>27</td></tr><tr><td>Emerging strain 2 (%)</td><td>13</td><td>18</td><td>23</td><td>21</td><td>32</td></tr><tr><td>Emerging strain 3 (%)</td><td>1</td><td>4</td><td>3</td><td>5</td><td>3</td></tr><tr><td>Emerging strain 4 (%)</td><td>10</td><td>10</td><td>24</td><td>20</td><td>24</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Approximated density $z(x)$ of the host plants in the study area. The density was normalized, so that $z(x) = z(x_1, x_2) = 0$ corresponds to an absence of cucurbit plants and $z(x) = 1$ to the maximum density. The axes $x_1$ and $x_2$ correspond to Lambert93 coordinates (in km). The white regions are non-terrestrial habitats (water). Land use data were not available in the gray regions; the host plant density was then computed by interpolation.</figcaption></figure>

**Landscape.** To approximate the density of WMV host plants over the study area, we used 2006 land use data (i.e. BD Ocsol 2006 PACA and LR) produced by the CRIGE PACA (http://www.crige-paca.org/) and the Association SIG-LR (http://www.siglr.org/association/la-structure.html). Based on satellite images, land use is determined at a spatial resolution of 1/50,000 using an improved three-level hierarchical typology derived from the European Corine Land Cover nomenclature. Here we used the third hierarchical level of the BD Ocsol typology (i.e. 42 land use classes) to classify the entire study area in three habitats: (1) WMV-susceptible crops, (2) habitats unfavorable to WMV host plants (e.g. forests, industrial and commercial units...) and, (3) non-terrestrial habitat (i.e. water). The proportion of WMV-susceptible crops was then computed within all cells of a raster covering the entire study area, with a spatial resolution of 1.4 × 1.4 km². These proportions were used to approximate host plant density $z(x)$, which was normalized, so that $z(x) = 0$ corresponds to an absence of host plants and $z(x) = 1$ to the maximum density of host plants (Fig. 1).

**Mechanistic-statistical model.** The general modeling strategy is based on a mechanistic-statistical approach<sup>12,22,33</sup>. This type of approach combines a mechanistic model describing the dynamics under investigation with a probabilistic model conditional on the dynamics, describing how the measurements have been collected. This method that has already proved its theoretical effectiveness in determining dispersal parameters using simulated genetic data<sup>12</sup> aims at bridging the gap between the data and the model for the determination of virus dynamics.

Here, the mechanistic part of the model describes the spatio-temporal dynamics of the virus strains, given the model parameters (demographic parameters, introduction dates/sites). This allows us to compute the expected proportions of the five types of virus strains (CS and ES<sub>1</sub> to ES<sub>4</sub>) at each date and site of observation. The probabilistic part of the mechanistic-statistical model describes the conditional distribution of the observed proportions of the virus strains, given the expected proportions. Using this approach, it is then possible to derive a numerically tractable formula for the likelihood function associated with the model parameters.

**Population dynamics.** The model is segmented into two stages: (1) the intra-annual stage describes the dispersal and growth of the five virus strains during the summer epidemics on cucurbit crops, and the competition between them, during a period ranging from May 1st (noted $t = 0$) to September 30 (noted $t = t_f, t_f = 153$ days); (2) the inter-annual stage describes the winter decay of the different strains when no crops are present and the virus overwinters in weeds. We denote by $c^n(t, x)$ and $e^n(t, x)$ the densities of classical strain (CS) and emerging strains (ES<sub>k</sub>, $k = 1, \dots, 4$), at position $x$ and at time $t$ of year $n$.
**Dynamics of the classical strain before the first introduction events.** Before the introduction of the first emerging strain, the intra-annual dynamics of the population CS is described by a standard diffusion model with logistic growth: $\partial_t c^n = D\Delta c^n + rc^n(z(\mathbf{x}) - c^n)$. Here, $\Delta$ is the Laplace 2D diffusion operator (sum of the second derivatives with respect to coordinate). This operator describes uncorrelated random walk movements of the viruses, with the coefficient $D$ measuring the mobility of the viruses (e.g.,$^{26,37}$). The term $rz(\mathbf{x})$ is the intrinsic growth rate (i.e., growth rate in the absence of competition) of the population, which depends on the density of host plants $z(\mathbf{x})$ and on a coefficient $r$ (intrinsic growth rate at maximum host density). Under these assumptions, the carrying capacity at a position $\mathbf{x}$ is equal to $z(\mathbf{x})$, which means that the population densities are expressed in units of the maximum host population density. The model is initialized by setting $c^{1980}(0, \mathbf{x}) = (1 - m_c)z(\mathbf{x})$, where $m_c$ is the winter decay rate of the CS (see the description of the inter-annual stage below). In other terms, we assume that the CS density is at the carrying capacity in 1979, i.e., 5 years after its first detection and 20 years before the first detections of ESs$^{38}$.

**Introduction events.** The ESs are introduced during years noted $n_k \ge 1981$, at the beginning of the intra-annual stage (other dates of introduction within the intra-annual stage would lead—at most—to a one-year lag in the dynamics). Their densities are 0 before introduction: $e_k^n = 0$ for $n < n_k$. Once introduced, the initial density of any ES is assumed to be 1/10th of the carrying capacity at the introduction point (other values have been tested without much effect, see Supplementary Fig. S1), with a decreasing density as the distance from this point increases:

$$e_k^n(0, \mathbf{x}) = \frac{z(\mathbf{x})}{10} \exp\left(-\frac{\|\mathbf{x} - \mathbf{X}_k\|^2}{2\sigma^2}\right),$$

where $\mathbf{X}_k$ is the location of introduction of the strain $k$. In our computations, we took $\sigma = 5$ km for the standard deviation.

**Intra-annual dynamics after the first introduction event.** Intra-annual dynamics were described by a neutral competition model with diffusion (properties of this model have been analyzed in [54]):

$$\left\{ \begin{aligned} \partial_t c^n(t, \mathbf{x}) &= D\Delta c^n + rc^n \left( z(\mathbf{x}) - c^n - \sum_{i=1}^4 e_i^n(t, \mathbf{x}) \right), \\ \partial_t e_k^n(t, \mathbf{x}) &= D\Delta e_k^n + re_k^n \left( z(\mathbf{x}) - c^n - \sum_{i=1}^4 e_i^n(t, \mathbf{x}) \right), \end{aligned} \right.$$

for $t = 0... t_f$ and for all introduced emerging strains, i.e. all $k$ such that $n \ge n_k$. We assume reflecting boundary conditions, meaning that the population flows vanish at the boundary of the study area, due to truly reflecting boundaries (e.g., sea coast in the southern part of the site) or symmetric inward and outward fluxes$^{26}$. In addition, in order to limit the number of unknown parameters, and in the absence of precise knowledge on the differences between the strains, we assume here that the diffusion, competition and growth coefficients are common to all the strains during the intra-annual stage (see the discussion for more details on this assumption).

**Inter-annual dynamics.** The population densities at time $t = 0$ of year $n$ are connected with those of year $n-1$, at time $t = t_f$, through the following formulas:

$$\left\{ \begin{aligned} c^n(0, \mathbf{x}) &= (1 - m_c)c^{n-1}(t_f, \mathbf{x}) & \text{for } n \ge 1981 \\ e_k^n(0, \mathbf{x}) &= (1 - m_e)e_k^{n-1}(t_f, \mathbf{x}) & \text{for } n \ge n_k + 1 \end{aligned} \right.$$

with $m_c$ and $m_e$ the winter decay rates of the CS and ESs strains (note that $m_e$ is common to all of the ESs). Estimation of CS and ES decay rates provides an assessment of the competitive advantage of one type of strain vs the other.

**Numerical computations.** The intra-annual dynamics were solved using Comsol Multiphysics time-dependent solver, which is based on a finite element method (FEM). The triangular mesh which was used for our computations is available as Supplementary Fig. S2.

**Probabilistic model for the observation process.** During the years $n = 2004, ..., 2008$, $I_n$ observations were made (see Section *Observations* above and Table 1). They consist in counting data, that we denote by $C_i$ and $E_{k,i}$ for $k = 1, ..., 4$ and $i = 1, ..., I_n$, corresponding to the number of samples infected by the CS and ESs strains, respectively, at each date of observation and location $(t_i, \mathbf{x}_i)$. Note that these dates and locations depend on the year of observation $n$. More generally, the above quantities should be noted $C_i^n, E_{k,i}^n, t_i^n, \mathbf{x}_i^n$; for simplicity, the index $n$ is omitted in the sequel, unless necessary.

We denote by $V_i = C_i + \sum_{k=1}^4 E_{k,i}$ the total number of infected samples observed at $(t_i, \mathbf{x}_i)$. The conditional distribution of the vector $(C_i, E_{1,i}, E_{2,i}, E_{3,i}, E_{4,i})$, given $V_i$ can be described by a multinomial distribution $\mathcal{M}(V_i, \mathbf{p}_i)$ with $\mathbf{p}_i = (p_i^1, p_i^2, p_i^3, p_i^4)$ the vector of the respective proportions of each strain in the population at $(t_i, \mathbf{x}_i)$. We chose to work conditionally to $V_i$ because the sample sizes are not related to the density of WMV.
**Statistical inference.** Unknown parameters. We denote by $\Theta$ the vector of unknown parameters: the diffusion coefficient $D$, the intrinsic growth rate at maximum host density $r$, the winter decay rates ($m_c, m_e$) and the locations ($x_k \in \mathbb{R}^2$) and years ($n_k$) of introduction, for $k = 1, \dots, 4$. Thus $\Theta \in \mathbb{R}^{16}$.

**Computation of a likelihood.** Given the set of parameters $\Theta$, the densities $c^n(t, \mathbf{x}|\Theta)$ and $e_k^n(t, \mathbf{x}|\Theta)$ can be computed explicitly with the mechanistic model for population dynamics. Thus, at a given year $n$, at ($t_i, x_i$), the parameter $\mathbf{p}_i$ of the multinomial distribution $\mathcal{M}(V_i, \mathbf{p}_i)$ writes:

$$p_i^c(\Theta) = \frac{c^n(t_i, \mathbf{x}_i|\Theta)}{c^n(t_i, \mathbf{x}_i|\Theta) + \sum_{i=1}^4 e_i^n(t_i, \mathbf{x}_i|\Theta)}, \quad p_i^e(\Theta) = \frac{e_k^n(t_i, \mathbf{x}_i|\Theta)}{c^n(t_i, \mathbf{x}_i|\Theta) + \sum_{i=1}^4 e_i^n(t_i, \mathbf{x}_i|\Theta)}$$

The probability $P(C_i, E_{1,i}, E_{2,i}, E_{3,i}, E_{4,i}|\Theta, V_i)$ of the observed outcome $C_i, E_{1,i}, E_{2,i}, E_{3,i}, E_{4,i}$ is then

$$P(C_i, E_{1,i}, E_{2,i}, E_{3,i}, E_{4,i}|\Theta, V_i) = \frac{(V_i)!}{C_i! \prod_{k=1}^4 E_{k,i}!} (p_i^c(\Theta))^C_i \prod_{k=1}^4 (p_i^e(\Theta))^{E_{k,i}}$$

Assuming that the observations during each year and at each date/location are independent from each other conditionally on the virus strain proportions, we get the following formula for the likelihood:

$$\mathcal{L}(\Theta) = \prod_{n=2004}^{2008} \prod_{i=1, \dots, I_n} P(C_i, E_{1,i}, E_{2,i}, E_{3,i}, E_{4,i}|\Theta, V_i).$$

**A priori constraints on the parameters.** By definition and for biological reasons, the parameter vector $\Theta$ satisfies some constraints. First, $D \in (10^{-4}, 10)$ km$^2$/day, $r \in (0.1, 1)$day$^{-1}$, and $m_c, m_e \in \{0, 0.1, 0.2, \dots, 0.9\}$, (see Supplementary Note S7 for a biological interpretation of these values). Second, we assumed that the locations of introductions $X_k$ belong to the study area. To facilitate the estimation procedure, the points $X_k$ were searched in a regular grid with 20 points (see Supplementary Fig. S3), and the dates of introduction $n_k$ were searched in {1985, 1990, 1995, 2000}.

**Inference procedure.** Due to the important computation time (4 min in average for one simulation of the model on an Intel(R) Core(R) CPU i7-4790 @ 3.60 GHz), we were not able to compute an a posteriori distribution of the parameters in a Bayesian framework. Rather, we used a simulated annealing algorithm to compute the maximum likelihood estimate (MLE), i.e., the parameter $\Theta^*$ which leads to the highest log-likelihood. This is an iterative algorithm, which constructs a sequence $(\Theta_j)_{j \ge 1}$ converging in probability towards a MLE. It is based on an acceptance-rejection procedure, where the acceptance rate depends on the current iteration $j$ through a "cooling rate" ($\alpha$). Empirically, a good trade-off between quality of optimization and time required for computation (number of iterations) is obtained with exponential cooling rates of the type $T_0\alpha^j$ with $0 < \alpha < 1$ and some constant $T_0 \gg 1$ (this cooling schedule was first proposed in$^{39}$). Too rapid cooling ($\alpha \ll 1$) results in a system frozen into a state far from the optimal one, whereas too slow cooling ($\alpha \approx 1$) leads to important computation times due to very slow convergence. Here, we ran 6 parallel sequences with cooling rates $\alpha \in \{0.995, 0.999, 0.9995\}$. For this type of algorithm, there are no general rules for the choice of the stopping criterion [HenJac03], which should be heuristically adapted to the considered optimization problem. Here, our stopping criterion was that $\Theta_j$ remained unchanged during 500 iterations. The computations took about 100 days (CPU time).

**Confidence intervals and goodness-of-fit.** To assess the model's goodness-of-fit, 95% confidence regions were computed for the observations $(C_i, E_{1,i}, E_{2,i}, E_{3,i}, E_{4,i})$ at each date/location ($t_i, \mathbf{x}_i$), and for each year of observation. The confidence regions were computed by assessing the probability of each possible outcome of the observation process, at each date/location, based on the computed proportions $p_i = (p_i^c, p_i^{e_1}, p_i^{e_2}, p_i^{e_3}, p_i^{e_4})$, corresponding to the output of the mechanistic model using the MLE $\Theta^*$ and given the total number of infected samples $V_i$. Then, we checked if the observations belonged to the 95% most probable outcomes.

## Results

**Convergence and goodness-of-fit.** As expected, the highest likelihood was obtained with a slow cooling rate ($\alpha = 0.9995$). The corresponding MLE, denoted by $\Theta^* = (D^*, r^*, m_c^*, m_e^*, X_1^*, X_2^*, X_3^*, X_4^*, n_1^*, n_2^*, n_3^*, n_4^*)$ is presented in Table 2. In average, 96% of the observations fell within the 95% confidence regions, indicating that the model fits the data well (277 observations over a total of 289; 94% in 2004, 98% in 2005, 96% in 2006, 96% in 2007 and 95% in 2008). We also note that the likelihood function is peaked at $\Theta^*$ in the sense that any perturbation in one component of $\Theta^*$ leads to lower likelihood (see Supplementary Note S8 for likelihood-ratio based confidence intervals and Supplementary Fig. S4 for more details on the profile of the likelihood function). This suggests that the MLE $\Theta^*$ is close to the actual maximizer of the likelihood function.

**Parameter values.** As shown in Table 2, the MLE corresponds to the same date of introduction for $ES_1, ES_2, ES_3$ whereas $ES_4$ has been introduced five years later. Note that a same date of introduction does not mean a same date of detection: depending on local conditions, some strains may establish and spread faster than others (as observed below for $ES_4$).
<table><caption>**Table 2.** Maximum likelihood estimates.</caption><thead><tr><th>Biological parameter</th><th>D*</th><th>r*</th><th>m*<sup>1</sup></th><th>m*<sup>2</sup></th></tr></thead><tbody><tr><td>Value</td><td>0.44 km<sup>2</sup> day<sup>-1</sup></td><td>0.31 day<sup>-1</sup></td><td>0.5 year<sup>-1</sup></td><td>0 year<sup>-1</sup></td></tr><tr><td>Date of introduction</td><td>n<sup>*</sup><sub>1</sub> (ES1)</td><td>n<sup>*</sup><sub>2</sub> (ES2)</td><td>n<sup>*</sup><sub>3</sub> (ES3)</td><td>n<sup>*</sup><sub>4</sub> (ES4)</td></tr><tr><td>Value</td><td>1990</td><td>1990</td><td>1990</td><td>1995</td></tr><tr><td>Site of introduction</td><td>X<sup>*</sup><sub>1</sub></td><td>X<sup>*</sup><sub>2</sub></td><td>X<sup>*</sup><sub>3</sub></td><td>X<sup>*</sup><sub>4</sub></td></tr><tr><td>Value (Lambert 93, km)</td><td>(926, 6369)</td><td>(926, 6369)</td><td>(758, 6369)</td><td>(758, 6294)</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Proportions of the classical and emerging strains in the landscape: data and simulations. The colors of the shaded regions indicate which strain is *the most prevalent*. The red regions correspond to the CS strain; light blue and blue: ES₁, ES₂ (these two strains have the same density, only ES₁ is represented); green: ES₃; pink: ES₄. The pie charts describe the relative proportions of the strains found in the data (same color legend). The white crosses on the 2004 panel represent the estimated sites of introduction. The simulation results presented here correspond to the middle of the intra-annual stage (2nd week of June), and were obtained with the MLE $\Theta^*$.</figcaption></figure>

Regarding the sites of introduction, the MLE indicates that ES₁ and ES₂ have been introduced at the north-eastern corner of the study area, ES₃ in the northwest and ES₄ in the southwest (see the white crosses in the 2004 panel of Fig. 2; see also Supplementary Fig. S3 for more details on the corresponding likelihood). Three introduction points (ESs 1, 2 and 3) were estimated at the edge of the study site, indicating that the introductions may have occurred outside of the study area.

Regarding the biological parameters, the winter decay rate of the CS strain was much higher than that of the ESs strains (0.5 vs. 0), reflecting a high competitive advantage of the ESs. The value $D^* = 0.44 \text{ km}^2 \text{ per day of the diffusion coefficient indicates that a virus travels about 1.2 km per day in average during the growing season of cucurbit crops (see Supplementary Note S7). The growth rate of } 0.31 \text{ day}^{-1} \text{ corresponds to an increase by a factor } e^{0.31} \approx 1.3 \text{ each day, in the absence of competition.}$

*Strain distributions.* Figure 2 depicts the most prevalent strains at each position in the landscape, during four of the five years of observation (2008 is presented in Supplementary Fig. S5), obtained by solving our model with the MLE $\Theta^*$, together with the data. We graphically note a good agreement between the positions of the observed strains and the distributions obtained with the model: ES₁ and ES₂ tend to be distributed in the eastern part of the study area, ES₃ in the northwestern part and ES₄ in the southwestern part, while the CS strain tends to
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Simulated proportions of the classical and emerging strains in the landscape: before and after the observation window. The simulation results presented here correspond to the middle of the intra-annual stage (2nd week of June), and were obtained with the MLE $\Theta^*$. The colors of the shaded regions indicate which strain is the most prevalent. The red regions correspond to the classical strain; blue: $ES_1$, $ES_2$ (these two strains have the same density, only $ES_1$ is represented); green: $ES_3$; pink: $ES_4$.</figcaption></figure>

be progressively confined to the central part of the study area. In 2007, $ES_3$ seems to be more prevalent according to the model than suggested by the observations, probably due to its introduction site, which is far from the observation sites (see the first panel in Fig. 2). Note that, with the deterministic framework used here, as $ES_1$ and $ES_2$ share the same date and position of introduction, and the same parameter values, their distributions are completely equal; thus only $ES_1$ is represented in the figures.

The spatial distributions of the different strains at each year where one of the emerging strains becomes locally more prevalent are depicted in Fig. 3. Although ESs 1, 2 and 3 have been introduced at the same date, their dynamics are influenced by local conditions: ESs 1 and 2 become the most prevalent at least in one part of the study area in 1996 (6 years after their introduction), $ES_3$ in 2001 (11 years after introduction) and $ES_4$ in 2002 (7 years after its introduction). Thus, despite the neutrality assumption, the heterogeneity of the landscape leads to different durations of the establishment stage. The full timeline of the dynamics of the different strain proportions in the landscape, from the first estimated introduction date of an emerging strain (1990) to 2019 is available as Supplementary Fig. S6. Since 2008, due to their competitive advantage (modelled here as a lower winter decay rate), the ESs replaced the CS, which is not anymore the most prevalent strain, whatever the position in the study area, 18 years after the first introduction. Before saturation, the spread rates of the ESs are about 5 km/year (estimated as the slope of $\sqrt{\text{invaded area}}/\pi$, the invaded area corresponding to virus densities > carrying capacity/100). Then, the distribution of the ESs remains almost at equilibrium until the last year of simulation, which is a consequence of the neutrality assumption (equal fitness of all the ESs) (Fig. 3; last panel).

*Average proportions in the study area and effect of the CS on the ESs.* To get a quantitative insight into the replacement of the CS by the ESs, we computed the relative global proportion of each strain by integrating the simulation results (with the parameters corresponding to the MLE $\Theta^*$) over the study area (Fig. 4, panel (A)). Before the first introduction in 1990, the classical strain represents 100% of the infections. In 2010, it represents only 10% of the infections. This decline, which was already visible in the 2004–2008 data$^{34,35}$, is well-captured by the model, though with a slight advance. These discrepancies between the predicted proportions and the data are probably due to the positions of the observation sites, which are concentrated at the center of the domain, where the CS is more prevalent (see Fig. 2). In order to understand the effect of the presence of a resident CS strain on the emerging ESs strains, we compared the dynamics presented in Fig. 4, panel (A) with a hypothetical scenario describing the dynamics of the ESs in the absence of CS. For this, we used the MLE $\Theta^*$, to simulate
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 4.** Estimated average proportions of the classical and emerging strains in the study area. (A) The plain lines correspond to the simulated proportions and the red crosses correspond to the proportions of CS in the data. (B) Simulated proportions of the ESs obtained by assuming that the CS is absent. In both cases, the parameter values correspond to the MLE $\Theta^*$. Note that the curves corresponding to the ESs 1 and 2 are superimposed.</figcaption></figure>

the hypothetical dynamics of the ESs assuming that the CS density is 0. The results are depicted in Fig. 4, panel (B). We observe a very fast convergence to an equilibrium, compared to the situation where the CS is present. Additionally, the last introduced ES ($ES_4$) cannot establish, and $ES_3$ which was confined in an unfavorable region in the presence of the CS, reaches more favorable regions, leading to a higher proportion. Thus, the competition with the CS alters the outcome of the competition between the ESs, and seems to promote the diversity of the ESs by slowing down the overall dynamics.

## Discussion

In this work, we developed a reaction-diffusion model to describe the spatial dynamics of invasion of a resident population inhabiting a spatially structured environment by newly introduced variants. Using a mechanistic-statistical framework, we confronted the model to a dataset recording the invasion by several emerging genetic variants of a resident population of WMV and succeeded in (1) estimating the dates and places of successful introduction of each emerging variant as well as parameters related to growth and dispersal, (2) reconstructing the invasion by the new variants from their introduction sites, (3) establishing a competitive advantage of the new variants as compared to the resident population and (4) predicting the fate of each variant. Simulations with the optimal parameter values showed an adequate fit, proving that the model is able to reproduce the observed spatial dynamics despite the strong mechanistic constraints of the model structure and the strong spatial censorship of the dataset, i.e. missing data. We used fraction frequencies in the dataset rather than counting data. We argue that such observations are more robust to heterogeneities in the sampling conditions as they do not require a standardized observation protocol. Similarly, frequency data should be more robust to heterogeneities in the susceptibility among host genotypes, provided that all of the strains are influenced in the same manner. As stated by⁴⁰, abundance data can only be used if the count-proportion, i.e. the ratio between expected count and population size, can safely be assumed to be constant, or if factors affecting variation in the count-proportion can be identified and then accommodated through parametric modeling.

The estimations suggest that three of the four emerging strains have been introduced at approximately the same date, while the fourth one was introduced 5 years later (the model considered only 5-yr intervals of introductions because of constraints on computation time). Despite the neutrality assumption that we made between the emerging strains, we observed different durations of the establishment stage: while ESs 1 and 2 became locally the most prevalent strains only 6 years after their introduction, $ES_3$ displayed delayed dynamics since it became locally the most prevalent strain 11 years after its introduction. In comparison, $ES_4$ was introduced 5 years later than $ES_3$ but became locally the most prevalent strain at about the same date. The low prevalence of the $ES_3$ in the dataset could be explained by a lower fitness of this strain, e.g. higher winter decay rate, lower growth rate or weaker competitiveness. Our results indicate that this pattern can also be observed with a neutrality assumption, as a result of the joint effects of the local composition of the landscape, and of the position of the sampling sites, far away from the introduction site. Indeed, in the area where $ES_4$ was first observed, cucurbits crops are very frequent with high connectivity between crops, whereas the area where $ES_3$ was found is patchier. Similarly, our results indicate that the overall prevalence of the CS strain in the study area has been slightly overestimated in the data, due to sampling sites concentrated in the regions where it is indeed the most prevalent strain. The reconstructed dynamics of the five strains therefore underline the importance of estimating jointly the places and dates of the introduction and the other ecological parameters as well as the importance of considering the spatial structure of the sampling design.
Based on WMV-infected samples collected in Southeastern France for more than 30 years, ES₁ was first detected in 1999³⁸ whereas the other ESs were observed in 2002–2004³⁵, i.e. 9–14 years after the estimated introduction dates. Such a lapse between the introduction of a plant pathogen and its first detection is consistent with estimations obtained for other plant viruses²¹,⁴¹⁻⁴³. ESs strains have been detected in several European and Mediterranean countries (¹⁶), and in the USA⁴⁴, in the few years following their description in France, and their prevalence in these countries seems to increase even if few time series data are available. The reasons for these almost simultaneous emergences in distant countries and variable environments are not fully understood. WMV being considered so far as not seed-transmitted, the ESs strains have probably been disseminated through long-distance exchanges of plant material³⁵.

In addition to the dates and places of successful introduction, our model provides estimation of ecological parameters in natura. In particular, the diffusion parameter, measuring the mobility of the viruses (or, more precisely, of their aphid vectors) was estimated, leading to a value 0.44km² per day for WMV. Among plant viruses, there are few estimates of diffusion coefficients, and most estimate rather focus on the speed of range expansion, which can be more directly derived from observations. For instance, an average speed of 33 km/yr and 13 km/yr was estimated for the leafhopper-transmitted wheat streak virus⁴⁵ and the whitefly-transmitted East African cassava mosaic virus⁴⁶ respectively, to be compared with the spread rate of 5 km/year found here. For non-persistently aphid-transmitted viruses like WMV and the other potyviruses, the insect remains viruliferous after probing on an infected plant for only a few minutes to hours¹⁶, and the dispersal distance is supposed to be limited, even if in exceptional climatic conditions, the potyvirus maize dwarf mosaic virus was transmitted by viruliferous aphids over more than 1000 km⁴⁷. Estimations for the potyvirus plum pox virus indicated that 50% of the infectious aphids leaving an infected plant land within about 90 m, while about 10% of flights terminate beyond 1 km⁴⁸. Here, a diffusion coefficient of 0.44 km² per day corresponds to a mean dispersal distance of $\sqrt{\pi D} \approx 1.2$ km after 1 day (see Supplementary Note S7), which seems in agreement with these data.

Another critical parameter that was estimated is the winter decay rate. Most studies focus on the epidemic period but off-season dynamics can be crucial to understand demography and genetic diversity³⁴,⁴⁹. We assumed here that emerging strains differed from the resident population only through this parameter. Indeed, no obvious differences in host range, including both cultivated hosts and weeds, have been found between CS and ES strains, whereas ESs strains were found to be better transmitted than CS ones from some weeds infected by both CS and ESs⁵⁰. This could contribute to more efficient transfer from weeds to crops at the beginning of the growing season, leading to a lower winter decay of ESs strains. In our modelling approach we found that ESs strains were able to invade the resident population because of a lower decay during winter. Nevertheless, as all of the other parameters in the model of population dynamics have been set to the same value for the CS and the ESs, all of the competitive advantage of the ESs can be expressed only through the decay rates, explaining these differences, and the unrealistic 0-decay rate of the ESs. The unrealistic null winter decay rate that we found for the ESs suggests that, contrarily to our assumptions, competitive advantage of ESs probably occurs also during summer. Consistently, there is no efficient cross-protection between CS and ESs strains⁵⁰, but²³ found that superinfection by ESs strains of a plant already infected by CS is easier than the opposite situation. In our model, we also make a neutrality assumption between the ESs that differ only by their introduction dates/sites. This assumption may not completely reflect the complexity of the interactions between viral strains and the biological variability between and within molecular groups⁵⁰. Our model could be extended to consider non-neutral competition between the strains and heterogeneous competitive advantages depending on the host genotype. However, such an approach would bring in a new difficulty in that it would involve a larger number of unknown parameters leading to identifiability issues and underdetermination of the model by the data. As is often the case, there is a trade-off between searching for a more realistic description and models that can be trusted because they still match the data satisfactorily with a small set of parameters to identify.

In 2008, the ESs have replaced the CS, which was not anymore the most prevalent strain, whatever the position in the study area, 18 years after the first ESs introduction. Moreover, this competitive advantage of the ESs is expected to lead to the total replacement of the CS by the ESs in about 25 years, i.e. in 2015. These results are consistent with current knowledge: new observations carried out in 2016 and 2017 showed than the classical CS strain is no more detectable⁵¹. Besides the disappearance of CS strains, the surveys performed in 2016–2017 revealed a complex and dynamic situation that fitted partially with the model. ES₃ was detected in only one location in the southwestern part of the study area, confirming its low dispersal and probable low fitness. As predicted by the model, ES₁ and ES₂ were present in the Eastern part of the area and ES₄ was present in all the study area. However, it was found to be more prevalent than ES₁ and ES₂ even in the Eastern part, suggesting that its fitness is higher than ES₁ and ES₂. New variants, not detected in 2004–2008, were also observed in 2016–2017, and some of them presented a high prevalence in all the area. Deep sequencing of two genomic regions revealed, contrary to the 2004–2008 situation⁵², a high prevalence of recombinants among these new strains⁵¹, blurring the distinction between molecular groups based on CP sequences only. As in 2004–2008, landscape heterogeneity seemed to affect virus dispersal⁵¹.

In a more general perspective, this study shows how mechanistic approaches can be used to infer the historical dynamics of invasive genotypes or species from initial introduction. These approaches enable considering hypothetical scenarios, to get a better understanding of the impact of the biological interactions on the overall dynamics. Here, in agreement with theoretical results in⁵³, the simulations without the CS strain showed that its presence promotes a higher diversity among the emerging strains, by altering the outcome of the competition between the ESs, and by slowing down the overall dynamics, thus reducing founder effects. Another advantage of using a mechanistic-statistical approach, compared to a correlative approach, is that the parameter values bring some insight into biological processes and life history traits, (e.g. the diffusion coefficient is related to dispersal ability). Good knowledge of the parameter values, especially for the biological parameters, will be helpful for future modelling, either with reaction-diffusion models or with other approaches such as stochastic diffusion
models, which share some common parameters with reaction-diffusion models (e.g.,¹²). The method developed in this work is computationally very costly. We plan to develop much faster methods, based on deterministic optimization algorithms and analytic descriptions of the gradients of the likelihood. Introducing specific fitness parameters, besides winter decay rate, for the different groups will also help to better understand the effect of interactions between variants on the evolution of viral populations.

Received: 13 July 2020; Accepted: 16 February 2021

Published online: 29 March 2021

## References

1. Seebens, H. et al. No saturation in the accumulation of alien species worldwide. *Nat. Commun.* **8**, 1–9 (2017).

2. Charudattan, R. Biological control of weeds by means of plant pathogens: significance for integrated weed management in modern agro-ecology. *Biocontrol* **46**, 229–260 (2001).

3. Anderson, P. K. et al. Emerging infectious diseases of plants: pathogen pollution, climate change and agrotechnology drivers. *Trends Ecol. Evol.* **19**, 535–544 (2004).

4. Patil, B. L. Plant viral diseases: Economic implications. In *Reference Module in Life Sciences* (Elsevier, 2020).

5. Crowl, T. A., Crist, T. O., Parmenter, R. R., Belovsky, G. & Lugo, A. E. The spread of invasive species and infectious disease as drivers of ecosystem change. *Front. Ecol. Environ.* **6**, 238–246 (2008).

6. Sakai, A. K. et al. The population biology of invasive species. *Ann. Rev. Ecol. Syst.* **32**, 305–332 (2001).

7. Martinetti, D. & Soubeyrand, S. Identifying lookouts for epidemio-surveillance: application to the emergence of Xylella fastidiosa in france. *Phytopathology* **109**, 265–276 (2019).

8. Elena, S. F., Fraile, A. & Garcia-Arenal, F. Evolution and emergence of plant viruses. In Maramorosch, K. & Murphy, F. (eds.) *Advances In Virus Research, Vol 88*, vol. 88 of *Advances in Virus Research*, 161–191 (2014).

9. Mcleish, M. J., Fraile, A. & Garcia-Arenal, F. Trends and gaps in forecasting plant virus disease risk. *Ann. Appl. Biol.* **176**, 102–108 (2020).

10. Rimbaud, L. et al. Using sensitivity analysis to identify key factors for the propagation of a plant epidemic. *R. Soc. Open Sci.* **5**, 171435 (2018).

11. Abboud, C., Bonnefon, O., Parent, E. & Soubeyrand, S. Dating and localizing an invasion from post-introduction data and a coupled reaction-diffusion-absorption model. *J. Math. Biol.* **79**, 765–789 (2019).

12. Roques, L., Walker, E., Franck, P., Soubeyrand, S. & Klein, E. K. Using genetic data to estimate diffusion rates in heterogeneous landscapes. *J. Math. Biol.* **73**, 397–422 (2016).

13. Soubeyrand, S. et al. Inferring pathogen dynamics from temporal count data: the emergence of Xylella fastidiosa in france is probably not recent. *New Phytol.* **219**, 824–836 (2018).

14. Lefeuvre, P. et al. The spread of tomato yellow leaf curl virus from the middle east to the world. *PLoS Pathog.* **6**, e1001164 (2010).

15. Moriones, E., Praveen, S. & Chakraborty, S. Tomato leaf curl New Delhi virus: an emerging virus complex threatening vegetable and fiber crops. *Viruses Basel* **9**, 264 (2017).

16. Lecoq, H. & Desbiez, C. Viruses of cucurbit crops in the mediterranean region: An ever-changing picture. In Loebenstein, G. & Lecoq, H. (eds.) *Viruses and Virus Diseases of Vegetables in the Mediterranean Basin*, vol. 84 of *Advances in Virus Research*, 67–126 (2012).

17. Desbiez, C., Moury, B. & Lecoq, H. The hallmarks of “green” viruses: Do plant viruses evolve differently from the others?. *Infect. Genet. Evol.* **11**, 812–824 (2011).

18. Drummond, A. J., Pybus, O. G., Rambaut, A., Forsberg, R. & Rodrigo, A. G. Measurably evolving populations. *Trends Ecol. Evol.* **18**, 481–488 (2003).

19. Lemey, P., Rambaut, A., Drummond, A. J. & Suchard, M. A. Bayesian phylogeography finds its roots. *PLoS Comput. Biol.* **5**, e1000520 (2009).

20. Pybus, O. G. et al. Unifying the spatial epidemiology and molecular evolution of emerging epidemics. *Proc. Natl. Acad. Sci. USA* **109**, 15066–15071 (2012).

21. Rakotomalala, M. et al. Comparing patterns and scales of plant virus phylogeography: Rice yellow mottle virus in madagascar and in continental Africa. *Virus Evol.* 5, vez023 (2019).

22. Soubeyrand, S. & Roques, L. Parameter estimation for reaction-diffusion models of biological invasions. *Popul. Ecol.* **56**, 427–434 (2014).

23. Fabre, F., Chadouf, J., Costa, C., Lecoq, H. & Desbiez, C. Asymmetrical over-infection as a process of plant virus emergence. *J. Theor. Biol.* **265**, 377–388 (2010).

24. Perefarres, F. et al. Frequency-dependent assistance as a way out of competitive exclusion between two strains of an emerging virus. *Proc. R. Soc. B Biol. Sci.* **281**, 20133374 (2014).

25. Beatriz M. A. & Lopez-Moya, J. When viruses play team sports: Mixed infections in plants. *Phytopathology* **110**, 29–48 (2020)

26. Turchin, P. *Quantitative Analysis of Movement: Measuring and Modeling Population Redistribution in Animals and Plants* (Sinauer, 1998).

27. Hastings, A. et al. The spatial spread of invasions: new developments in theory and evidence. *Ecol. Lett.* **8**, 91–101 (2005).

28. Shigesada, N. & Kawasaki, K. *Biological Invasions: Theory and Practice* (Oxford Series in Ecology and Evolution, Oxford: Oxford University Press, 1997).

29. Berestycki, H., Hamel, F. & Roques, L. Analysis of the periodically fragmented environment model: I—species persistence. *J. Math. Biol.* **51**, 75–113 (2005).

30. Roques, L., Soubeyrand, S. & Rousselet, J. A statistical-reaction-diffusion approach for analyzing expansion processes. *J. Theor. Biol.* **274**, 43–51 (2011).

31. Berliner, L. M. Physical-statistical modeling in geophysics. *J Geophys Res* **108**, 8776 (2003).

32. Soubeyrand, S., Laine, A.-L., Hanski, I. & Penttinen, A. Spatiotemporal structure of host-pathogen interactions in a metapopulation. *Am. Nat.* **174**, 308–320 (2009).

33. Soubeyrand, S., Neuvonen, S. & Penttinen, A. Mechanical-statistical modeling in ecology: from outbreak detections to pest dynamics. *Bull Math Biol* **71**, 318–338 (2009).

34. Lecoq, H., Wipf-Scheibel, C., Nozeran, K., Millot, P. & Desbiez, C. Comparative molecular epidemiology provides new insights into zucchini yellow mosaic virus occurrence in france. *Virus Res.* **186**, 135–143 (2014).

35. Desbiez, C., Joannon, B., Wipf-Scheibel, C., Chandeysson, C. & Lecoq, H. Emergence of new strains of watermelon mosaic virus in South-Eastern France: evidence for limited spread but rapid local population shift. *Virus Res.* **141**, 201–208 (2009).

36. Moury, B., Desbiez, C., Jacquemond, M. & Lecoq, H. Genetic diversity of plant virus populations: towards hypothesis testing in molecular epidemiology. *Adv. Virus Res.* **67**, 49–87 (2006).

37. Murray, J. D. *Mathematical Biology* (Third edition, Interdisciplinary Applied Mathematics 17 (Springer, 2002).
38. Desbiez, C., Costa, C., Wipf-Scheibel, C., Girard, M. & Lecoq, H. Serological and molecular variability of watermelon mosaic virus (genus Potyvirus). *Adv. Virol.* **152**, 775–781 (2007).

39. Kirkpatrick, S., Gelatt, C. D. & Vecchi, M. P. Optimization by simulated annealing. *Science* **220**, 671–680 (1983).

40. Link, W. A. & Sauer, J. R. Estimating relative abundance from count data. *Aust. J. Stat.* **27**, 83–97 (1998).

41. Harkins, G. W. et al. Dating the origins of the maize-adapted strain of maize streak virus, MSV-A. *J. General Virol.* **90**, 3066 (2009).

42. Yasaka, R. et al. The timescale of emergence and spread of turnip mosaic potyvirus. *Sci. Rep.* **7**, 1–14 (2017).

43. Hajizadeh, M., Gibbs, A. J., Amirnia, F. & Glasa, M. The global phylogeny of plum pox virus is emerging. *J. Gen. Virol.* **100**, 1457–1468 (2019).

44. Vincelli, P. & Seebold, K. Report of a watermelon mosaic potyvirus strain in Kentucky undetected by ELISA. *Plant Health Progress* **10**, 47 (2009).

45. Monjane, A. L. et al. Reconstructing the history of maize streak virus strain a dispersal to reveal diversification hot spots and its origin in southern Africa. *J. Virol.* **85**, 9623–9636 (2011).

46. De Bruyn, A. et al. East African cassava mosaic-like viruses from Africa to Indian ocean islands: molecular diversity, evolutionary history and geographical dissemination of a bipartite begomovirus. *BMC Evol. Biol.* **12**, 228 (2012).

47. Zeyen, R., Stromberg, E. & Kuehnast, E. Long-range aphid transport hypothesis for maize dwarf mosaic virus: history and distribution in Minnesota, USA. *Ann. Appl. Biol.* **111**, 325–336 (1987).

48. Pleydell, D. R. et al. Estimation of the dispersal distances of an aphid-borne virus in a patchy landscape. *PLoS Comput. Biol.* **14**, e1006085 (2018).

49. Mcleish, M., Sacristan, S., Fraile, A. & Garcia-Arenal, F. Scale dependencies and generalism in host use shape virus prevalence. *Proc. R. Soc. B Biol. Sci.* **284**, 20172066 (2017).

50. Lecoq, H. et al. Search for factors involved in the rapid shift in watermelon mosaic virus (WMV) populations in south-eastern France. *Virus Res.* **159**, 115–123 (2011).

51. Desbiez, C. et al. Distribution and evolution of the major viruses infecting cucurbitaceous and solanaceous crops in the French Mediterranean area. *Virus Res.* **286**, 198042 (2020).

52. Desbiez, C., Joannon, B., Wipf-Scheibel, C., Chandeysson, C. & Lecoq, H. Recombination in natural populations of watermelon mosaic virus: new agronomic threat or damp squib?. *J. Gen. Virol.* **92**, 1939–1948 (2011).

53. Roques, L., Hosono, Y., Bonnefon, O. & Boivin, T. The effect of competition on the neutral intraspecific diversity of invasive species. *J. Math. Biol.* **71**, 465–489 (2014).

54. Roques, L., Garnier, J., Hamel, F., & Klein, E. K. Allee effect promotes diversity in traveling waves of colonization. *Proc Natl Acad Sci USA* **109** (23), 8828–8833 (2012).

## Acknowledgements

This work was funded by INRAE grant "MEDIA".

## Author contributions

All authors conceived the model; L.R. and J.P. wrote the first draft; L.R. carried out the numerical computations; all authors discussed the results and contributed to the final manuscript.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1038/s41598-021-86314-y.

**Correspondence** and requests for materials should be addressed to L.R.

**Reprints and permissions information** is available at www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_9.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2021