# Contrasting phylogeography of two Western Palaearctic fish parasites despite similar life cycles.

Marie-Jeanne Perrot-Minnot, Marta Špakulová, Rémi Wattier, Petr Kotlík, Serdar Düşen, Ali Aydoğdu, Christelle Tougard

► To cite this version:

Marie-Jeanne Perrot-Minnot, Marta Špakulová, Rémi Wattier, Petr Kotlík, Serdar Düşen, et al. Contrasting phylogeography of two Western Palaearctic fish parasites despite similar life cycles. Journal of Biogeography, 2018, 45 (1), pp.101-115. 10.1111/jbi.13118. hal-01683001

HAL Id: hal-01683001

https://hal.science/hal-01683001v1

Submitted on 11 Mar 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Contrasting phylogeography of two Western Palaearctic fish
parasites despite similar life cycles

Marie-Jeanne Perrot-Minnot¹*, Marta Špakulová², Rémi Wattier¹, Petr Kotlík³, Serdar Düşen⁴,
Ali Aydoğdu⁵, Christelle Tougard⁶

1 Univ. Bourgogne Franche-Comté, UMR CNRS 6282 Biogéosciences, 6 bld. Gabriel, 21000 Dijon, France

2 Institute of Parasitology, Slovak Academy of Sciences, Hlinkova 3, 04001 Košice, Slovakia.

3 Institute of Animal Physiology and Genetics, The Czech Academy of Sciences, Laboratory of Molecular Ecology, Rumburská 89, 277 21 Liběchov, Czech Republic.

4 Pamukkale University, Department of Biology, Kinikli Campus, Kinikli, 20017 Denizli-Turkey

5 Uludag University, Faculty of Veterinary Medicine, Department of Aquatic Animal Diseases, TR-16059 Gorukle - Turkey

6 Institut des Sciences de l'Evolution de Montpellier, UMR CNRS/UM/EPHE 5554, IRD 226, CIRAD 117, Université de Montpellier, Place Eugène Bataillon, CC065, 34095 Montpellier cedex 05, France.

<footnote>* Author for correspondence</footnote>

<footnote>Short running head: Phylogeography of acanthocephalan parasites in the Western Palaearctic</footnote>

<footnote>Word count inclusive of abstract, main text and references: 7534</footnote>

<footnote>Estimated number of journal pages required by the figures and tables: 4.5</footnote>
# 1 ABSTRACT

2

# 3 Aim

4 We used comparative phylogeography of two intestinal parasites of freshwater fish to test
5 whether similarity in life cycle translates into concordant phylogeographic history. The thorny-
6 headed worms *Pomphorhynchus laevis* and *P. tereticollis* (Acanthocephala) were formerly
7 considered as a single species with a broad geographical and host range within the Western
8 Palaearctic.

9 **Location** Central and eastern parts of Northern Mediterranean area, Western and Central
10 Europe, Ponto-Caspian Europe.

11 **Methods** A mitochondrial marker (COI) was sequenced for 111 *P. laevis* and 50 *P. tereticollis*
12 individuals and nuclear ITS1 and ITS2 sequences were obtained for 37 *P. laevis* and 21 *P.
13 tereticollis*. Genetic divergence, phylogenetic relationships and divergence time were estimated
14 for various lineages within each species, and their phylogeographic patterns were compared to
15 known palaeogeographic events in Western Palaearctic. Biogeographic histories of each species
16 were inferred.

17 **Results** The two species show very different phylogeographic patterns. Five lineages were
18 identified in *P. laevis*, partially matching several major biogeographical regions defined in the
19 European riverine fish fauna. The early stages of *P. laevis* diversification occurred in the peri-
20 Mediterranean area, during the Late Miocene. Subsequent expansion across Western Europe
21 and Russia was shaped by dispersal and vicariant events, from Middle Pliocene to Middle
22 Pleistocene. By contrast, *P. tereticollis* has differentiated more recently within the Western and
23 Central parts of Europe, and shows weak geographical and genetic structuring.

24 **Conclusion.** Our study highlights weak to moderate similarity in the phylogeographic pattern
25 of these acanthocephalan parasites compared to their amphipod and fish hosts. The observed
differences in the timing of dispersion and migration routes taken may reflect the use of a range
of final hosts with different ecologies and dispersal capabilities. By using a group
underrepresented in phylogeographic studies, our study is a valuable contribution to revealing
the biogeography of host-parasite interactions in continental freshwaters.

**Keywords** : amphipod, British islands, comparative phylogeography, Cyprinidae, Danube, helminth, Mediterranean, Messinian Salinity Crisis, *Pomphorhynchus*, Ponto-Caspian.

33
34 INTRODUCTION

35

36 Comparative phylogeography is a powerful approach to unravel the effects of biogeography
37 and evolutionary history on the genetic structure of co-distributed species. More specifically,
38 the life history and ecology of species provide a comparative framework for testing whether
39 their phylogeographic structure is more impacted by palaeogeographic and palaeoclimatic
40 events, or by biotic factors (Papadopoulou & Knowles, 2016). This is particularly exemplified
41 in parasitic taxa, where differences in life history traits such as transmission mode, life cycle,
42 and host specificity can result in discordant genetic structure in response to similar historical
43 events (van der Mescht et al., 2015, and references therein). Yet, comparative phylogeographic
44 studies on parasitic species are relatively scarce when considering that metazoan parasites
45 comprise approximately 40% of all living metazoan species (Dobson et al., 2008). Most of
46 these studies investigate host-parasite co-phylogeny, and predict weak host-parasite co-
47 evolutionary signal if host associations play little role in the dispersal and diversification of
48 parasites (Hoberg & Brooks, 2008) as, for instance, in parasites with long-lived free-living
49 stages, in generalist multi-host parasites, or in parasites with a complex life cycle involving
50 intermediate host(s). By contrast, few studies have compared the phylogeographic patterns of
51 multiple parasite species having similar life cycle.

52 The aim of our study is to compare the phylogeographies of two Western Palaearctic
53 freshwater parasite species sharing a similar life cycle and ecology, and occupying overlapping
54 ranges. Continental-scale biogeography studies have demonstrated that the palaeogeographic
55 history of continental waters within the Western Palaearctic is complex, with numerous
56 dispersal and diversification opportunities for freshwater and brackish-water species,
57 determined by the alternations of isolation and interconnection among hydrographic basins
58 (several references in Costedoat et al., 2006; Levy et al., 2009; Dubut et al., 2012; Appendix
S1). Several regions (e.g. the peri-Mediterranean, Aegean, and Ponto-Caspian regions, and the Danubian system) have been identified as important sources of diversification during a time frame spanning from the Oligocene-Miocene to the Quaternary post-glacial period (see Appendix S1). By contrast, macroparasites of fish and crustaceans have been used only exceptionally for historical biogeography reconstructions.

We address whether two species of parasites with complex life cycles, *Pomphorhynchus laevis* (Müller, 1776) and *P. tereticollis* (Rudolphi, 1809) (thorny-headed worms, Acanthocephala), responded in a similar way to the palaeogeography of continental waters within the Western Palaearctic. Because of previous taxonomic confusion (Appendix S1 in Supporting Information), *P. tereticollis* has been recorded as *P. laevis* in most parasitological surveys (Špakulová et al., 2011; Vardić-Smrzlić et al., 2015), therefore leading us to use here the general term *P. laevis s.l.* wherever taxonomic distinction between these species could not be made. *Pomphorhynchus laevis s.l.* is potentially a valuable parasite model in the phylogeography of continental waters because it exhibits a large geographical distribution, and uses a broad range of freshwater and brackish-water fish species as final hosts and gammaridean amphipods as intermediate hosts (Kennedy, 2006; Špakulová et al., 2011; Vardić-Smrzlić et al., 2015)-(Appendix S1). Although the precise distribution and host range of each of the two *Pomphorhynchus* species is unclear, they appear to occur in sympatry across Europe, with verified records from France and Slovakia, where they also share the same crustacean hosts and some of the fish hosts (Špakulová et al., 2011; Perrot-Minnot, unpublished data).

We expected that the expansion and differentiation of their intermediate amphipod host and final fish host, which began as early as the Late Eocene-Middle Miocene (Kotlík & Berrebi, 2002; Levy et al., 2009; Perea et al., 2010; Hou et al., 2011; Dubut et al., 2012; Buonerba et al., 2015; Mamos et al., 2016; Appendix S1), has provided ecological opportunities for their own dispersal and differentiation. Our study relies on a range-wide sampling both in terms of
hosts and geography, to identify evolutionary lineages for *P. laevis* and *P. tereticollis*, and to
compare their respective biogeographic histories. Mitochondrial cytochrome *c* oxidase I (COI)
and internal transcribed spacers (ITS1 and ITS2) sequences were used for species assignment
and estimation of divergence time, while only COI was used to infer intraspecific phylogenetic
relationships between lineages, and plausible biogeographic scenarios accounting for the
genetic diversity within each species. We compared the observed phylogeographic pattern of
each species to the major biogeographical units identified in European riverine fish fauna based
on a cluster analysis of fish community composition (Reyjol *et al.*, 2007). We also related these
patterns to known geological events and palaeoenvironmental conditions in Western
Palaearctic, starting from the Messinian Salinity Crisis (MSC) in the Mediterranean around 6
Myr ago. At this time, near-desiccation of the Mediterranean Sea promoted the cutting of deep
fluvial canyons, connections of drainage basins, and lake or river captures (Orszag-Serber,
2006; Rouchy & Caruso, 2006), thereby offering dispersal opportunities for freshwater and
brackish-water species. The phylogeographic patterns uncovered were also related to the
expansion of freshwater systems during Pliocene (Neubauer *et al.*, 2015), and the major
Quaternary glaciations during Early and Middle Pleistocene (1.7 to 0.12 Myr).

We tested the hypothesis that these co-occurring species with similar life cycle and
ecologies shared similar biogeographic histories. On the other hand, contrasting
phylogeographic patterns may reveal ecological divergence or historical contingency, and
would call for a thorough assessment of their respective host specificities. By comparing the
phylogenetic pattern of *P. laevis* and *P. tereticollis* to that of their gammaridean and fish hosts,
we address the relative contribution of the intermediate and final hosts to their dispersal and
diversification. Thus, our study makes a contribution to revealing the historical biogeography
of host-parasite interactions in continental freshwaters of the Western Palaearctic.
109 MATERIALS AND METHODS

110

111 Data acquisition

112 One hundred and sixty one *P. laevis s.l.* individuals were collected from 56 localities across
113 Western Palaearctic (Fig. 1; Table S1.1.), and fixed and stored in 95% ethanol. DNA was
114 extracted following a standard CTAB-proteinase K and phenol-chloroform protocol (Perrot-
115 Minnot, 2004), from either a small piece of adult worm or from whole larval cystacanth. A
116 portion of the COI and the ITS + 5.8S genes were PCR-amplified and purified according to
117 Perrot-Minnot (2004) (Table S1.2.). COI sequences and ITS + 5.8S sequences were obtained
118 by direct sequencing of purified PCR products (MWG, Germany, or Macrogen Inc., Seoul,
119 South Korea and Amsterdam, The Netherlands). Pseudogenes were occasionally amplified with
120 the COI general primers in *P. laevis* samples, as confirmed by the cloning of PCR products that
121 yielded multiple sequences. We therefore designed internal primers specific to the authentic
122 mtDNA COI gene, and used them to analyse some *P. laevis* samples from the Balkans and Italy
123 (Table S1.2).

124 The majority of amplicons of both genes were sequenced in both directions, and were
125 manually aligned with BIOEDIT software (Hall, 1999). Seven of the previously collected
126 sequences of *P. tereticollis* used for the taxonomic revision of this species (Špakulová *et al.*,
127 2011) have been included in the dataset (accession numbers: ITS, JF06705; COI; JF706706,
128 JN695504, JN695505 to JN695508).

129

130 Phylogenetic inference

131 Phylogenetic reconstructions were performed using two probabilistic methods: a maximum
132 likelihood method (ML) implemented in PHYML 3.0 (Guindon et al., 2010) and a Bayesian
133 approach (BA) in MRBAYES 3.1.2 (Ronquist & Huelsenbeck, 2003). Best-fitting models of
sequence evolution for each marker were determined using MRMODELTEST 2.3 (Nylander, 2004). ML and BA analyses were conducted under the general time-reversible (GTR) model (Yang, 1994) with a proportion of invariable sites (I) and a gamma distribution (G) for COI, and under the HKY model (Hasegawa *et al.*, 1985) +I+G for ITS.

Node robustness was estimated with ML bootstrap percentages (BP) after 1000 pseudo-replicates, whereas Bayesian posterior probabilities (PP) were obtained from 50% majority rule consensus trees, after discarding the first 25,000 trees as burn-in. Three independent runs of five incrementally heated Markov chains Monte Carlo (MCMC) were performed with trees sampled every 100th generation for 5,000,000 generations.

143

## 144 Species tree and divergence time estimation

The species tree for *P. laevis* and *P. tereticollis* and for geographically delimited populations within each species (Fig. 1) was estimated using MCMC for the multispecies coalescent model implemented in STARBEAST2 (Ogilvie *et al.*, 2017). The method infers the species tree from multiple genes sampled from multiple individuals from each species (Ogilvie *et al.*, 2017). The COI (158 sequences) and ITS (52 sequences) genes were therefore included in the analysis as independent loci. Because STARBEAST2 makes the assumption that there is no gene flow following the population divergence, COI haplotype Pl46, likely representing a recent migration from the Adriatic to Tyrrhenian population (Fig. 2a), was excluded from the input datasets.

Six models corresponding to three molecular clock assumptions (strict, relaxed lognormal or relaxed exponential) and two priors for the species tree (Yule and Birth-Death), within which the gene trees follow the multispecies coalescent prior, were compared using the Akaike's information criterion through MCMC (AICM; Raftery *et al.*, 2007) in order to determine which combination best fits our data. Each analysis was repeated with the same priors

158
but with a different random seed to ensure MCMC convergence and thus robustness of the analysis. Five independent runs of 100,000,000 generations each were then performed with the best-fitting clock and speciation model combination, with a burn-in stage of 10%. Data from repeated runs were combined with LOGCOMBINER 2.4.0 (Bouckaert et al., 2014) and a maximum clade credibility consensus tree with mean node heights was generated using TREEANNOTATOR v2.4.0 (Bouckaert et al., 2014).

The absence of acanthocephalan (or even helminth) fossils complicates construction of a time-calibrated tree for *P. laevis* s.l. Furthermore, significant rate variation among invertebrate phyla limits the reliability of an arthropod COI clock (1.4% to 2.8% sequence divergence/Myr) (Knowlton & Weigt 1998; Solà *et al.*, 2013 and ref. therein), for inferring divergence times within another taxa (Thomas *et al.*, 2006). Alternatively, palaeogeographic events of known age can be used to calibrate a molecular clock (Solà *et al.*, 2013). Preliminary phylogenetic analysis in *P. laevis* revealed a basal split between Italian lineages and all other lineages. One of the major palaeogeographic events in the Mediterranean basin was the MSC and its associated changes in the extent and connections of freshwater bodies. The date of 6.05 ± 0.09 Myr (6.14–5.96 Myr, Rouchy & Caruso, 2006) for the isolation of the Mediterranean Sea and the divergence of the most recently diverged Italian population of *P. laevis* (Adriatic; Pl_L4) was therefore used as the calibration point.

## Genetic characterization of *Pomphorhynchus laevis* and *P. tereticollis* populations

Due to the limited intraspecific variation of the ITS sequences these analyses were performed with the COI dataset only. A Kimura-2-parameter distance was chosen for the estimation of the genetic divergence within and between *P. laevis* and *P. tereticollis* lineages using MEGA v5.2.2 (Tamura *et al.*, 2011). The haplotype number (*nh*), nucleotide (*π*) and haplotype (*h*) diversities and the average number of nucleotide differences between two sequences (*k*) were obtained for
each lineage using DNASP 5.10.01b (Librado & Rozas, 2009). The demographic history of lineages (stability or expansion) was investigated using three neutrality tests and mismatch indices (see Appendix S2 for details). Finally, population structures were evaluated with a median-joining network with NETWORK v4.5.1.6 (http://www.fluxus-engineering.com/sharenet.htm; Bandelt *et al.*, 1999) for the haplotypes of *P. laevis* and *P. tereticollis*.

## Biogeographic analysis

To infer the geographical origin of *P. laevis* and *P. tereticollis*, the distribution of both species were subdivided into, respectively, 13 and five geographical areas, corresponding to the main drainages and/or recipient seas (Table S1.1.). Considering the complex history of the Danubian system and its importance in most phylogeographic studies, four sub-areas were distinguished: Carpathian, Pannonian, and Balkan areas, and the Danube itself. Biogeographic inferences were obtained with the programs Statistical Dispersal-Vicariance Analysis (S-DIVA; Yu *et al.*, 2010, 2014) and Bayesian Binary MCMC (BBM) as implemented in the RECONSTRUCT ANCESTRAL STATE IN PHYLOGENIES software (RASP 3.0; Yu *et al.*, 2014, 2015). S-DIVA analyses were performed with Bayesian trees obtained for phylogenetic reconstructions. The maximum number of ancestral areas was set to two. In BBM analyses, MCMC was carried out using 10 chains (temperature 0.1) for 100,000 cycles with state sampled every 100th cycle, and a burn-in stage of 100. The Jukes-Cantor model +G was used, while the maximum number of ancestral areas and the root distribution were set, respectively, to 2 and as null. In both cases, the possible ancestral ranges were thus obtained at each node of S-DIVA and BBM trees.

## RESULT
209 Phylogenetic inference

210 Among the 161 sequences of *P. laevis* s.l., 111 belong to *P. laevis* and 50 to *P. tereticollis*,
211 based on the interspecific differences at ITS and COI reported in Špakulová *et al.* (2011). The
212 gene tree constructed for 58 concatenated sequences of ITS1 and ITS2 confirmed the reciprocal
213 monophyly of *P. tereticollis* and *P. laevis* (Appendix S2, Fig. S2.1.).

214 The COI alignment represents 614 positions for both *P. laevis* and *P. tereticollis* with,
215 respectively, 172 and 25 phylogenetically informative sites. For both species, BA analysis (Fig.
216 2a and 2b) and ML (Fig. S2.2) provided congruent tree topologies. When included in the same
217 analysis (results not shown), COI sequences of *P. laevis* and *P. tereticollis* present a reciprocal
218 monophyly with high support (BP=100%; PP= 1.00). In *P. laevis*, five highly supported clades
219 (Pl_L1 to Pl_L5) are identified (Fig. 2a and S2.2a), suggesting strong geographic structuring of
220 the populations. Clades Pl_L4 and Pl_L5 are two basal lineages found in Italian populations,
221 one with Tyrrhenian (Pl_L5), and one with Tyrrheno-Adriatic (Pl_L4) distribution. Clade
222 Pl_L3 is found in the Aegean and south-western Anatolian parts of Turkey. Clade Pl_L2 is
223 distributed in Ponto-Caspian Europe, and includes haplotypes from the Volga (Pl_L2b) and
224 Danube and Vistula rivers (Pl_L2a and c). Finally, clade Pl_L1 has a wide geographic
225 distribution covering three major biogeographic areas: Eastern peri-Mediterranean, Ponto-
226 Caspian, and Western Europe. It comprises haplotypes from the Balkan tributaries of the
227 Danube (Pl_L1b), from an Aegean river (Pl_L1c), from Pannonian-Western Carpathians and
228 Western Europe (Pl_L1a), and from Eastern Carpathians (Pl_L1b, Pl_L1d) (Fig. 2a). By
229 contrast, the COI gene tree of *P. tereticollis* shows much shallower structure with only one
230 well-supported clade (Pt_L1) including haplotypes from Western and Central Europe, while
231 the remaining haplotypes (collectively labelled as Pt_L2) are distributed in Western and Ponto-
232 Caspian Europe (Fig. 2b and S2.2b).

233
234 **Molecular dating**

235 AICM values suggested that the relaxed exponential clock and the Yule model were
236 significantly more likely given our COI and ITS datasets than the other priors. Consequently,
237 the time-calibrated species tree for *P. laevis* and *P. tereticollis* was estimated with these models
238 (Fig. 3). With TRACER v.1.6 (Rambaut *et al.*, 2014), Markov chain convergence was ascertained
239 by visual inspection and comparison of the traces and posterior distributions, and the effective
240 sample size was >200 for all parameters.

241 The species tree (Fig. 3) is largely concordant with the COI gene tree estimated by ML
242 and BA (Figs 1 and 2). The calibration point gave an estimate of 0.0201 (0.0199 – 0.0202)
243 substitution/bp/Myr, similar to the substitution rate calibrated for invertebrate mitochondrial
244 genes (1.4 to 2.8%: Knowlton & Weigt, 1998; Solà *et al.*, 2013). The analyses suggested that
245 the split between *P. laevis* and *P. tereticollis* took place during the Late Miocene (8.40 Myr).
246 The major geographic populations of *P. laevis* began to diverge during Late Miocene (from
247 6.82 Myr). However, the divergence within groups of populations characterised by the two
248 major COI clades (Pl_L1 and Pl_L2) seems related to Middle (0.50 Myr) and Early (1.27 Myr)
249 Pleistocene, respectively (Fig. 3). By contrast, the populations of *P. tereticollis* appear to have
250 diverged much later, in Middle Pleistocene (0.42 Myr) (Fig. 3).

251

252 **Genetic diversity, demographic history and population structure**

253 In *P. laevis*, percentage divergence is from 1.2 ± 0.3% to 2.4 ± 0.4% within clades, and from
254 10.5 ±1.2% to 20.3 ± 2.1% between clades (Table 1a). The genetic divergence is 1.6 ± 0.3%
255 within *P. tereticollis*, and reaches 24.4 ± 1.8% between the two species (Table 1b). For both
256 species, the nucleotide (π) and haplotype (h) diversities are rather homogeneous, whereas the
257 average number of nucleotide differences between two sequences is relatively high (Table 2).
258 The three neutrality tests and mismatch distribution reject population expansion for the five *P*.
259 *laevis* lineages and for *P. tereticollis* (Table S2.3.).

260 For both species, the median-joining network showed the same general pattern as the
gene tree, with five *P. laevis* and two *P. tereticollis* haplogroups separated by 125 mutational
261
262 steps (Fig. 4). *Pomphorhynchus laevis* groups are separated by 47 to 83 mutational steps, while
263 the two *P. tereticollis* groups are separated by only 7 mutational steps (Fig. 4).

264 Although the demographic analysis performed for each lineage of *P. laevis* rejects
265 population expansion, the star-like network and the polytomic tree topology of the Danubian
266 sub-lineage Pl_L2a suggest population expansion (Fig. 4). By contrast, the three basal peri-
267 Mediterranean lineages in Italy and Turkey, and the derived Volga sub-lineage could have
268 experienced a loss of diversity, as suggested by the number of mutational steps and missing
269 haplotypes within each lineage (Fig. 4).

270

271 **Biogeographic analyses**

272 For *P. laevis*, the geographic origin of ancestral populations is poorly resolved both in S-DIVA
273 and BBM analyses, as evidenced in widespread distribution range and occurrence frequencies
274 less than 10% (nodes 120, 115, 108 in Fig. S2.3). Both analyses generally suggest, however,
275 that the ancestral range encompassed the peri-Mediterranean area (BBM) or the peri-
276 Mediterranean to Ponto-Caspian area (S-DIVA). The ancestral population of the Turkish and
277 all non-Mediterranean populations was distributed in several alternative biogeographic areas,
278 all anchored in the Eastern peri-Mediterranean (node 108 in Fig. S2.3).

279 The ancestral population of all non-Mediterranean lineages could have been in the
280 Danubian-Carpathian area within Ponto-Caspian Europe (EH in Table 3; node 103 in Fig. S2.3).
281 This population gave origin to a Danubian population, ancestral to the Ponto-Caspian lineage
282 (Pl_L2), and a Carpathian-Western Europe population, ancestral to the widely distributed
283 Ponto-Caspian and Western Europe lineage (Pl_L1) (H & EC in Table 3; nodes 81-80 and 102
284 in Fig. S2.3). Genetic differentiation within the former occurred in the Danubian–Caspian Sea
285 area, resulting in one sub-lineage in the Danube and Vistula rivers (Pl_L2a), and the other
286 strictly Caspian (Volga, Pl_L2b), with no evidence of subsequent mixing (Fig. 2a and S2.3).
287 Genetic differentiation within the Carpathian–Western Europe ancestral population occurred in
288 a widespread area, and resulted in an Eastern Carpathian sub-lineage (Pl_L1d), a Western
289 Carpathian-Pannonian-Western European sub-lineage (Pl_L1a) and an Eastern Carpathian-
290 Northern Aegean sub-lineage (Pl_L1b, c) (Fig. 2a; Table 3; nodes 102 and 91 in Fig. S2.3). The
291 latter sub-lineage experienced at least two biogeographic events: the splitting of a Northern
292 Aegean lineage (Pl_L1c), and the dispersal towards the Balkans (Pl_L1b) (Table 3; nodes 87
293 and 85 in Fig. S2.3, respectively). The Pannonian-Western European sub-lineage further
dispersed towards its western margins, to the Atlantic river drainage (Pl_L1a; Fig. S2.3).
294

295 For *P. tereticollis*, S-DIVA and BBM analyses suggest one possible ancestral range in
296 Western Europe, more specifically in the Rhone drainage basin, possibly reaching as far as the
297 Ponto-Caspian area (Table 3; node 59 in Fig. S2.4). The geographic expansion of this ancestral
298 population resulted in the formation of a Central and Western European lineage (Pt_L1), and a
299 Western and Ponto-Caspian European lineage (Pt_L2) (node 56 in Fig. S2.4). The Central and
300 Western European lineage subsequently underwent repeated dispersal, vicariance, and
301 admixture events between the Baltic Sea and the British Isles (Fig. 2b and S2.4). The Western
302 and Ponto-Caspian European lineage also experienced a complex history of dispersal,
303 vicariance, admixture and extinction events within an area comprising the Rhone and Rhine
304 systems, the Carpathians, and the British Isles (Table 3; nodes 55 and 50 in Fig. S2.4).
305

306

306 DISCUSSION

307
The present study aimed at comparing the range-wide phylogeographical structure of two Western Palaearctic acanthocephalan parasites with a similar complex life cycle. Our results supports the evolutionary relationships between *P. laevis* and *P. tereticollis* inferred from species-specific features (Špakulová *et al.*, 2011). Phylogenetic and biogeographic analysis, and molecular dating using the Messinian Salinity Crisis (6 Myr) to age-calibrate the phylogeny, revealed contrasted biogeographic histories between these species. Their distinct differentiation pattern at both geographic and time scales highlights their potential contribution to our understanding of the Western Palaearctic historical biogeography, as detailed below.

*The origin and comparative biogeographic structuring of P. laevis and P. tereticollis*

*Pomphorhynhus laevis* and *P. tereticollis* diverged during the Late Miocene from a common ancestral lineage, possibly within the peri-Mediterranean area. Despite the low resolution of biogeographic analysis at deep nodes, three lines of evidence support this geographic origin: (i) in the haplotypic network, the closest relative of *P. tereticollis* to *P. laevis* is from the Rhone drainage, and is related to Turkish *P. laevis* haplotypes, (ii) in *P. tereticollis*, the highest number of haplotypes is found in the Rhone drainage, and (iii) the basal lineages of *P. laevis* are found in the Central and Eastern Mediterranean area.

In *P. laevis*, five lineages of Pleistocene age with 10.5% to 20.3% divergence were identified. They are distributed in the Central and Eastern peri-Mediterranean and throughout most of the Ponto-Caspian and Western Europe, but with almost completely non-overlapping geographical distribution except in the Danubian system. Such high divergence and phylogeographic structure provide evidence for multiple populations with a history of long-term persistence and separation. By contrast, *P. tereticollis* has no clear biogeographic structuring, as evidenced in a basal polytomy restricted mainly to Western, Central, and Ponto-
Caspian Europe. The overall genetic diversity of *P. tereticollis* is comparable to the within-lineage genetic diversity of *P. laevis*.

*In and Out of the Mediterranean Sea: the pre-Quaternary evolutionary history of P. laevis*

The two oldest genetic lineages of *P. laevis* find their origin in Late Miocene in the Italian Peninsula, one strictly Tyrrhenian and one Tyrrheno-Adriatic. According to a first scenario of colonization of Southern Europe (Tsigenopoulos *et al.*, 2003; Levy *et al.*, 2009; Gante, 2011), dispersal from North Africa during the Messinian lacustrine phase of the Mediterranean Sea could account for the colonization of the Western Mediterranean. Under this 'southern Sea dispersal scenario', the Italian peninsula would have been colonized during the Lago Mare stage (6 Myr). However, the oldest lineage to all *P. laevis* has here a pre-Messinian origin, and the two basal Italian lineages are too divergent to both originate from this short-time event. Alternatively, the Italian peninsula could have been colonized prior to the intensification of Alpine orogenesis during Miocene, through river capture across north-eastern Europe (Perdices *et al.*, 2003; Perea *et al.*, 2010). The Adriatic and Tyrrhenian populations could have subsequently diverged from a widespread ancestral population, at a time coincident with the emergence of the Northern Apennine Chain in Late Miocene; the flowing of Tyrrhenian and Adriatic rivers through opposite side of the Apennines would have precluded later exchange, as suggested for the speciation of endemic fluvio-lacustrine *Barbus* species (Buonerba *et al.*, 2015). Our data, however, do not allow us to favour one scenario over the other, and both a 'southern route' and a 'northern route' have been proposed for the colonization of the Mediterranean by leuciscine freshwater fish (Durand *et al.*, 2003; Perdices *et al.*, 2003; Levy *et al.*, 2009; Perea *et al.*, 2010; Gante, 2011).

The third distinct Mediterranean lineage is specific to the Turkish populations of the
358 North Aegean-West Anatolian area, and is basal in the phylogenetic tree relative to the non-359 Mediterranean lineages. Its geographic origin is unclear given insufficient resolution of the 360 biogeographic analysis. However, the Lago Mare stage of the MSC has possibly allowed the 361 dispersal of freshwater hosts of *P. laevis* within the eastern part of the Mediterranean basin. The 362 presence of contiguous freshwater bodies at a regional scale, suspected from biogeographic 363 studies and faunal assemblages (Solá *et al.*, 2013; Neubauer *et al.*, 2015), could have promoted 364 the dispersal and subsequent isolation of freshwater populations between adjacent areas 365 (Tsigenopoulos *et al.*, 2003; Perea *et al.*, 2010; Levy *et al.*, 2009), or at a circum-Mediterranean 366 scale (Durand *et al.*, 2003; Gante, 2011; Buonerba *et al.*, 2015). Additionally, conditions for 367 dispersal for freshwater organisms were less restricted within the Eastern Mediterranean than 368 in the Western Mediterranean (Orszag-Sperber, 2006; Levy *et al.*, 2009). For instance, several 369 cyprinid fish species have dispersed in the Eastern Mediterranean during the Lago Mare 370 between 5.5 and 5.3 Myr ago, and subsequently differentiated after the refilling of the 371 Mediterranean Sea (Durand *et al.*, 2003; Dubut *et al.*, 2012). This is consistent with the 372 molecular dating obtained for the split between the Turkish population and the non-373 Mediterranean populations around 3.17 Myr, that postdates the Lago Mare event. 374

374 The Late Miocene marine flooding of the Mediterranean Basin and the final settling of 375 geomorphological barriers during Miocene (Alps in Northern Italy) and Pliocene (Turkish 376 straits and Marmara Sea) have precluded any subsequent peri-Mediterranean dispersal, and 377 promoted the independent evolution of these three populations. Their long-term persistence in 378 isolation is reflected in their low contribution to the recent genetic pool of *P. laevis* in Western 379 and Ponto-Caspian Europe. This result is in line with the high level of endemism of the 380 Mediterranean freshwater ichthyofauna (Reyjol *et al.*, 2007; Perea *et al.*, 2010; Gante 2011). 381 To confirm this pattern, more samples need to be analysed from other peri-Mediterranean 382 localities. During Mid-Pliocene, the expansion and differentiation of one or several populations
ancestral to the two non-Mediterranean lineages matches the expansion of landmasses and
associated networks of new freshwater lacustrine-riverine habitats (Neubauer *et al.*, 2015).

*Pleistocene expansion and differentiation of P. laevis and P. tereticollis throughout Europe*

As for many Western Palaearctic taxa during the Pleistocene, *P. laevis* and *P. tereticollis* have
experienced postglacial range expansion alternating with more or less severe range contraction
and fragmentation. In *P. laevis*, the Danubian-Caspian populations (Pl_L2) began to diversify
during the Pleistocene around 1.27 Myr ago, possibly within the northern Black Sea-Caspian
Sea region. This scenario agrees with the existence of Ponto-Caspian refugium in several
freshwater fishes (e.g. Durand *et al.*, 1999; Perdices *et al.*, 2003; Kotlík *et al.* 2004; Seifertová
*et al.*, 2012). The pattern of Pleistocene differentiation between the Western and Ponto-Caspian
European populations (Pl_L1) appears to be geographically more widespread. Our molecular
dating and ancestral biogeographic analysis suggest several expansion and differentiation
events around 0.50 Myr ago from a Carpathian-Western Europe ancestral population, south-
eastward into several sub-lineages covering the Eastern Carpathian, Balkan and North Aegean
areas, and westward into a Pannonian-Western Europe area. Given the history of these
Pleistocene lineages, the occurrence within the Danubian system of two differentiated
populations (Pl_L2 and Pl_L1), in the Danube itself and in its tributaries respectively, is
puzzling. One scenario that remains to be tested is the replacement of Pl_L1 in the Danube by
Pl_L2 from the Vistula or Volga Rivers along with introduced Ponto-Caspian amphipods and/or
fish hosts. Replacement followed by rapid population expansion could account for the star-like
network and polytomic gene tree of this Danubian population.

For *P. tereticollis*, the haplotype network, molecular dating and biogeographic analysis
suggest an ancestral population in the Rhone basin, which began to diversify around 0.42 Myr
ago. The maintenance of ancestral polymorphism in the Rhone drainage is evidenced in the highest haplotypic diversity found in its tributaries at a very local scale. One population descending from the Rhone ancestral population appears to have dispersed northward to the North and Baltic Sea, and the other one eastward to the Rhine and the Carpathians. Dispersal from a Rhone ancestral population to Rhine river could have occurred when the contemporary Alpine Rhine was diverted northward to flow into the Rhine drainage during Middle Pleistocene (Dubut *et al.*, 2012, and ref. therein). Occasional inter-basin exchanges within the Rhine/Rhone corridor prevented Pleistocene differentiation between these basins. On the other hand, considering repeated and transient connections between the Rhine, Rhone and Danube drainages during Pleistocene (op. cit.), the restricted Ponto-Caspian distribution of *P. tereticollis* to the Carpathians within this lineage is puzzling, and raises issues about either sample representativeness, ecological constraints, or Pleistocene extinction in Ponto-Caspian Europe.

Pomphorhynchus tereticollis and the biogeographic history of the British Isles

In *P. tereticollis*, several waves of Pleistocene colonization from continental Europe left their genetic imprint in the British Isles, with the presence of the two COI lineages Pt_L1 and Pt_L2. Two routes of dispersal from continental Europe accounted for the introduction of haplotypes of both lineages into the British Isles during the past 0.42 Myr, a northern route across the North Sea, and a southern route through the English Channel. The northern route across the Baltic-North Seas allowed the introduction of Pt_L1 haplotypes. The lack of geographic structure within Pt_L1 suggests a recent colonization or recurrent dispersal opportunities since the Middle Pleistocene between the Northern British Isles and Northwestern Europe, despite prolonged periods of glaciation. At the time, the Fennoscandinavian ice sheet (FIS) could have
433

forced northern rivers from the current Baltic Sea drainage to flow into the North Sea and
connect with Scottish continental waters. The southern route through the English Channel
allowed the introduction of Pt_L2 haplotypes from continental Rhine-Rhone populations. Here
again, the weak genetic structure within this lineage suggests that introduction into Southern
England has a recent origin or, alternatively, has occurred repeatedly since the Middle
Pleistocene. This latter interpretation is congruent with the activity of the 'Fleuve Manche'
palaeoriver over the last 0.35 Myr, connecting rivers from southern England and Western
Europe, including the Thames, Seine and Rhine, via the English Channel (Toucanne et al.,
2009). The presence of the two lineages of P. tereticollis validates the description of distinct
strains of P. laevis s.l. in the British Isles (O'Mahony et al., 2004; Kennedy, 2006), here
identified as P. tereticollis.

434

435

436

437

438

439

440

441

442

443

444

445

*Host-related biogeographic history*

446

447

Overall, the phylogeographic pattern of these parasites is not directly comparable to that of host
species, with respect to the geographic distribution and timing of diversification.

*Pomphorhynchus laevis* shows higher co-occurrence of divergent genetic lineages across
Western Europe and Ponto-Caspian area, and more recent divergence during Pleistocene,
compared to freshwater amphipod hosts (Hou *et al.*, 2011; Mamos *et al.*, 2016), and to
rheophilic fish hosts (Kotlík & Berrebi, 2002; Costedoat *et al.*, 2006; Šedivá *et al.*, 2008).

448

449

450

451

452

453

454

455

456

Compared to its contemporary main fish hosts, *P. laevis* shows more ancient differentiation
than its common host, *Barbus barbus* (Kotlík & Berrebi, 2001; Kotlík *et al.*, 2004), and a
distinct Pleistocene history compared to the European chub, *Squalius cephalus* (Seifertová *et
al.*, 2012), but with a similar biogeographic pattern (Durand *et al.*, 2000). Interestingly, lowland
warm-adapted fish species responded differently to Pleistocene glaciation cycles compared to

457
458 cold adapted species from headwaters, with higher dispersal and/or extinction events (Šedivá
459 et al., 2008; Dubut et al., 2012; Seifertová et al., 2012). The use of a range of final hosts with
460 different ecology may thus explain the mixed pattern of the timing and migration route observed
461 in *P. laevis*. The most striking difference compared to all freshwater fish studied in the Western
462 Palaearctic is the co-occurrence, within the Danubian system, of two lineages, one in the
463 Danube itself, the other in its tributaries. The ichthyofauna and amphipod fauna from the
464 Danube have been profoundly changed with the introduction of Ponto-Caspian species in recent
465 times. Replacement by the Ponto-Caspian *P. laevis* may have occurred with the expansion of
466 the barbel *B. barbus* from Black Sea refugium into the Danube around the end of the Pleistocene
467 (Kotlík et al., 2004), along with Ponto-Caspian amphipod hosts. The use of distinct intermediate
468 hosts in the Danube and its tributaries could further constitute an ecological barrier to lineage
469 mixture within the Danube basin, a hypothesis that deserves further investigation.

470 In *P. tereticollis*, host strains have been reported in the British Isles (O'Mahony *et al*.,
471 2004; Kennedy, 2006), that could match the two phylogeographic lineages identified here. The
472 Central and Western Europe lineage likely correspond to the marine-estuarine strain infecting
473 flatfish, plaice and founder and the Irish strain infecting trout, and the Western and Ponto-
474 Caspian Europe lineage to the strain found in cyprinids. We hypothesize that only the Irish
475 marine strain, parasitizing anadromous cold-adapted fish hosts with long distance dispersal,
476 could have reached the Northern part of the British Isles from the Baltic and North Seas during
477 Pleistocene. This scenario is further supported by the literature records of this strain infecting a
478 salinity-tolerant (*Gammarus duebeni*) and an estuarine (*Gammarus zadacchi*) gammarid
479 species as intermediate hosts (Guillen-Hernandez & Whitfield, 2001, and ref. herein).

480

481 *Conclusion*

482
Our phylogeographic study provides unique insights into the contrasted histories and distributions of two species of acanthocephalan parasites of fish. The phylogeographic patterns of *P. laevis* and *P. tereticollis* are concordant with known palaeogeographic events within the Western Palaearctic. They are also weakly to moderately concordant with the phylogeographic pattern of their amphipod and fish hosts, as expected for multihost parasites with a complex life cycle. A more extensive sampling scheme, with respect to host range, within-population variation, and biogeographic area, should be implemented in the future, to allow full understanding of the relative contribution of historical and ecological factors to the biogeographic history of *P. laevis* and *P. tereticollis*. The possible existence of cryptic species within *P. laevis*, as suggested by the extent of genetic divergence between lineages and divergence time estimates, should also be addressed. We expect that future research on these subjects will further enhance the value of the acanthocephalan parasites as models for studying the historical biogeography of host-parasite interactions in continental freshwaters.

496

497

## 498 Acknowledgements

The following persons made a significant contribution to the present study by kindly providing samples of *Pomphorhynchus*: Federica Berrilli and Bahram Sayyaf Dezfuli (Italy), Celia Holland (Ireland), Yuri Kvach (Ukraine), Milen Nachev and Bernd Sures (Germany), and Markéta Ondračková (Czech Republic). We would like to thank three anonymous referees for constructive comments on earlier drafts of the manuscript. The voucher specimens of *Pomphorhynchus* from Isikli lake, Denizli (Turkey), were deposited in the Biology Department, Faculty of Sciences and Arts, Pamukkale University, Denizli, Turkey (PAU-HELM-1-5/2007 and PAU-HELM-1-4/2008) (EMBL accession numbers LN994949 and LN994950). All other voucher specimens are deposited in a collection of the UMR-CNRS Biogeosciences, Université
Bourgogne Franche-Comté (Dijon, France) as EtOH-preserved worms. M.-J. P.-M. received funding support from the ANR, project PARADIV (ANR-07-BLAN-0209). P.K. received funding support from the project EXCELLENCE CZ.02.1.01/0.0/0.0/15_003/0000460 OP RDE and the institutional support RVO:67985904. Phylogenetic and molecular dating analyses were performed through the technical facilities of the Platform Montpellier Bioinformatics Biodiversity of the "Institut des Sciences de l'Evolution de Montpellier" (Centre Méditerranéen de l'Environnement et de la Biodiversité, Montpellier, France).

515

## 516 References

517

Bandelt, H.J., Forster, P. & Röhl, A. (1999). Median-joining network for inferring intraspecific phylogenies. *Molecular Biology and Evolution*, **16**, 37-48.

520

Bouckaert, R., Heled, J., Kühnert, D., Vaughan, T., Wu, C.H., Xie, D., Suchard, M.A., Rambaut, A. & Drummond, A.J. (2014) BEAST 2: A Software Platform for Bayesian Evolutionary Analysis. *PLoS Computational Biology*, **10**, e1003537.

524

Buonerba, L., Zaccara, S., Delmastro, G.B., Lorenzoni, M., Salzburger, W. & Gante, H.F. (2015). Intrinsic and extrinsic factors act at different spatial and temporal scales to shape population structure, distribution and speciation in Italian *Barbus* (Osteichthyes: Cyprinidae). *Molecular Phylogenetics and Evolution*, **89**, 115-129.

529

Costedoat, C., Chappaz, R., Barascud, B., Guillard, O. & Gilles, A. (2006). Heterogeneous colonization pattern of European Cyprinids, as highlighted by the dace complex (Teleostei: Cyprinidae). *Molecular Phylogenetics and Evolution*, **41**, 127-148.

532
533
534 Dobson, A., Lafferty, K.D., Kuris, A.M., Hechinger, R.F., & Jetz, W. (2008). Homage to Linneaeus: How many parasites? How many hosts? *Proceedings of the National Academy of Sciences, 105* (Supp.1), 11482-11489.

537
538 Dubut, V., Fouquet, A., Voisin, A., Costedoat, C. Chappaz, R. & Gilles, A. (2012). From Late Miocene to Holocene: processes of differentiation within the *Telestes* genus (Actinopterygii: Cyprinidae). *PLoS ONE*, 7(3), e34423.

541
542 Durand, J.D., Persat, H. & Bouvet, Y. (1999). Phylogeography and postglacial dispersion of the chub (*Leuciscus cephalus*) in Europe. *Molecular Ecology*, 8, 989-997.

544
545 Durand, J.D., Ünlü, E., Doadrio, I., Pipoyan, S. & Templeton, A.R. (2000). Origin, radiation, dispersion and allopatric hybridization in the chub, *Leuciscus cephalus*. *Proceedings of the Royal Society of London B*, **267**, 1687-1697.

548
549 Durand, J.D., Bianco, P.G., Laroche, J. & Gilles, A. (2003). Insight into the origin of endemic Mediterranean ichthyofauna: phylogeography of *Chondrostoma* genus (Teleostei, Cyprinidae). *Journal of Heredity*, **94**, 315-328.

552
553 Gante, H.F. (2011). Diversification of circum-Mediterranean barbels. In: Grillo, O., Venora, G. (Eds.), *Changing Diversity in Changing Environment*. InTech, pp. 283-298.

555
Guillen-Hernandez, S. & Whitfield, P.J. (2001). A comparison of freshwater and marine/estuarine strains of *Pomphorhynchus laevis* occurring sympatrically in flounder, *Platichthys flesus*, in the tidal Thames. *Journal of Helminthology*, **75**, 237-243.

Guindon, S., Dufayard, J.F., Lefort, V., Anisimova, M., Hordijk, W. & Gascuel, O. (2010). New algorithms and methods to estimate maximum-likelihood phylogenies: assessing the performance of PhyML 3.0. *Systematic Biology*, **57**, 307-321.

Hall T (1999) BioEdit: a user-friendly biological sequence alignment editor and analysis program for Windows 95/98/NT. *Nucleic Acids Symposium Series*, **41**, 95-98.

Hasegawa, M., Kishino, H. & Yano, T. (1985). Dating the human-ape splitting by a molecular clock of mitochondrial DNA. *Journal of Molecular Evolution*, **22**, 160-174.

Hoberg, E.P. & Brooks, D.R. (2008). A macroevolutionary mosaic: episodic host-switching, geographical colonization and diversification in complex host-parasite systems. *Journal of Biogeography*, **35**, 1533-1550.

Hou, Z.E., Sket, B., Fiser, C. & Li, S. (2011). Eocene habitat shift from saline to freshwater promoted Tethyan amphipod diversification. *Proceedings of the National Academy of Sciences of the United States of America*, **108**, 14533-14538.

Kennedy, C.R. (2006). Ecology of the Acanthocephala. Cambridge University Press, Cambridge, UK.
581 Knowlton N., Weigt, L.A. (1998). New dates and new rates for divergence across the Isthmus of Panama. *Proceedings of the Royal Society of London B, **265***, 2257-2263.

582
583
584 Kotlík, P. & Berrebi, P. (2001). Phylogeography of the barbel (*Barbus barbus*) assessed by mitochondrial DNA variation. *Molecular Ecology, **10***, 2177-2185.

586
587 Kotlík, P. & Berrebi, P. (2002). Genetic subdivision and biogeography of the Danubian
588 rheophilic barb *Barbus petenyi* inferred from phylogenetic analysis of mitochondrial DNA
589 variation. *Molecular Phylogenetics and Evolution, **24***, 10-18.

590
591 Kotlík, P., Bogutskaya, N.G. & Ekmekci, F.G. (2004). Circum Black Sea phylogeography of
592 *Barbus* freshwater fishes: divergence in the Pontic glacial refugium. *Molecular Ecology, **13*,
593 87-95.

594
595 Levy, A., Doadrio, I. & Almada, V.C. (2009). Historical biogeography of European leuciscins
596 (Cyprinidae): Evaluating the Lago Mare dispersal hypothesis. *Journal of Biogeography, **39***,
597 55-65.

598
599 Librado, P. & Rozas, J. (2009). DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. *Bioinformatics, **25***, 1451-1452.

601
602 Mamos, T., Wattier, R., Burzyński, A. & Grabowskiand, M. (2016). The legacy of a vanished
603 sea: a high level of diversification within a European freshwater amphipod species complex
604 driven by 15 My of Paratethys regression. *Molecular Ecology, **25***, 795-810.

605
Neubauer, T.A., Harzhauser, M., Kroh, A., Georgopoulou, E. & Mandic, O. (2015). A gastropod-based biogeographic scheme for the European Neogene freshwater systems. *Earth-Science Reviews*, **143**, 98–116.

Nylander, J.A.A. (2004). MrModeltest v2. Program distributed by the author. Evolutionary Biology Centre, Uppsala University.

Ogilvie, H.A. Bouckaert, R.R. & Drummond, A. (2017). StarBEAST2 brings faster species tree inference and accurate estimates of substitution rates. *Molecular Biology and Evolution*, **msx126**, doi: 10.1093/molbev/msx126.

O'Mahony, E.M., Bradley, D.G., Kennedy, C.R. & Holland, C.V. (2004). Evidence for the hypothesis of strain formation in *Pomphorhynchus laevis* (Acanthocephala): an investigation using mitochondrial DNA sequences. *Parasitology*, **129**, 341–347.

Orszag-Sperber, F. (2006). Changing perspectives in the concept of “Lago-Mare” in Mediterranean Late Miocene evolution. *Sedimentary Geology*, **188–189**, 259–277.

Papadopoulou, A. & Knowles, L.L. (2016). Toward a paradigm shift in comparative phylogeography driven by trait-based hypotheses. *Proceedings of the National Academy of Sciences*, **113**, 8018–8024.

Perea, S., Böhme, M., Zupančič, P., Freyhof, J., Šanda, R., Özuluğ, M., Abdoli, A. & Doadrio, I. (2010). Phylogenetic relationships and biogeographical patterns in Circum-Mediterranean
630 subfamily Leuciscinae (Teleostei, Cyprinidae) inferred from both mitochondrial and nuclear data. *BMC Evolutionary Biology*, **10**, 265.

632
633 Perdices, A., Doadrio, I., Economidis, P.S., Bohlen, J. & Banarescu, P. (2003). Pleistocene effects on the European freshwater fish fauna: double origin of the cobitid genus *Sabanejewia* in the Danube basin (Osteichthyes: Cobitidae). *Molecular Phylogenetics and Evolution*, **26**, 289-299.

637
638 Perrot-Minnot, M.-J. (2004). Larval morphology, genetic divergence, and contrasting levels of host manipulation between forms of *Pomphorhynchus laevis* (Acanthocephala). *International Journal for Parasitology*, **34**, 45-54.

641
642 Raftery, A, Newton, M, Satagopan, J, & Krivitsky, P. (2007). Estimating the integrated likelihood via posterior simulation using the harmonic mean identity. In: Bernardo, J.M., Bayarri, M.J., Berger, J.O. (Eds.), Bayesian statistics., Oxford University Press, New York, pp. 1-45.

646
647 Rambaut, A, Suchard, M.A., Xie, D., & Drummond, A.J. (2014). Tracer version 1.6, Available at http://beast.bio.ed.ac.uk/Tracer.

649
650 Reyjol, Y., Hugueny, B., Pont, D., Bianco, P. G., Beier, U., Caiola, N., Casals, F., Cowx, I., Economou, A., Ferreira, T., Haidvogl, G., Noble, R., de Sostoa, A., Vigneron, T. & Virbickas, T. (2007). Patterns in species richness and endemism of European freshwater fish. *Global Ecology and Biogeography*, **16**, 65-75.

654
Ronquist, F. & Huelsenbeck, J.P. (2003). MrBayes3: Bayesian phylogenetic inference under mixed models. *Bioinformatics*, **19**, 1572-1574.

Rouchy, J.M. & Caruso, A. (2006). The Messinian Salinity Crisis in the Mediterranean basin: a reassessment of the data and an integrated scenario. *Sedimentary Geology*, **188-189**, 35-67.

Šedivá, A., Janko, K., Slechtová, V., Kotlík, P., Simonović, P., Delic, A. & Vassilev, M. (2008). Around or across the Carpathians: colonization model of the Danube basin inferred from genetic diversification of stone loach (*Barbatula barbatula*) populations. *Molecular Ecology*, **17**, 1277-1292.

Seifertová, M., Bryja, J., Vyskočilová, M., Martínková, N. & Šimková, A. (2012). Multiple Pleistocene refugia and postglacial colonization in the European chub (*Squalius cephalus*) revealed by combined use of nuclear and mitochondrial markers. *Journal of Biogeography*, **39**, 1024-1040.

Solà, E., Sluys, R., Gritzalis, K. & Riutort, M. (2013). Fluvial basin history in the northeastern Mediterranean region underlies dispersal and speciation patterns in the genus *Dugesia* (Platyhelminthes, Tricladida, Dugesiidae). *Molecular Phylogenetics and Evolution*, **66**, 877-888.

Špakulová, M., Perrot-Minnot, M.-J., & Neuhaus, B. (2011). Resurrection of *Pomphorhynchus tereticollis* (Rudolphi, 1809) (Acanthocephala: Pomphorhynchidae) based on new morphological and molecular data. *Helminthologia*, **48**, 268-277.
680 Tamura, K., Peterson, D., Peterson, N., Stecher, G., Nei, M., & Kumar, S. (2011). MEGA5: 681 Molecular Evolutionary Genetics Analysis using Maximum Likelihood, Evolutionary Distance, 682 and Maximum Parsimony Methods. *Molecular Biology and Evolution*, **28**, 2731-2739.

683 Thomas, J.A., Welch, J.J., Woolfit, M. & Bromham, L. (2006). There is no universal molecular 684 clock for invertebrates, but rate variation does not scale with body size. *Proceedings of the 685 National Academy of Sciences of the United States of America*, **103**, 7366-7371.

687 Toucanne, S., Zaragosi, S., Bourillet, J.F., Cremer, M., Eynaud, F., Van Vliet-Lanoë, B., 688 Penaud, A., Fontanier, C., Turon, J.L., Cortijo, E. & Gibbard, P.L. (2009). Timing of massive 689 ‘Fleuve Manche’ discharges over the last 350 kyr: insights into the European ice-sheet 690 oscillations and the European drainage network from MIS 10 to 2. *Quaternary Science Reviews*, 691 **28**, 1238-1256

693 Tsigenopoulos, C. S., Durand, J.-D., Ünlü, E. & Berrebi, P. (2003). Rapid radiation of the 694 Mediterranean *Luciobarbus* species (Cyprinidae) after the Messinian Salinity Crisis of the 695 Mediterranean Sea, inferred from mitochondrial phylogenetic analysis. *Biological Journal of 696 the Linnean Society*, **80**, 207-222.

698 van der Mescht, L., Matthee, S. & Matthee, C.A. (2015). Comparative phylogeography between 699 two generalist flea species reveal a complex interaction between parasite life history and host 700 vicariance: parasite-host association matters. *BMC Evolutionary Biology*, **15**, 105. DOI 701 10.1186/s12862-015-0389-y.

703
704 Vardić-Smrzlić, I., Valić, D., Kapetanović, D., Filipović Marijić, V. Gjurčević, E. & Teskeredžić, E. (2015). *Pomphorhynchus laevis* (Acanthocephala) from the Sava River basin: New insights into strain formation, mtDNA-like sequences and dynamics of infection. *Parasitology International*, **64**, 243-250.

708
709 Yang, Z. (1994). Estimating the pattern of nucleotide substitution. *Journal of Molecular Evolution*, **39**, 105-111.

711
712 Yu, Y., Harris, A.J. & He, X.J. (2010). S-DIVA (Statistical Dispersal-Vicariance Analysis): a tool for inferring biogeographic histories. *Molecular Phylogenetics and Evolution*, **56**, 848-850.

714
715 Yu, Y., Harris, A.J. & He, X.J. (2014). RASP (Reconstruct Ancestral State in Phylogenies)
716 3.0. Available at http://mnh.scu.edu.cn/soft/blog/RASP.

717
718 Yu, Y., Harris, A.J., Blair, C. & He, X.J. (2015). RASP (Reconstruct Ancestral State in Phylogenies): A tool for historical biogeography. *Molecular Phylogenetics and Evolution*, **87**, 46-49

721
722 SUPPORTING INFORMATION

723 Additional Supporting Information may be found in the online version of this article:

724
725 Appendix S1 Detailed information on *P. laevis* and *P. tereticollis* taxonomy and host range, on host phylogeography; and additional sampling information and methods

727 Appendix S2: Demographic analysis, ITS analysis (phylogenetic tree) and COI maximum-likelihood tree, and biogeographic analysis
728
729 Appendix S3: Supplementary discussion on an alternative scenario based on calibration point
730 at Pl_L5, on possible human-mediated introduction, and on the putative significance of
731 pseudogenes

732
733 **Data Accessibility:**
734 All new COI and ITS sequences were deposited in the EMBL database 'European Nucleotide
735 Archive' under the accession numbers LN994840 to LN995000, and LN995001 to LN995058,
736 respectively, and are available at http://www.ebi.ac.uk/ena/data/view/LN994840-LN995058.

737
738 **Biosketch**
739 Marie-Jeanne Perrot-Minnot is an evolutionary ecologist, whose main interest is the evolution
740 of parasites with complex life cycle. The authors gather their expertise in parasite systematics
741 (M.S.), phylogeography (P. K., C. T., R. W.) and parasite ecology (A. A., S. D.), to better
742 understand the historical biogeography of these two ecologically important parasite species.
743 Author contributions: A.A., S.D., P.K., M-J.P.-M., and M.S. provided samples; M-J.P.-M. and
744 R.W. gathered the dataset; C.T. performed the analysis, M-J.P.-M., C.T., P.K. and R. W.
745 contributed to data interpretation; M-J.P.-M. and C.T. wrote a first draft; and all authors
746 contributed to writing.

747
748 Editor: Michelle Gaither

749
Table 1: Genetic distance based on COI sequences within and between (a) *Pomphorhynchus laevis* and (b) *P. tereticollis* lineages. The species used as outgroup for *P. laevis* and *P. tereticollis* were *P. tereticollis* and *P. laevis*, respectively. Genetic distance between lineages is given below the diagonal, and standard error above the diagonal.

753

<table><thead><tr><th rowspan="2">(a)</th><th rowspan="2">Lineage</th><th colspan="2">Distance within</th><th colspan="6">Distance between lineages</th></tr><tr><th>lineage (SE)</th><th>Pl_L1</th><th>Pl_L2</th><th>Pl_L3</th><th>Pl_L4</th><th>Pl_L5</th><th>Outgroup</th></tr></thead><tbody><tr><td rowspan="5">P. laevis</td><td>Pl_L1</td><td>0.024 (0.004)</td><td></td><td>0.011</td><td>0.012</td><td>0.015</td><td>0.019</td><td>0.022</td></tr><tr><td>Pl_L2</td><td>0.023 (0.004)</td><td>0.107</td><td></td><td>0.014</td><td>0.015</td><td>0.021</td><td>0.023</td></tr><tr><td>Pl_L3</td><td>0.017 (0.004)</td><td>0.105</td><td>0.117</td><td></td><td>0.014</td><td>0.019</td><td>0.020</td></tr><tr><td>Pl_L4</td><td>0.012 (0.003)</td><td>0.118</td><td>0.148</td><td>0.125</td><td></td><td>0.016</td><td>0.020</td></tr><tr><td>Pl_L5</td><td>0.015 (0.003)</td><td>0.193</td><td>0.203</td><td>0.186</td><td>0.181</td><td></td><td>0.022</td></tr><tr><td>Outgroup</td><td></td><td>0.017 (0.003)</td><td>0.236</td><td>0.266</td><td>0.217</td><td>0.231</td><td>0.230</td><td></td></tr></tbody></table>

<table><thead><tr><th>(b)</th><th>Distance within<br>lineage (SE)</th><th>Distance between Pt<br>and outgroup</th></tr></thead><tbody><tr><td><i>P. tereticollis</i></td><td>0.016 (0.003)</td><td>0.018</td></tr><tr><td>Outgroup</td><td>0.142 (0.011)</td><td>0.244</td></tr></tbody></table>

754

755
<table><caption>Table 2: Genetic diversity indices from COI sequences for each lineage of *Pomphorhynchus laevis* and for *P. tereticollis*.</caption><thead><tr><th>Species</th><th>Lineage</th><th>n<sup>a</sup></th><th>nh<sup>b</sup></th><th>π (sd)<sup>c</sup></th><th>h<sup>d</sup></th><th>k<sup>e</sup></th></tr></thead><tbody><tr><td rowspan="5">P. laevis</td><td>Pl_L1</td><td>50</td><td>22</td><td>0.021 (0.001)</td><td>0.922</td><td>13.049</td></tr><tr><td>Pl_L2</td><td>39</td><td>19</td><td>0.014 (0.003)</td><td>0.830</td><td>8.669</td></tr><tr><td>Pl_L3</td><td>8</td><td>5</td><td>0.011 (0.002)</td><td>0.857</td><td>6.964</td></tr><tr><td>Pl_L4</td><td>9</td><td>7</td><td>0.013 (0.001)</td><td>0.917</td><td>7.778</td></tr><tr><td>Pl_L5</td><td>5</td><td>5</td><td>0.017 (0.004)</td><td>1.000</td><td>10.500</td></tr><tr><td>P. tereticollis</td><td></td><td>50</td><td>28</td><td>0.015 (0.001)</td><td>0.968</td><td>9.220</td></tr></tbody></table>
Table 3: Biogeographic history inferred from S-DIVA analysis based on COI sequences for (a) *Pomphorhynchus laevis* lineages (except for peri-Mediterranean ones), (b) *P. tereticollis* lineages. The reconstructed history includes the divergence time estimates (from the analysis with StarBEAST2) and geographic distribution of the most recent common ancestor (MRCA), and the biogeographic scenario leading to derived populations (disp. dispersal, vic. vicariance, ext. extinction) and its associated probability (> 17%). Abbreviations for biogeographic units are detailed in Table S1.1. (a): B, Pannonian area of Ponto-Caspian Europe; C, Western Europe; E, Carpathian area of Ponto-Caspian Europe; F, Balkan area of Ponto-Caspian Europe (Danubian southern tributaries); G, Eastern peri-Mediterranean (North Aegean-Anatolian); H, river Danube in Ponto-Caspian Europe; J, Caspian part of Ponto-Caspian Europe; K, Eastern peri-Mediterranean (North Aegean and West Anatolian). (b): A, British Isles; B, Baltic Sea; C, Western Europe; D, North Sea; E, Carpathian area. The symbol ‘x’ combines two or more biogeographic units covered by a derived population; na: not available

776
<table><caption>(a)</caption><thead><tr><th>P. laevis ancestral population (MRCA to derived populations)</th><th>Time to MRCA</th><th>Ancestral biogeographic distribution (% of alternative biogeographic area)</th><th>Biogeographic events</th><th>Derived populations or lineages</th><th>Ref. node (1)</th></tr></thead><tbody><tr><td>Ponto-Caspian x Western Europe</td><td>2.55 Myr</td><td>EH - BH - CH (33:33:33)</td><td>EH -> ECH -> EC | H (P = 0.17)</td><td>disp. (1) vic. (1)</td><td>Carpathian x Western Europe // Danubian (Pl_L1, Pl_L2)</td><td>103</td></tr><tr><td>Danubian x Caspian</td><td>1.27 Myr</td><td>HJ (100)</td><td>HJ -> H | J (P=1)</td><td>vic. (1)</td><td>Danube, Vistula // Volga (Pl_L2a, Pl_L2b)</td><td>79-81</td></tr><tr><td>Carpathian x Western Europe</td><td>0.22 Myr</td><td>EC - EB (50:50)</td><td>EC -> BEC -> E | BC (P = 0.25)</td><td>disp. (1) vic. (1)</td><td>Carpathian // Pannonian x Western Europe</td><td>102</td></tr><tr><td>Carpathian</td><td>0.50 Myr</td><td>E (100)</td><td>E -> E ^E -> E ^EG -> E | EG (P = 1)</td><td>disp. (1)</td><td>Eastern Carpathian // Eastern Carpathian x North Aegean (Pl_L1d, Pl_L1b + c)</td><td>91</td></tr><tr><td>Eastern Carpathian x North Aegean</td><td>0.50 Myr</td><td>EG (100)</td><td>EG -> EGF -> EF | G (P = 1)</td><td>disp. (1) vic. (1)</td><td>Eastern Carpathian x Balkans // North Aegean (Pl_L1b, Pl_L1c)</td><td>87</td></tr><tr><td>Pannonian x Western Europe</td><td>na</td><td>BC - B (50:50)</td><td>BC -> B ^BC -> B | BC (P=0.25)</td><td>disp. (1)</td><td>Pannonian // Pannonian x Western Europe (within Pl_L1a)</td><td>101</td></tr><tr><td>Eastern Carpathian x Balkans</td><td>na</td><td>EF (100)</td><td>EF -> E | F (P=1)</td><td>vic. (1)</td><td>Eastern Carpathian // Balkans (within Pl_L1b)</td><td>85</td></tr></tbody></table>

<footnote>(1) Ref. node in Fig. S2.3</footnote>
(b)

<table><thead><tr><th><i>P. tereticollis</i><br>ancestral population<br>(MRCA to derived<br>populations)</th><th>Time to<br>MRCA</th><th>Ancestral<br>biogeographic<br>distribution<br>(% of alternative<br>biogeographic area)</th><th>Biogeographic events</th><th>Derived populations or lineages</th><th>Ref.<br>node<br>(2)</th></tr></thead><tbody><tr><td>ancestral to all<br>lineages</td><td>0.42 Myr</td><td>C (100)</td><td>C-&gt;C^C-&gt;CE^C-&gt;C|CE<br>(P=0.6)</td><td>disp. (1)</td><td>Western Europe (Rhone drainage) //<br>Western Europe x Carpathian</td><td>59</td></tr><tr><td>Western Europe x<br>Carpathian</td><td>0.42 Myr</td><td>CE or C (61:39)</td><td>CE-&gt;CE^C-&gt;CEB^C-&gt;<br>CB|CE (P=0.34)</td><td>disp. (2)</td><td>North-western Europe (Pt_L1) //<br>Western and Central Europe (Pt_L2)</td><td>56</td></tr><tr><td>North-western<br>Europe</td><td>0.29 Myr</td><td>CB (100)</td><td>CB-&gt;C|B (P=1) and<br>B-&gt;B^B-&gt;AB^B-&gt;AB|B<br>(P=1)</td><td>disp. (1)<br>vic. (1)</td><td>Western Europe (Rhone drainage) //<br>North-western Europe (Baltic sea,<br>British islands)</td><td>44, 41</td></tr><tr><td>Western and Central<br>Europe</td><td>na</td><td>CE or E (52:48)</td><td>CE-&gt;E-&gt;E^E-&gt;EA^E-&gt;E|AE<br>(P=0.23)</td><td>disp. (2)<br>vic. (1)<br>ext. (1)</td><td>Eastern Carpathian // Southern<br>Carpathian - British Islands</td><td>55</td></tr><tr><td>Central Carpathian -<br>British Islands</td><td>na</td><td>AE or CE or DE<br>(40:35:25)</td><td>AE-&gt;EAC-&gt;E|AC (P=0.17)</td><td>disp. (1)<br>vic. (1)</td><td>Eastern Carpathian // Western and<br>North Europe, British Islands<br>(within Pt_L2)</td><td>50</td></tr></tbody></table>

778

779

<footnote>(2) Ref. node in Fig. S2.4</footnote>
<figure><img src="image_2.png" /></figure>

781
782 Figure 1: Geographical location of sampling sites for the fish acanthocephalan parasites, *Pomphorhynchus laevis* and *P. tereticollis*, within the Western Palaearctic. Outlined shapes are major biogeographic regions identified from the composition of European riverine fish fauna by Reyjol *et al*. (2007), and colours refer to sub-units within each region. For each sampling site, the number refers to the locality, and the symbol to the affiliation of haplotypes to one of the two *Pomphorhynchus* species (circle, *P. laevis*; triangle, *P. tereticollis*) or both (square) based on sequencing data. Full details on the samples (biogeographical sub-unit, locality, host species, haplotype and accession numbers) are given in Table S1.1.
783
784
785
786
787
788
789
790
791
792
793
<figure><img src="image_3.png" /></figure>
<figure><img src="image_4.png" /><figcaption>Figure 2: Bayesian tree reconstructed from COI sequences of (a) *Pomphorhynchus laevis*, (b) *P. tereticollis*, and the species used as outgroup (*P. tereticollis* and *P. laevis*, respectively). Branch support values are maximum-likelihood bootstrap (≥ 50%) percentages and Bayesian posterior probabilities (≥ 0.70). Haplotypes are indicated by numbers as given in Table S1.1. Colour chart is corresponding to biogeographic sub-units as defined in Fig. 1 and Table S1.1.</figcaption></figure>

799

800

801

802

803

804

805
<figure><img src="image_5.png" /><figcaption>Figure 3: Species tree chronogram for *Pomphorhynchus laevis* and *P. tereticollis* estimated by STARBEAST2 based on COI and ITS data. Numbers at nodes are for divergence times (Myr) and their 95% confidence intervals estimated from a geologic calibration point, 6.05 ± 0.09 Myr (Messinian Salinity Crisis, Rouchy & Caruso, 2006). Black circles are for nodes supported with PP < 0.70 meaning that molecular dating at these nodes should be taken with caution.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 4: Median-joining network of the COI haplotypes from *Pomphorhynchus laevis* and *P. tereticollis*. Line lengths in the network reflect the number of mutational changes (in brackets when >1) between haplotypes, and the size of the circles is proportional to the frequencies of the represented haplotype. Black dots represent hypothetical missing or unsampled ancestral haplotypes. Haplotypes are indicated by numbers as given in Table S1.1. Colour coding corresponds to geographic areas as defined in Fig. 1 and Table S1.1.</figcaption></figure>