# Phylodynamics and phylogeography of watermelon mosaic virus: Multiple local invasion routes in southern France and recombination-driven limits to global analysis

Patrick Hoscheit, Cecile Desbiez

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

► To cite this version:

Patrick Hoscheit, Cecile Desbiez. Phylodynamics and phylogeography of watermelon mosaic virus: Multiple local invasion routes in southern France and recombination-driven limits to global analysis. Infection, Genetics and Evolution, 2025, 129, pp.105732. 10.1016/j.meegid.2025.105732. hal-04973376

HAL Id: hal-04973376

https://hal.inrae.fr/hal-04973376v1

Submitted on 7 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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>

Distributed under a Creative Commons Attribution 4.0 International License
# Phylodynamics and phylogeography of watermelon mosaic virus: Multiple local invasion routes in southern France and recombination-driven limits to global analysis

Patrick Hoscheit<sup>a,*</sup>, Cécile Desbiez<sup>b</sup>

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

## ARTICLE INFO

### Keywords:

Phylodynamics

Plant viruses

WMV

Molecular epidemiology

## ABSTRACT

Watermelon mosaic virus (WMV) is a major plant pathogen, infecting over 170 plant species, including cucurbits and legumes. Though mostly propagated locally by aphids in a non-persistent manner, long-range dispersal can occur through human-induced plant or vector movements. Understanding patterns of local and global spread of WMV is crucial to help formulate adequate control strategies. We used phylodynamic methods based on partial and whole-genome sequences collected in France between 2000 and 2017 to reconstruct the introduction of new lineages in the past 30 years and their subsequent diffusion in the country. We identified at least 11 different introduction events, hailing from different parts of the global diversity of WMV, highlighting the critical role international exchanges play in the spread of plant pathogens. For three of these lineages, we estimated the time and location of their introduction in the mid-1990s in the south of France and the speed at which they spread in this specific landscape. We also showed that the highly recombinogenic nature of WMV, as with most potyviruses, makes the use of whole genomes necessary to classify these viruses on a global scale and must be taken into consideration to reconstruct viral evolutionary history. Our results demonstrate how genomic sequencing of plant viruses can help reconstruct specific viral outbreaks and understand global circulation patterns of plant pathogens.

## 1. Introduction

Plant viruses represent significant threats to nearly all crop productions. In the absence of curative methods against them, their control relies on prophylaxis in order to prevent virus or vector introductions, control of the viral vectors (insects, mites, nematodes...) in crops, and breeding for plant resistance. In cucurbits, due to the intensification of cultivation practices and long-distance plant material exchanges, the viral pathosystem is constantly changing, with potentially dramatic agronomic consequences (Lecoq and Desbiez, 2012; Desbiez et al., 2020; Desbiez, 2020). Understanding patterns of local and global virus spread is a key step to help formulate adequate control strategies. Among cucurbit-infecting viruses, watermelon mosaic virus (WMV) is an important pathogen worldwide, particularly in temperate and Mediterranean climates. Its rapid population changes in the last decades also make it a very good model to study the evolution of viral populations and the key factors for pathogen emergence. WMV has a broad host

range, including more than 170 species, including carrot, vanilla, ginseng, and numerous weeds that can constitute virus reservoirs (Desbiez and Lecoq, 2021). WMV is a potyvirus with a single-stranded positive-stranded RNA genome of about 10 kb encapsidated in rod-shaped virions circa 700 nm long. The genome is translated as a large polyprotein, cleaved by three viral proteases in ten functional proteins, with an eleventh protein expressed through translational frameshift (Lecoq and Desbiez, 2012). WMV is transmitted non-persistently by more than 30 aphid species, which accounts for its efficient within-field dissemination but not for long-distance exchanges, the interaction between the virus and its aphid vector being very transient.

In France, WMV has been present for more than 50 years. However, molecular epidemiology studies revealed that significant population changes have taken place in the last 25 years: the “classic” (CL) strains present before 2000 were replaced within a few years by “emerging” (EM) strains, probably in relation to a higher fitness of the EM strains at some steps of the infection cycle (Desbiez et al., 2009; Fabre et al., 2010;

<footnote><sup>a</sup> Université Paris-Saclay, INRAE, MaIAGE, 78350 Jouy-en-Josas, France</footnote>

<footnote><sup>b</sup> INRAE, Pathologie Végétale, Montfavet, France</footnote>

<footnote>* Corresponding author at: INRAE - Unité MaIAGE, Domaine de Vilvert, 78352, Jouy-en-Josas Cedex, France.
E-mail address: patrick.hoscheit@inrae.fr (P. Hoscheit).</footnote>
Lecoq et al., 2011). The appearance of EM strains was correlated with more severe symptoms on leaves and fruits of zucchini squash (Desbiez et al., 2009). At least four subgroups of EM strains, identified using their capsid protein (CP) coding region and named EM1 to EM4, were observed between 2000 and 2008 in Southeastern France (Desbiez et al., 2009). There was no obvious difference in fitness or severity between the different EM (Lecoq et al., 2011), but preliminary data suggest that some EM isolates, particularly EM4, could overcome resistance determinants present in the melon germplasm, what may compromise their efficiency and durability if they are used in breeding programs (Lecoq and Desbiez, unpublished). Partial or complete genome sequencing of French isolates collected between 2004 and 2008 showed that recombinants between CL and EM strains, or between different EM strains, were observed in the few years following the introductions but did not emerge durably in crops (Desbiez et al., 2011). Using a reaction-diffusion model coupled with a probabilistic observation model, based on data of EM1 to EM4 subgroup prevalence and geographic distribution from 2004 to 2008 surveys in southern France, provided estimations for the dates and places of introduction of the four subgroups in that area as well as a reconstruction of the dynamics of each of them (Roques et al., 2021). Ten years later, at least seven more groups of strains were observed, and the EM1 to EM4 subgroups displayed very different fates, from successful emergence -as for EM4- to nearly complete extinction for EM3 (Desbiez et al., 2020; Roques et al., 2021). The diffusion-based statistical approach was not applicable to the 2016–2017 dataset since the time course was only two years, and the sampling depth was lower than in 2004–2008. On the other hand, partial and complete sequences of the different variants data were made available. We thus applied phylogeographic approaches, taking into consideration the high frequency of recombination among WMV isolates, to infer the virus global evolution and its dissemination history in France.

## 2. Methods

### 2.1. Whole-genome sequences

We searched for complete genome sequences of WMV on Genbank, which yielded 125 sequences, to which we added three previously unpublished sequences from France, obtained by Sanger sequencing of RT-PCR products as described in (Desbiez and Lecoq, 2008) (see Supplementary Information, Table 1). For 121 of these 128 sequences, a collection date could be found either in the Genbank metadata or the accompanying publications. We manually removed two sequences from the dataset due to 100 % identity with a previously published sequence, as well as three highly divergent sequences (isolates hollyhock (Genbank accession number MK217416), WMV-GZca (MF418043), and WMV-Pg (KX926428)). This yielded a final dataset of 116 complete and dated sequences, varying in length from 9858 nt to 10,145 nt. We aligned the sequences against the reference sequence NC_006262 using ViralMSA (Moshiri, 2021).

### 2.2. Phylogenetic analysis

We used ModelFinder (Kalyaanamoorthy et al., 2017) to select the best-fitting substitution model according to the Bayesian Information Criterion (BIC). Using this model, we constructed a maximum-likelihood tree with IQTREE (Nguyen et al., 2015), which we subsequently used to check for the presence of temporal signal in the set of complete sequences, using TempEst v1.5.3 (Rambaut et al., 2016).

### 2.3. Phylogenetic network analysis and subsampling

A graphical phylogenetic network analysis was performed using Splitstree6 (Huson and Bryant, 2006), with the HKY85 (Hasegawa et al., 1985) distance to compute branch lengths. The Neighbor-Net method (Bryant and Moulton, 2004) was used, and the ensuing network was

<caption>**Table 1**</caption>

Characteristics of 12 recombination-free segments of the WMV genome, as estimated by a phylodynamic molecular clock analysis.

<table><thead><tr><th>Segment</th><th>Positions</th><th>Length (nt)</th><th>Coding region(s)</th><th>Estimated TMRCA (median)</th><th>CL-EM split</th></tr></thead><tbody><tr><td>1</td><td>890-1492</td><td>603</td><td>P1 - HC-Pro</td><td>1657 (95 % CI [1497-1780])</td><td>Yes</td></tr><tr><td>2</td><td>3029-3430</td><td>402</td><td>P3</td><td>1654 (95 % CI [1425-1814])</td><td>No</td></tr><tr><td>3</td><td>4132-4435</td><td>304</td><td>CI</td><td>1582 (95 % CI [1301-1774])</td><td>No</td></tr><tr><td>4</td><td>4664-4914</td><td>251</td><td>CI</td><td>1625 (95 % CI [1361-1797])</td><td>No</td></tr><tr><td>5</td><td>4915-5201</td><td>286</td><td>CI</td><td>1650 (95 % CI [1432-1805])</td><td>No</td></tr><tr><td>6</td><td>5674-5924</td><td>251</td><td>CI - 6 K2</td><td>1640 (95 % CI [1396-1811])</td><td>Yes</td></tr><tr><td>7</td><td>5983-6437</td><td>455</td><td>Nla-VPg</td><td>1589 (95 % CI [1300-1775])</td><td>Yes</td></tr><tr><td>8</td><td>6438-6740</td><td>303</td><td>Nla-VPg - VPg-Pro</td><td>1553 (95 % CI [1227-1772])</td><td>Yes</td></tr><tr><td>9</td><td>6741-7157</td><td>417</td><td>VPg-Pro</td><td>1604 (95 % CI [1340-1779])</td><td>Yes</td></tr><tr><td>10</td><td>7471-7769</td><td>298</td><td>NIb</td><td>1765 (95 % CI [1603-1863])</td><td>Yes</td></tr><tr><td>11</td><td>7936-8202</td><td>267</td><td>NIb</td><td>1773 (95 % CI [1633-1871])</td><td>Yes</td></tr><tr><td>12</td><td>9044-9654</td><td>611</td><td>CP</td><td>1748 (95 % CI [1589-1861])</td><td>No</td></tr></tbody></table>

displayed using the Splits Network method. We then used this splits tree to select optimal subsets of sequences with sizes varying from 10 to 110 in order to maximise phylogenetic split diversity using the PDA software (Chernomor et al., 2015). We selected the optimal subset of sequences such that 95 % of split diversity is conserved, then added sequences to ensure that all 11 groups identified in France in 2016–2017 (Desbiez et al., 2020), as well as the reference sequence WMV-Fr (accession number AY437609), were represented. This yielded a dataset of 62 sequences.

### 2.4. Recombination analysis

We checked the 62 sequences for signals of recombination using RDP5 (Martin et al., 2021), including the RDP, GENECONV, Chimaera, MaxChi and 3SEQ methods for detection of recombination events, as well as BootScan and SiScan for secondary checking of putative events detected by other methods. Only events confirmed by at least four methods, at a significance level of $p = 10^{-5}$ were included in the final set of selected events. We used a breakpoint distribution plot to check for recombination hot- and cold spots along the genome of WMV and extracted sequence alignments for recombination-free regions of the genome. We kept the twelve regions with length > 250 nt for a combined phylogenetic analysis using BEAST1.10.4, fitting all segments with a common molecular clock. We performed the analysis using a strict and a (uncorrelated lognormal) relaxed clock, selecting the best-fitting model using stepping-stone analysis. We finally built a tanglegram representing the differences in the estimated MCC time trees using Espalier (Rasmussen and Guo, 2023).

### 2.5. Global diversity analysis

We performed a recombination-aware phylodynamic analysis of the 62 complete WMV sequences using BEAST2.7.4 (Bouckaert et al., 2014) with the Recombination package (Müller et al., 2022). We used a fixed value of $3.0 \times 10^{-4}$ substitutions/site/year for the molecular clock, consistent with the value estimated from the nonrecombining segments (see Results) and ran four independent MCMC chains, sampling every 10,000 steps. After excluding the first 10 % of the samples of each chain
for MCMC burn-in, the chains were combined using LogCombiner v2.7.4. This procedure yielded a posterior sample of 20,000 recombination networks. This sample was summarised as a maximum clade credibility (MCC) network using the internal tool from the *Recombination* package. We also extracted the local trees corresponding to every position in the genome, enabling us to compute statistics such as the time to the Most Recent Common Ancestor (MRCA) of the 62 sequences.

## 2.6. Partial genome sequences

We selected all samples collected in the south of France between 2004 and 2017 (in two separate campaigns 2004–2008 and 2016–2017) for which both a fragment from the cylindrical inclusion (CI) gene and a fragment from the CP gene were unambiguously available. When identical samples were collected with the same sampling date and location, we included only one sequence for each such triple (sequence, date, location). Sequences from 2004 to 2008 were assigned the “commune” (district) in which they were collected, whereas sequences from 2016 to 2017 included GPS coordinates of the plant from which they were sampled. This yielded 191 sequences from the EM1 group, 79 from the EM2 group, and 147 from the EM4 group. There were not enough observations of the EM3 group in 2016–2017 (only a single observation in 2016) to meaningfully assess its spread over the whole observation period. We also assigned to each sequence a file containing the boundary shape of the commune in which it was sampled to serve as uniform prior for its sampling location (see “Local phylogeographic analysis”). The shapes were obtained from the public GEOFLA database (Nguyen et al., 2015) in its 2015 version.

## 2.7. Local phylogeographic analysis

We performed a phylogeographic analysis of the 417 partial sequences using BEAST1.10.4, linking all three clades through a common molecular clock and a common substitution model. We used a strict clock model with a lognormal prior (mean $4.0 \times 10^{-4}$ substitutions/site/year), consistent with what was estimated for potyviruses (Gibbs et al., 2008; Gibbs et al., 2020) and a HKY + Γ4 substitution model. We chose a relaxed random walk with Cauchy-distributed rates to reflect possible inhomogeneity in diffusion rates along the branches of the tree. We used a uniform prior for each sequence location within the boundaries of the commune from which it was sampled, reflecting uncertainty about exact sampling locations (Dellcour et al., 2020). We combined eight independent runs of the MCMC chain, each with $5 \times 10^8$ states, sampled every 50,000 steps, checking that they converged to the same posterior. We analysed the posterior sample of trees using TreeAnnotator v1.10.4 to obtain 80 % HPD (highest posterior density) estimates of root location and MRCA times for all three clades.

We also assessed the presence of a temporal signal in this dataset by performing a date randomisation test (Ramsden et al., 2008; Duchène et al., 2015) by randomly reassigning the sampling dates to the sequences and running the phylogeographic analysis five times to obtain a clock rate estimate in the absence of phylogenetic temporal signal, then comparing the range of clock estimates to the randomized estimates. We tested for confounding between the temporal and the genetic structure of the data using a Mantel permutation test (Murray et al., 2016). Based on the MCC tree estimated from the data with the true dates, we defined clusters from the sampling dates, and performed date-randomizations between clusters in order to correct for the biases induced by confounding. Following (Duchène et al., 2015), we applied both the CR1 and the more stringent CR2 test. The CR1 test consists in comparing the mean clock rate estimate with the correct dates to the 95 % credible intervals of the same statistic estimated from the date-randomized samples. The CR2 test checks whether the 95 % credible interval of the non-randomized analysis overlaps with the 95 % credible intervals of the date-randomized analysis.

# 3. Results

## 3.1. Global diversity

Analysis of the temporal signal present in the dataset of all complete genomes using TempEst (see Supplementary Information) clearly showed clock-like behaviour (root-to-tip regression coefficient $R^2 = 0.37$), although the estimated clock rate was far too high ($4.6 \times 10^{-3}$ substitutions/site/year) compared to previous estimates obtained for potyviruses and RNA viruses in general (Gao et al., 2020; Kawakubo et al., 2021). This probably reflects the high levels of recombination present in the dataset, which tends to bias maximum-likelihood tree estimates, on which root-to-tip regression methods are based. Indeed, a graphical analysis of the splits tree (Fig. 1) shows extensive recombination in the WMV genome, with no clear cluster structure apparent. This is consistent with the high number of recombination events identified by RDP5. Indeed, 105 potential events were flagged by at least one detection method, of which 82 were detected by at least four methods. None of the sequences were considered non-recombinant by RDP5, making a classical phylogenetic analysis discarding recombinant sequences impossible, as was already noticed in previous work (Ben Mansour et al., 2023).

The recombination breakpoint plot (Fig. 2) shows recombination occurring in all ten genetic regions, with six hotspots located respectively in the P1 gene, in the HC-Pro gene, the 3' terminal end of the P3 gene, and in the NIB gene.

Splitting the genome into recombination-free regions yielded twelve segments longer than 250 nt, with a total length of 4448 nt, a little under half of the whole genome (Table 1).

A joint phylogenetic analysis of these segments using a common molecular clock model showed strong incongruence of local trees, even in immediately adjacent regions (Fig. 3). A path-sampling/steppingstone analysis strongly favoured the relaxed clock model over the strict clock (log-marginal likelihood -17,560 vs -17,583). The mean clock rate was estimated at $2.60 \times 10^{-4}$ substitutions/site/year (95 % CI [1.93 $\times 10^{-4}$, 3.29 $\times 10^{-4}$]). The median time to the most recent common ancestor varied between the segments, with the oldest estimated around 1553 and the youngest around 1773, albeit with wide credible intervals (Table 1).

The incongruence of local trees makes the determination of migration patterns difficult: for instance, the lineage of FMF00-LL1, first identified in France in 2000, clusters with two sequences found in South Korea (TA-om3 and GC-wm1) for nine out of the twelve local trees (Table 2), thus possibly indicating a south-east Asian origin for this part of the genome. However, for segment 6, located across the CI and 6 K2 genes, FMF00-LL1 was most closely related to the RKG2 lineage found in India in 2014. This pattern, while still broadly consistent with an Asian origin for this specific lineage, suggests a complex history of migrations and recombinations, which the currently available dataset cannot fully recover. Only one of the 11 French lineages had a consistent closest related lineage across all twelve segments (EM170501, related to the American lineage TX113).

The Bayesian recombination analysis of 62 complete genomes confirms the high levels of recombination in this dataset, with a median of 293 recombination events identified in the ancestral history of the sample. The estimated median recombination rate (conditionally on the fixed substitution rate of $3.0 \times 10^{-4}$ substitutions/site/year) was $5.69 \times 10^{-6}$ recombinations per site per year (95 % CI [3.72 $\times 10^{-6}$, 7.61 $\times 10^{-6}$]), which translates to around 0.05 recombinations per lineage per year. This rate is comparable to the estimated recombination rate for coronaviruses (Kalyaanamoorthy et al., 2017). Analysis of the MCC recombination network (Fig. 4) shows pervasive recombination with no clear clade structure.

Even though the MCC network as a whole is statistically well-supported, with many branches having high posterior probability, the local trees at given positions of the genome have low phylogenetic
<figure><img src="image_4.png" /><figcaption>**Fig. 1.** Splits tree of the 116 complete genome sequences of WMV. The 11 new lineages observed in France between 2004 and 2017 are marked in red. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Top: Distribution of recombination hotspots (in red) and recombination-free regions longer than 250 nt (in orange). Bottom: Recombination breakpoint plot of the 62 subsampled WMV complete sequences, computed using RDP5 (Martin et al., 2021). The solid black line represents the number of recombination breakpoints within a sliding 200 nt window. The dark grey and light grey shaded regions represent 95 % and 99 % confidence intervals for the local cold/hotspot tests. The check marks on top of the figure represent estimated recombination breakpoints, and the dotted lines are 95 % and 99 % confidence thresholds for the global cold/hotspot tests (a precise description of these quantities can be found in the reference manual of RDP5). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

support, indicated by very low posterior probabilities, especially for the deepest branches in the tree. Thus, while the complete Bayesian analysis under the coalescent with recombinations yields statistically well-

supported posteriors for relevant parameters of interest, such as recombination rates or divergence dates between given sets of sequences, there is a significant amount of topological uncertainty. This
<figure><img src="image_6.png" /><figcaption>Fig. 3. Tanglegram representing relationships between consecutive maximum clade credibility (MCC) trees of 12 recombination-free segments of the WMV genome.</figcaption></figure>

<table><caption>**Table 2**

Characteristics of the 11 lineages identified in France between 2004 and 2017.</caption><thead><tr><th>Group</th><th>Accession number</th><th>Date of first sample</th><th>Closest related lineage (segments)</th></tr></thead><tbody><tr><td>Gr1</td><td>EU660581 (FMF00-LL1)</td><td>2000</td><td>RKG2 (6), TA-om3 (3-5, 7, 9-12), GC-wm1 (8-12)</td></tr><tr><td>Gr2</td><td>EU660583 (FMF03-141)</td><td>2003</td><td>TX29 (1-6, 9, 12), RKG2 (7-8, 11)</td></tr><tr><td>Gr3</td><td>EU660586 (FBR04-37)</td><td>2004</td><td>TX113 (1-2), WMV-CHN (3-4, 8-9), Henan (5)</td></tr><tr><td>Gr4</td><td>EU660585 (C05-270)</td><td>2005</td><td>Yeongam4-2 (2), RKG2 (3), TX29 (4), RKG (8), NJ-me10 (10), TX113 (1)</td></tr><tr><td>Gr5</td><td>MN914158 (C08-451)</td><td>2008</td><td>WMV-Smecko (1, 3, 5-10, 12), NJ-me10 (4, 6), WMV_1_SDE1 FF (2)</td></tr><tr><td>Gr6</td><td>MN914159 (Cg15-5055)</td><td>2015</td><td>Urumqi (2, 8-9, 11), WMV-WS (2-3, 6-10)</td></tr><tr><td>Gr7</td><td>MN914160 (EM160093)</td><td>2016</td><td>WMv-Li (4-5, 6), Henan (6-7), Sangju2_2014 (8), MeWM7 (9-10), CN-sm2-2 (1-2)</td></tr><tr><td>Gr8</td><td>PP505905 (EM170143)</td><td>2017</td><td>{WMV-WS, Urumqi} (1-2, 5, 9-11), Urumqi (3-4, 7)</td></tr><tr><td>Gr9</td><td>PP505906 (EM170465)</td><td>2017</td><td>TX113 (1, 8), Urumqi (2-5, 7, 12), WMV-CHN (9-11)</td></tr><tr><td>Gr10</td><td>PP505907 (EM170501)</td><td>2017</td><td>TX113 (1-12)</td></tr><tr><td>Gr11</td><td>EU660587 (C06-188)</td><td>2006</td><td>CHI87-620 (2, 12), IR02-54 (1, 3, 5, 9), CHI02-481 (8), TURK91 (4), VE10-099 (6), RKG (11)</td></tr></tbody></table>

uncertainty reflects the statistical difficulty of the task due to the enormous size of phylogenetic network space, as well as the challenging nature of this specific dataset, which was already apparent from the highly reticulated nature of the splits tree estimated earlier (Fig. 1). For example, we assessed the support of the deep split between CL and EM

sequences, originally described based on a 267–273 nt fragment of the N-terminal part of the CP gene (not covered by any of the twelve non-recombining regions of the genome analysed in this work). Interestingly, this split also appears when considering full genomes, using either maximum-likelihood phylogenetic reconstruction (Rabadán and Gómez, 2023) or full Bayesian analysis under the assumption of no recombination (results not shown). Using the WMV-Fr isolate as basal representative of the CL strains and the four reference sequences (FMF00-LL1, FMF03-141, FBR04-37 and C05-270) for the EM strains, we computed the posterior support of the {WMV-Fr|FMF00-LL1,FMF03-141, FBR04-37,C05-270} split in the local trees at positions across the genome, every 10 nt. We found consistently low support (<25 %) with some local variability according to the position (Fig. 5). However, seven of the MCC trees estimated from the nonrecombining segments were compatible with this split (Table 1), indicating a more complex evolutionary relationship between the strains identified as “classical” and those identified as “emergent”.

Analysing the local trees at each position in the genome, we also computed the posterior distribution of the time to the MRCA (TMRCA) of the whole population, comparing it to the times obtained by analysing the twelve nonrecombining segments (Fig. 6). Overall, we observed that TMRCAs estimated from the joint analysis of the twelve segments were comparable to those obtained from the whole genome analysis, with the median MRCA estimate from the recombination-aware analysis within the credibility intervals of the joint analysis for 9 out of 12 segments. The global analysis estimates the MRCA of the sample of 62 sequences to have occurred about 370 years before the present (the most recent sequence in the sample having been collected in 2020), with wide credible intervals, especially in the 3'-terminal region of the P1 gene.

## 3.2. Spread of EM1, EM2 and EM4 in the south of France

The combined phylogeographic analysis showed similar MRCA times
<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Maximum clade credibility network (drawn as a cladogram, using Icytree), obtained from a posterior sample of networks, representing the ancestral recombination history of 62 complete WMV genomes.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Posterior support for the CL-EM split according to the position in the genome.</figcaption></figure>

for EM1, EM2 and EM4 clades (Fig. 7), namely 1995 (95 % CI [1988,1999]), 1994 (95 % CI [1987,1999]) and 1994 (95 % CI [1988,1999]). These estimates are consistent with the analysis from (Roques et al., 2021), who found introduction dates of 1990, 1990 and 1995, respectively, using a maximum-likelihood framework based on a spatiotemporal diffusion model with interacting strains. The posterior density for the molecular clock rate had a median of $3.41 \times 10^{-4}$ sub-stitutions/site/year (95 % CI [2.7 $\times 10^{-4}$,4.2 $\times 10^{-4}$]). This estimate is consistent with the relaxed-clock estimates obtained from the non-recombinant segments of whole genomes. While the Mantel test revealed significant confounding between the genetic and temporal structure of the data ($p = 0.001$), the dataset still passed the clustered

date randomisation test in both the CR1 and the more stringent CR2 version (Duchène et al., 2015; Murray et al., 2016), indicating the presence of sufficient temporal signal to make meaningful phylogenetic analyses (Supplementary Information).

Our analysis showed likely introductions of EM1 in the north of the Bouches-du-Rhône department, EM2 in Vaucluse and EM4 in the eastern part of Gard (Fig. 8). Uncertainty regarding the root location of EM2, in particular, was high due to the low number of available sequences. The estimated dispersal rate (Fig. 7) was markedly lower for EM2 (median 5.05 km/year, 95 % CI [3.7,6.6]) than for EM1 (median 6.9 km/year, 95 % CI [5.4, 8.4]) and EM4 (median 6.8 km/year, 95 % CI [5.3,8.6]), possibly reflecting lower intrinsic fitness of the EM2 clade as well as less
<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Time to the most recent common ancestor of the whole sample of 62 sequences, as estimated from the global analysis using the coalescent with re-combinations (red) and from the tree-based analysis of the twelve nonrecombining segments (green). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Plots of the posterior distributions of time to the most recent common ancestor (left) and diffusion rate in km/year (right).</figcaption></figure>

favourable landscape composition. Indeed, EM2 was not detected in
limited surveys performed after 2017 (data not shown).

**4. Discussion**

In this work, we investigated the diversity of WMV sequences found
in France over two sampling periods, between 2004 and 2008 and
2016–2017. We showed that, like many plant viruses and especially
potyviruses, WMV is recombining at high rates, which makes it neces-
sary to use methods that explicitly take this aspect into account when
estimating phylogenetic trees representing the ancestral history of a
sample. In particular, we showed that the division of WMV strains into

"classical" (CL) and "emergent" (EM) groups, which was initially based
on phylogenetic analysis of a fragment of the CP gene, is not well-
supported over the whole genome. This raises the question of the defi-
nition of a robust strain taxonomy for viruses with high recombination
rates, such as WMV.

Methods aiming to reconstruct the entire Ancestral Recombination
Graph (ARG) from genomic data have been developed for human data
(Rasmussen et al., 2014; Kelleher et al., 2019; Mahmoudi et al., 2022)
but come at a significant computational cost. This is also true for the
recombination-aware inference framework introduced in (Müller et al.,
2022), based on Hudson's coalescent with recombination model. In this
work, we used this method to obtain phylodynamical estimates of
<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Median locations and 80 % posterior density areas of the root location of EM1 (orange), EM2 (green) and EM4 (purple) clades in the south of France. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

quantities such as the overall recombination rate for WMV or the time to the most recent common ancestor of the 62 samples in our study. We found that, while such low-dimensional statistics were robustly estimated from our data, the topology of the ARG was affected by high uncertainty. This is consistent with results in the literature (Hayman et al., 2023; Brandt et al., 2024) showing that, at least theoretically, ARG topologies might not be identifiable from genome data as soon as the sample size exceeds a few dozen. Hence, the posterior distribution is spread out over a large set of phylogenetic networks while still yielding low uncertainty estimates for summary statistics.

This work highlights several statistical challenges related to phylogenetics of viruses undergoing frequent recombination. Indeed, current statistical frameworks used in Bayesian phylogenetics, such as *Recombination* (Müller et al., 2022) or *Bacter* (Vaughan et al., 2017), are limited to analysing a few dozen sequences. This precludes scaling up analyses such as the one presented in this work to larger samples, which will undoubtedly be available in coming years and, thus, limits the use of recombination-aware inference methods to datasets with limited or absent recombination. Areas for future improvement include the development of efficient parallelisation methods for likelihood computations, approximation schemes to the full coalescent with recombination, or the design of well-mixing operators in the space of phylogenetic networks, such as the “Kuhner move” proposed in (Mahmoudi et al., 2022). There are also significant challenges in visualising and summarising posterior distributions over the space of phylogenetic networks. Indeed, these are complex high-dimensional objects, and questions of biological interest (such as the existence of the CL-EM split investigated in this work) require software capable of integrating over the outputs of Bayesian methods.

The estimates of divergence times found in this work differ from the analysis made by Rabadán and Gomez (Rabadán and Gómez, 2023;

Rabadán, 2025) using the *nextstrain* sequence analysis framework (Hadfield et al., 2018), even though the latter analysis did not explicitly take recombination into account. Indeed, for the 62 sequences considered in this work, we estimate the time of the most recent common ancestor of around 370 years before the present across most of the genome (Fig. 6), corresponding to a date of around 1650 CE. By contrast Rabadán and Gomez (Rabadán and Gómez, 2023) found a date of 1474 CE (95 % confidence interval [1081,1619]), as displayed on the most recent version of the WMV Nextstrain tree). Using a subtree dating method, Ben Mansour et al. (Ben Mansour et al., 2023) found a divergence date of around 1014 CE. The discrepancy between these numbers and the more recent ones we found using the method from (Müller et al., 2022) can probably be ascribed to the confounding effects of recombination on temporal signals in viral genomic datasets.

Several works have used the spatial locations of viral samples to study the geographical spread of potyviruses, hence leveraging the additional signal contained in these metadata to further inform phylogenetic tree estimates. Most notably, Kawakubo et al. (Kawakubo et al., 2021) managed to recover patterns of Eurasian spread of turnip mosaic virus (TuMV) from a sample of 153 sequences, obtained by removing all identified recombinants. The same method was used to study global migration dynamics of potato virus Y (PVY) (Gao et al., 2020). Incorporating spatial metadata in recombination-aware phylodynamic analyses would be a significant step forward, enabling the replication of such analyses for highly recombinogenic viruses such as WMV, and also allowing for the use of all collected data instead of having to remove obvious recombinants from the datasets, thereby potentially biasing obtained estimates and increasing uncertainty.

It is striking that the 11 groups identified in France between 2004 and 2017 span almost the entire known variability of WMV. This likely reflects high levels of interconnectivity of plant trading networks
worldwide, in line with previous studies on WMV and other potyviruses (Kawakubo et al., 2021; Rabadán and Gómez, 2023), but a complete phylogeographical analysis taking recombination effects into account is out of reach of current methods, especially when considering the high topological uncertainty in the inferred phylogenetic network, precluding the use of tools such as SCAR (Guo et al., 2022).

We also showed how continuous Bayesian phylogeographical methods can be used to reconstruct the invasion dynamics of three strains, leveraging historical sequence data to obtain estimates of time and place of introduction, as well as the rate of spatial diffusion. This approach can be contrasted with results in (Roques et al., 2021), which used only the strain occurrence data to calibrate a spatiotemporal diffusion model. The introduction dates we find are close to the dates estimated using this model, but we were able to use the phylogenetic information present in the sample to obtain estimates of the location of the introductions that appear more realistic, as they are at the centre of the areas where the different EM strains were first observed.

## Data statement

All complete sequences used in this study are available on Genbank, including 3 previously unpublished sequences (accession numbers PP505905-PP505907). All other data are available upon request from the corresponding author PH. The data are not publicly available due to them containing information that could compromise the privacy of participants.

## CRediT authorship contribution statement

**Patrick Hoscheit:** Writing – review & editing, Writing – original draft, Validation, Formal analysis, Data curation, Conceptualization.
**Cécile Desbiez:** Writing – review & editing, Writing – original draft, Validation, Funding acquisition, Data curation, Conceptualization.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgments

We gladly acknowledge financial support from the INRAE metaprogram SuMCrop "Sustainable crop health management". We are also grateful to the Institut Français de Bioinformatique and the « IFB-core » platform, financed under the Programme d'Investissements d'Avenir funded by the Agence Nationale de la Recherche (RENABI-IFB ANR-11-INBS-0013 and MUDIS4LS ANR-21-ESRE-0048), for providing computing and storage resources.

## Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.meegid.2025.105732.

## Data availability

Data will be made available on request.

## References

* Ben Mansour, K., Gibbs, A.J., Komínková, M., Komínek, P., Brožová, J., Kazda, J., et al., 2023. Watermelon mosaic virus in the Czech Republic, its recent and historical origins. Plant Pathol. 72 (8), 1528-1538.

* Bouckaert, R.R., Heled, J., Kühnert, D., Vaughan, T.G., Wu, C.H., Xie, D., et al., 2014. BEAST 2: a software platform for Bayesian evolutionary analysis. PLoS Comput. Biol. 10 (4), e1003537.

* Brandt, D.Y.C., Huber, C.D., Chiang, C.W.K., Ortega-Del, Vecchio D., 2024. The promise of inferring the past using the ancestral recombination graph. Van Dorp L, editor. Genome Biol. Evol. 16 (2), evae005.

* Bryant, D., Moulton, V., 2004. Neighbor-net: an agglomerative method for the construction of phylogenetic networks. Mol. Biol. Evol. 21 (2), 255-265.

* Chernomor, O., Minh, B.Q., Forest, F., Klaere, S., Ingram, T., Henzinger, M., et al., 2015. Split diversity in constrained conservation prioritization using integer linear programming. Methods Ecol. Evol. 6 (1), 83-91.

* Dellicour, S., Lemey, P., Artois, J., Lam, T.T., Fusaro, A., Monne, I., et al., 2020. Incorporating heterogeneous sampling probabilities in continuous phylogeographic inference — application to H5N1 spread in the Mekong region. Schwartz R, editor. Bioinformatics 36 (7), 2098-2104.

* Desbiez, C., 2020. The never-ending story of cucurbits and viruses. Acta Hortic. 1294, 173-192.

* Desbiez, C., Lecoq, H., 2008. Evidence for multiple intraspecific recombinants in natural populations of watermelon mosaic virus (WMV, Potyvirus). Arch. Virol. 153 (9), 1749-1754.

* Desbiez, C., Lecoq, H., 2021. Watermelon mosaic virus and zucchini yellow mosaic virus (Potyviridae). In: Bamford, D.H., Zuckerman, M. (Eds.), Encyclopedia of Virology, (Fourth Edition) [Internet]. Academic Press, Oxford, pp. 862-870 [cited 2024 Sep 2] Available from. https://www.sciencedirect.com/science/article/pii/B9780128096338212478.

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

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

* Desbiez, C., Wipf-Scheibel, C., Millot, P., Berthier, K., Girardot, G., Gognalons, P., et al., 2020. Distribution and evolution of the major viruses infecting cucurbitaceous and solanaceous crops in the French Mediterranean area. Virus Res. 286, 198042.

* Duchene, S., Duchene, D., Holmes, E.C., Ho, S.Y.W., 2015. The performance of the date-randomization test in phylogenetic analyses of time-structured virus data. Mol. Biol. Evol. 32 (7), 1895-1906.

* Fabre, F., Chadœuf, J., Costa, C., Lecoq, H., Desbiez, C., 2010. Asymmetrical over-infection as a process of plant virus emergence. J. Theor. Biol. 265 (3), 377-388.

* Gao, F., Kawakubo, S., Ho, S.Y.W., Ohshima, K., 2020. The evolutionary history and global spatio-temporal dynamics of potato virus Y. Virus Evol. 6 (2), 1-10.

* Gibbs, A.J., Ohshima, K., Phillips, M.J., Gibbs, M.J., 2008. The prehistory of potyviruses: their initial radiation was during the dawn of agriculture. PLoS ONE 3 (6).

* Gibbs, A.J., Hajizadeh, M., Ohshima, K., Jones, R.A.C., 2020. The potyviruses: an evolutionary synthesis is emerging. Viruses 12 (2), 132.

* Guo, F., Carbone, I., Rasmussen, D.A., 2022. Recombination-aware phylogeographic inference using the structured coalescent with ancestral recombination. PLoS Comput. Biol. 18 (8), e1010422.

* Hadfield, J., Megill, C., Bell, S.M., Huddleston, J., Potter, B., Callender, C., et al., 2018. Nextstrain: real-time tracking of pathogen evolution. Bioinformatics 34 (23), 4121-4123.

* Hasegawa, M., Kishino, H., Yano, T., aki., 1985. Dating of the human-ape splitting by a molecular clock of mitochondrial DNA. J. Mol. Evol. 22 (2), 160-174.

* Hayman, E., Ignatieva, A., Hein, J., 2023. Recoverability of ancestral recombination graph topologies. Theor. Popul. Biol. 154, 27-39.

* Huson, D.H., Bryant, D., 2006. Application of phylogenetic networks in evolutionary studies. Mol. Biol. Evol. 23 (2), 254-267.

* Kalyaanamoorthy, S., Minh, B.Q., Wong, T.K.F., von Haeseler, A., Jermiin, L.S., 2017. ModelFinder: fast model selection for accurate phylogenetic estimates. Nat. Methods 14 (6), 587-589.

* Kawakubo, S., Gao, F., Li, S., Tan, Z., Huang, Y.K., Adkar-Purushothama, C.R., et al., 2021. Genomic analysis of the brassica pathogen turnip mosaic potyvirus reveals its spread along the former trade routes of the silk road. Proc. Natl. Acad. Sci. 118 (12), e2021221118.

* Kelleher, J., Wong, Y., Wohns, A.W., Fadil, C., Albers, P.K., McVean, G., 2019. Inferring whole-genome histories in large population datasets. Nat. Genet. 51 (9), 1330-1338.

* Lecoq, H., Desbiez, C., 2012. Chapter 3 - Viruses of cucurbit crops in the Mediterranean region: an ever-changing picture. In: Loebenstein, G., Lecoq, H. (Eds.), Advances in Virus Research [Internet]. Academic Press, pp. 67-126 [cited 2024 Sep 2]. (Viruses and Virus Diseases of Vegetables in the Mediterranean Basin; vol. 84). Available from. https://www.sciencedirect.com/science/article/pii/B9780123943149000038.

* Lecoq, H., Fabre, F., Joannon, B., Wipf-Scheibel, C., Chandeysson, C., Schoeny, A., et al., 2011. Search for factors involved in the rapid shift in watermelon mosaic virus (WMV) populations in South-Eastern France. Virus Res. 159 (2), 115-123.

* Mahmoudi, A., Koskela, J., Kelleher, J., Chan, Y., Ban Balding, D., 2022. Bayesian inference of ancestral recombination graphs. Kosakovsky Pond SL, editor. PLOS Comput. Biol. 18 (3), e1009960.

* Martin, D.P., Varsani, A., Roumagnac, P., Botha, G., Maslamoney, S., Schwab, T., et al., 2021. RDP5: a computer program for analyzing recombination in, and removing signals of recombination from, nucleotide sequence datasets. Virus Evol. 7 (1), veaa087.

* Moshiri, N., 2021. ViralMSA: massively scalable reference-guided multiple sequence alignment of viral genomes. Bioinformatics 37 (5), 714-716.

* Müller, N.F., Kistler, K.E., Bedford, T., 2022. A Bayesian approach to infer recombination patterns in coronaviruses. Nat. Commun. 13 (1), 4186.

* Murray, G.G.R., Wang, F., Harrison, E.M., Paterson, G.K., Mather, A.E., Harris, S.R., et al., 2016. The effect of genetic structure on molecular dating and tests for temporal signal. Gilbert M, editor. Methods Ecol. Evol. 7 (1), 80-89.
Nguyen, L.T., Schmidt, H.A., von Haeseler, A., Minh, B.Q., 2015. IQ-TREE: a fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies. Mol. Biol. Evol. 32 (1), 268–274.

Rabadán, P., 2025. Nextstrain for Watermelon Mosaic Virus [Internet] [cited 2025 Feb 21]. Available from. https://nextstrain.org/community/PiR92/WMV.

Rabadán, M.P., Gómez, P., 2023. Global phylodynamics of two relevant aphid-transmitted viruses in cucurbit crops: cucurbit aphid-borne yellows virus and watermelon mosaic virus. Phytopathol. Res. 5 (1), 53.

Rambaut, A., Lam, T.T.Y., Max Carvalho, L., Pybus, O.G., 2016. Exploring the temporal structure of heterochronous sequences using TempEst (formerly path-O-gen). Virus Evol. 2 (1), vew007.

Ramsden, C., Holmes, E.C., Charleston, M.A., 2008. Hantavirus evolution in relation to its rodent and insectivore hosts: no evidence for Codivergence. Mol. Biol. Evol. 26 (1), 143–153.

Rasmussen, D.A., Guo, F., 2023. Espalier: efficient tree reconciliation and ancestral recombination graphs reconstruction using maximum agreement forests. Syst. Biol. 72 (5), 1154–1170.

Rasmussen, M.D., Hubisz, M.J., Gronau, I., Siepel, A., 2014. Genome-wide inference of ancestral recombination graphs. PLoS Genet. 10 (5).

Roques, L., Desbiez, C., Berthier, K., Soubeyrand, S., Walker, E., Klein, E.K., et al., 2021. Emerging strains of watermelon mosaic virus in southeastern France: model-based estimation of the dates and places of introduction. Sci. Rep. 11 (1), 1–11.

Vaughan, T.G., Welch, D., Drummond, A.J., Biggs, P.J., George, T., French, N.P., 2017. Inferring ancestral recombination graphs from bacterial genomic data. Genetics 205 (2), 857–870.