# It's a small world for parasites: evidence supporting the North American invasion of European Echinococcus multilocularis

Maria A. Santa, Gérald Umhang, Claudia Klein, Danielle M. Grant, Kathreen E. Ruckstuhl, Marco Musiani, John S. Gilleard, Alessandro Massolo

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

► To cite this version:

Maria A. Santa, Gérald Umhang, Claudia Klein, Danielle M. Grant, Kathreen E. Ruckstuhl, et al. It's a small world for parasites: evidence supporting the North American invasion of European Echinococcus multilocularis. Proceedings of the Royal Society B: Biological Sciences, 2023, 290 (1994), pp.20230128. 10.1098/rspb.2023.0128. anses-04127178

HAL Id: anses-04127178

https://anses.hal.science/anses-04127178

Submitted on 13 Jun 2023

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>
PROCEEDINGS B

royalsocietypublishing.org/journal/rspb

Research

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

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

Check for
updates

**Cite this article:** Santa MA, Umhang G, Klein C, Grant DM, Ruckstuhl KE, Musiani M, Gilleard JS, Massolo A. 2023 It's a small world for parasites: evidence supporting the North American invasion of European *Echinococcus multilocularis*. *Proc. R. Soc. B* **290**: 20230128. https://doi.org/10.1098/rspb.2023.0128

Received: 17 January 2023

Accepted: 13 February 2023

**Subject Category:**

Ecology

**Subject Areas:**

ecology, health and disease and epidemiology,
genetics

**Keywords:**

parasite invasion, *Echinococcus multilocularis*,
North America, genetic characterization,
allochthonous strains

**Author for correspondence:**

Alessandro Massolo
e-mail: alessandro.massolo@unipi.it

Electronic supplementary material is available
online at https://doi.org/10.6084/m9 figshare.
c.6440232.

# It's a small world for parasites: evidence supporting the North American invasion of European *Echinococcus multilocularis*

Maria A. Santa<fnref n="1" />, Gerald Umhang<fnref n="2" />, Claudia Klein<fnref n="3" />,<fnref n="4" />, Danielle M. Grant<fnref n="3" />,<fnref n="5" />, Kathreen E. Ruckstuhl<fnref n="1" />, Marco Musiani<fnref n="1" />,<fnref n="6" />, John S. Gilleard<fnref n="7" /> and Alessandro Massolo<fnref n="8" />,<fnref n="9" />,<fnref n="10" />

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

MAS, 0000-0003-2512-1573; AM, 0000-0002-6333-4281

*Echinococcus multilocularis* (*Em*), the causative agent of human alveolar echino-
coccosis (AE), is present in the Holarctic region, and several genetic variants
deem to have differential infectivity and pathogenicity. An unprecedented out-
break of human AE cases in Western Canada infected with a European-like
strain circulating in wild hosts warranted assessment of whether this strain
was derived from a recent invasion or was endemic but undetected. Using
nuclear and mitochondrial markers, we investigated the genetic diversity of
*Em* in wild coyotes and red foxes from Western Canada, compared the genetic
variants identified to global isolates and assessed their spatial distribution to
infer possible invasion dynamics. Genetic variants from Western Canada
were closely related to the original European clade, with lesser genetic diversity
than that expected for a long-established strain and spatial genetic discontinu-
ities within the study area, supporting the hypothesis of a relatively recent
invasion with various founder events.

1. Introduction

In the current Anthropocene epoch, biological invasions, due to trans- and inter-
continental movements, are generating a global biotic homogenization, influen-
cing at the same time global patterns of disease [1,2]. Such invasions can lead
to the introduction of novel hosts and their parasites, opening up opportunities
for the emergence of diseases that can lead to outbreaks or the establishment of
new endemism [3]. Indeed, many aetiological agents of emerging infectious dis-
eases can be considered biological invaders [4]. As an example, the cestode
*Echinococcus multilocularis* (*Em*), the causative agent of alveolar echinococcosis
(AE), is considered an emerging or re-emerging disease in some regions, becom-
ing an important public health concern in many countries worldwide [5,6]. This
parasite is distributed across the northern hemisphere and is transmitted in sylva-
tic and semi-synanthropic cycles, involving small mammals as intermediate hosts

© 2023 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution
License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original
author and source are credited.

<footnote n="1">¹Department of Biological Sciences, Faculty of Science, University of Calgary, Calgary, Alberta, Canada</footnote>

<footnote n="2">²Nancy Laboratory for Rabies and Wildlife, National Reference Laboratory for *Echinococcus* spp., Wildlife Surveillance and Eco-epidemiology Unit, ANSES, Technopôle Agricole et Vétérinaire, Malzéville, France</footnote>

<footnote n="3">³Department of Clinical and Veterinary Clinical Sciences, Faculty of Veterinary Medicine, University of Calgary, Calgary, Alberta, Canada</footnote>

<footnote n="4">⁴Friedrich-Loeffler-Institut, Bundesforschungsinstitut für Tiergesundheit, Federal Research Institute for Animal Health, Neustadt, Germany</footnote>

<footnote n="5">⁵NORCE Climate, NORCE Norwegian Research Centre, Bjerknes Centre for Climate Research, Bergen, Norway</footnote>

<footnote n="6">⁶Dipartimento di Scienze Biologiche, Geologiche e Ambientali (BiGeA), University of Bologna, Italy</footnote>

<footnote n="7">⁷Department of Comparative Biology & Experimental Medicine, Host-Parasite Interactions (HPI) program, Faculty of Veterinary Medicine, University of Calgary, Calgary, Alberta, Canada</footnote>

<footnote n="8">⁸Ethology Unit, Department of Biology, University of Pisa, Pisa, Italy</footnote>

<footnote n="9">⁹UMR CNRS 6249 Chrono-environnement, Université Bourgogne Franche-Comté, Besançon, France</footnote>

<footnote n="10">¹⁰Faculty of Veterinary Medicine, University of Calgary, Calgary, Alberta, Canada</footnote>
(primarily arvicoline rodents) and canids as definitive hosts (e.g. red foxes (*Vulpes vulpes*), coyotes (*Canis latrans*) and dogs). In small mammals, the parasite invades and develops in internal organs (primarily the liver), being highly pathogenic. Therefore, the host either dies from the infection or is easily preyed on by the definitive hosts [6]. Conversely, wild canids only develop an intestinal infection, causing no substantial pathology [7]. Humans can act as dead-end hosts and develop AE by accidental ingestion of eggs via contaminated food, water or through contact with faeces of infected definitive hosts [8].

Although the geographical distribution of *E. multilocularis* is limited to the Northern hemisphere, its distribution and prevalence in wild and domestic hosts are increasing across its range, likely due to human activities [9–11]. The parasite has been reported in new areas previously considered non-endemic [12–15]. Furthermore, synanthropic hosts and domestic animals now have important roles in its transmission, contributing to the establishment of urban cycles [16,17]. Additionally, increased translocation of domestic and wild animals has generated new host species assemblages and introduced *E. multilocularis* (or particular genetic variants) into new areas [9,18,19].

Using concatenated sequences of three mitochondrial genes (*cox1*, *nad2* and *cob*), *Em* genetic variants have been grouped into four clades, mostly linked with their geographical origins: Asian, European, North American and one restricted to Mongolia [20]. However, Asian haplotypes have also been reported in Western Russia [21], Saint Lawrence Island (Alaska) [20] and, more recently, in central-eastern Europe (Poland) [22,23]. In addition, haplotypes clustering most closely with the North American clade have also been detected in Siberia [21]. In North America, the situation is similar, with genetic variants linked to geographically distant regions. In recent years, European-like haplotypes have been detected in Western Canada (Alberta, British Columbia, and Saskatchewan) in wild and domestic hosts [12,18,24–26]. Additionally, despite a historically low incidence of human AE in continental North America, with only two locally acquired cases ever reported (in 1923 and 1977, respectively), at least 17 cases have been described in the province of Alberta since 2013 [25,27], along with a case recently confirmed in the adjacent province of Saskatchewan [28], and more recently, two cases in the US [29,30]. Molecular characterization of six of the Albertan human cases implicated a European-like haplotype (ECA) that was also detected in local wild hosts [25,26]. Similarly, European-like haplotypes (different from those found in Canada) were identified in the two human cases in the US [29,30]. However, the distribution of European-like strains, their origin in North America, and potential interactions (e.g. competition or hybridization) with native North American genetic variants have not been elucidated. It has been suggested that the European strain could be more pathogenic and/or virulent than the North American strain, based on the historically low incidence of human cases of AE reported for the region. Consequently, it is of great importance to genetically characterize *E. multilocularis* in Western Canada, to assess the extent of the spread of European strains, their distribution in the main definitive hosts (coyotes and red foxes) and to elucidate potential sources of invasion of these strains.

In Canada, genetic characterization of *E. multilocularis* has been conducted primarily using mitochondrial (mtDNA) markers, yet these markers have relatively low variability. Conversely, the microsatellite EmsB, a multi-locus nuclear DNA marker, has greater discriminatory power than classic

mtDNA markers, and has been used to identify spatio-temporal characteristics of *E. multilocularis* transmission in several European countries, at continental, national and local geographical scales [31–34]. Additionally, the combination of mitochondrial and nuclear markers has greater discriminatory power to identify genetic profiles and detect potential introgression events among genetic variants [22].

With this study, we aimed to: (i) assess the genetic diversity of *E. multilocularis* in Western Canada, based on EmsB profiles and mtDNA haplotypes; (ii) investigate genetic relationships with *E. multilocularis* isolates from other geographical regions across the globe, to identify possible sources of invasion; and (iii) evaluate the spatial distribution of the detected genetic variants to elucidate the spread dynamics of the European-like strains, possible spatio-temporal scenarios of its invasion process, and the likely origins of these European strains in Western Canada.

## 2. Materials and methods

### (a) Parasite collection and DNA extraction

*Echinococcus multilocularis* specimens were collected from gastrointestinal (GI) tracts of red foxes and coyotes of either road-killed or trap-harvested animals (trapped for purposes independent of this study), collected between 2012 and 2017 in Western Canada. Trapped animals were obtained from licensed trappers with the collaboration of the Alberta Trappers Association. GI tracts were screened using a modification of the scraping, filtration and counting technique, to identify and collect *Echinococcus* spp. worms [35,36]. We analysed *Em* worms from 70 coyotes and 13 foxes from northern, central and southern Alberta (AB); four coyotes from north-west British Columbia (BC); and 10 coyotes from southeast Saskatchewan (SK). Extraction of DNA was performed on up to five individual worms per host using the Nucleospin 96 Tissue Kit (Macherey-Nagel, Germany) for samples processed in France (Anses Nancy Laboratory for Rabies and Wildlife) and the E.Z.N.A. MicroElute Genomic DNA Kit (Omega Bio-tek, US) for samples processed in Canada (University of Calgary, Faculty of Veterinary Medicine). Extraction was performed following the manufacturer's instructions, and DNA was stored at -20°C until processed.

### (b) Genetic characterization using mtDNA

Genetic diversity of *E. multilocularis* was characterized by sequencing the genes *nad2* (1068 bp), *cob* (882 bp) and *cox1* (1608 bp). We used one to five worms per host, depending on worm availability and quality (intactness), using primers and PCR parameters as previously described [20]. Sequences obtained from the same *Em* worms from our previous studies were also included [25,26]. The sequences obtained per genetic locus were concatenated and aligned in Geneious 10.0.9 (Biomatters Ltd, New Zealand) and compared to nucleotide databases using the NCBI Nucleotide BLAST tool to identify the strain/haplotype (https://blast.ncbi.nlm.nih.gov). To analyse genetic relationships and possible origins of haplotypes identified in Western Canada, we constructed a haplotype network, based on Hamming distance, including reported European, Asian, Mongolian and North American haplotypes (based on full sequences of *cob*, *nad2* and *cox1*) [20,24,25,29]. R statistical software (R Development Core Team, 2022) and the package *pegas 0.14* [37] were used to construct the haplotype network. A phylogenetic tree was built via partitioned Bayesian analysis using MrBayes 3.2.7 [38], considering three partitions, corresponding to *cob*, *nad2* and *cox1*, and using *Echinococcus granulosus sensu stricto* (G1) as the outgroup. Reversible-jump Markov chain Monte Carlo analyses were run for 1 million generations,
producing 10 000 trees, with the first 2,500 regarded as burn-in.
The tree was plotted using the R-package *phangorn* 2.3.1 [39].

### (c) Characterization of EmsB profiles

A fluorescent PCR assay was used to amplify the EmsB microsatellite, following a previous protocol [34]. Fragment analysis of PCR products was performed by capillary electrophoresis on an automatic sequencer (ABI Prism 310 and ABI 3500/3730; Life Technologies, CA). The resulting EmsB electropherograms were composed of several peaks between 209 and 241 bp. The size (base-pair length) and height (fluorescent signal intensity) of peaks present in each EmsB electropherogram were determined using GeneMapper 5.0 (Life Technologies, CA). Characterization of each EmsB profile was performed as described in the EmsB guidelines from the EWET database (https://ewet-db.univ-fcomte.fr/) [40]. Briefly, peaks below 10% of the sample's maximum peak height were considered artefacts and discarded. To normalize raw data, the height of each peak was divided by the sum of the height of all peaks of a given sample.

### (i) Emsb genotyping, clustering and ordination analysis

Genotyping and hierarchical cluster analysis were performed by calculating the Euclidean distance between profiles and using the average link clustering method UPGMA [41]. The uncertainty of clusters was tested using multi-scale bootstrap resampling (1000 bootstrap replicates), obtaining approximately unbiased *p*-values (au) [42]. A dendrogram was built using *E. granulosus sensu stricto* (G1) as an outgroup. A distance threshold of 0.08 (average genetic distance observed after three generations in the rodent *Meriones unguiculatus*) was used to identify unique EmsB profiles; any branching below this threshold was considered genetically identical [43]. Unique EmsB profiles represented by only one sample were excluded from the analysis as they could not be technically validated [44]. To understand the structure of the clusters detected and unravel relational patterns among the genetic profiles, we performed a non-metric multi-dimensional scaling (NMDS) analysis based on the Euclidean distance matrix, using two and three dimensions to visualize the genetic distance between samples and their arrangement in reduced dimensions [45].

To evaluate genetic relationships with genotyped *E. multilocularis* isolates from Europe (historical and peripheral endemic areas), Asia (Japan, China and Kyrgyzstan) and North America (Canada, Alaska), each unique EmsB profile identified was compared to a reference world collection of profiles from the EWET database (https://ewet-db.univ-fcomte.fr; updated until 2017) [40] and from profiles reported after 2017 [33,44]. Profiles of 1275 samples from 17 countries were first compared to profiles obtained from our study; thereafter, the five profiles with the lowest genetic distance to each of our profiles (P1-P16) were used to build a dendrogram. Hierarchical cluster analysis was performed using R package *pvclust* [46] and NMDS was done using Primer & Permanova + add-on, Version 6 (PRIMER-E Ltd.).

### (ii) Diversity analysis of EmsB profiles from Alberta

We evaluated the diversity of EmsB profiles in the province of Alberta (AB), the most extensively sampled area in our study. We divided this area into six geographical subregions: North-West (AB-NW), North-East (AB-NE), Central-West (AB-CW), Central-East (AB-CE), South-West (AB-SW) and South-East (AB-SE).

We used an integrated approach based on the framework of Hill numbers to assess: (i) sample completeness, (ii) asymptotic diversity estimates to infer true diversities, (iii) non-asymptotic standardization via rarefaction and extrapolation and (iv) evenness of each subregion (except for AB-NW due to low sample size), following the methodology described by Chao *et al.* [47]. First, to evaluate whether sample effort was sufficient and to assess the extent of undetected diversity, sample completeness was estimated for each

subregion based on orders $0 \ge q \le 2$, where *q* is a number that determines the sensitivity of the measures to 'species' (i.e. EmsB profiles) abundances. Second, asymptotic diversity estimates were calculated using a sample size-based rarefaction and extrapolation sampling method [48]. The sample was extrapolated to double the size of the observed sample. If the curve stabilized and levelled off, then the asymptotic estimates were used to infer true diversity. Here, Hill numbers for order $q \ge 0$ included the three most widely used diversity measures (i.e. Chao1 richness estimator, exponential Shannon diversity and inverse Simpson diversity, as special cases of orders $q=0, 1$ and $2$, respectively). Third, if data did not contain sufficient information to accurately infer the true diversity, this was inferred for a standardized sample coverage [49]. Here, we calculated non-asymptotic coverage-based rarefaction and extrapolation estimates for diversity orders $q=0, 1$ and $2$. The standardized sample coverage (i.e. equal fraction of an assemblage's individuals; $C_{max}$) was selected as the minimum among the coverage values for samples extrapolated to double the size of the reference sample. Lastly, the evenness (i.e. distribution of EmsB profiles abundances) was calculated for each subregion. A bootstrap method ($n=100$) was applied to obtain the associated 95% confidence intervals for all estimates.

Analysis of *β*-diversity was performed based on the Bray-Curtis dissimilarity index and complete linkage agglomerative clustering for hierarchical agglomeration [41]. Diversity of EmsB profiles in Alberta was compared to reported diversity in historical European endemic areas (i.e. Switzerland, Germany, Czech Republic, Austria and France) [31,50]. Diversity estimates were recalculated for these areas based on the standardized sample coverage ($C_{max}$). Analyses were performed using the R package *iNEXT-4steps* [47] and *vegan* [51].

### (iii) Analysis of the spatial distribution of EmsB profiles

Alberta is divided into 174 wildlife management units (WMUs). Samples were collected from 38 WMUs and the location of each host was assigned as a random point within a radius of 15 km (minimum radius of a WMU) from the centroid point of the WMU where it was collected. This methodology was used because the exact georeference (latitude, longitude) for all the hosts was not available. A Mantel test was used to assess the hypothesis of genetic isolation by geographical distance, comparing a matrix of the genetic distance between EmsB profiles, and a matrix of the geographical distances between samples based on Euclidean pairwise distances [52]. A Mantel correlogram was used to investigate the underlying structure of the relationship and to measure the correlation between each class of distances [41]. A distance-based redundancy analysis (dbRDA), based on Bray-Curtis dissimilarity matrix of EmsB profiles, was performed using the two spatial variables (latitude and longitude) to estimate the proportion of the genetic variation explained by spatial structures [53]. Statistical significance of the model and variables was evaluated with a permutation test. Analyses were done using the R package *vegan* [51] and Primer & Permanova+add-on, Version 6 (PRIMER-E Ltd.).

## 3. Results

### (a) Genetic relationships of mtDNA haplotypes from North America supported the hypothesis of multiple invasion events of European origin

To characterize the genetic diversity of *E. multilocularis* based on mtDNA, we analysed 96 concatenated sequences of *cox1*, *nad2* and *cob* genes in adult worms recovered from 13 foxes and 77 coyotes. Virtually all (93/96) were identified as the previously described European-like haplotypes ECA, EAB, ESK, ESK2 and BC1 [25,26], and only three as the North American
<figure><img src="image_6.png" /><figcaption>**Figure 1.** Genetic relationships among haplotypes of *Echinococcus multilocularis* from Western Canada and other historically endemic regions. (a) Mitochondrial haplotype network based on Hamming distance, including European, Asian and North American haplotypes previously reported (white circles) [20,24,25], compared to haplotypes from Western Canada detected in our study (grey circles). The network was constructed based on concatenated sequences of mitochondrial genes *nad2*, *cob* and *cox1*. The ECA haplotype (in red) is the most prevalent in Alberta and has been associated with the most recent cases of human AE in this province. (b) Phylogenetic tree inferred by partitioned Bayesian analysis performed on concatenated mitochondrial DNA, using *E. granulosus* (G1) as an outgroup. Values on tree nodes are Bayesian posterior probabilities. The 0.002 scale bar denotes genetic distance (nucleotide substitutions per site).</figcaption></figure>

haplotype N2 (in two coyotes and one fox). Furthermore, the per cent identity to previously described haplotypes was 100%, with no new SNPs identified. The haplotypes ECA, EAB and N2 were present only in Alberta; ESK2 was in Saskatchewan and Alberta; ESK only in Saskatchewan; and BC1 only in British Columbia. The ECA haplotype was the most prevalent, with a frequency of 78.1% (75/96), followed by ESK2 with a fre quency of 11.5% (11/96). Haplotypes BC1, EAB and N2 had a frequency of 3.1% (3/96) each, and only one sample was ident ified as ESK. The haplotype network showed that the European- like haplotypes found in Western Canada were closely related to the European clade, with few mutational steps (range, 1–6) between them, and that these were more similar to the European haplotype E4 (AB461395.1, AB461404.1, AB461414.1) than any other, suggesting a recent origin derived from the European clade. The most genetically distant haplotype was one from British Columbia (BC1), with six mutational steps to E4, fol lowed by one from South-East Alberta and Saskatchewan (ESK2). Conversely, the haplotypes ECA, EAB and ESK were closely related, differing by only one or two mutations, having the lowest number of mutational steps to E4 (figure 1a), indicat ing a possible stem from a single invasion event. The phylogenetic tree confirmed the close relationship between the European-like haplotypes in Western Canada and the original European clade. Within this clade, the three similar haplotypes (EAB, ECA and ESK) were in the same branch with the E4, whereas BC1 was in another (along with the haplotype SK1), and more deeply branched, indicating a distinct most recent common ancestor. These branching nodes were supported by the Bayesian posterior probabilities (figure 1b).

(b) Cluster and ordination analysis of EmsB profiles
from Western Canada confirms the invasion
hypothesis

Through fragment analysis of the multi-locus EmsB micro-
satellite, 204 individual EmsB profiles were computed and
technically validated [40] (doi:10.6084/m9.fishshare.16818607).

Comparing the similarity of these 204 profiles, we identified 16 unique EmsB profiles based on the 0.08 threshold average genetic distance. These profiles formed three significant clusters: the first with samples from AB only, the second with samples from AB and SK, and the third with samples from AB and BC (figure 2a).

Two profiles (P13 and P15; exclusively from AB) were the most frequent, each representing 32.8% of the worms. When comparing the two definitive hosts, profile P13 was present in 69% of foxes and 64% of coyotes, and was the most prevalent in both hosts, representing 56% and 45% of worms in foxes and coyotes, respectively. Profile P4, the third most common (11.8%), was present in AB and SK. Profile P11 was only in BC and AB, corresponding to 5.9% of the worms. The other profiles represented 1–3.4% of worms and were primarily from AB (electronic supplementary material, table S1).

When comparing our profiles with profiles from the EWET world database (1275 samples, from 17 countries) to evaluate genetic relationships with other *E. multilocularis* isolates, 13 of our 16 unique profiles (81.2%) clustered with isolates from Europe, with eight being genetically indistinguishable (based on the 0.08 threshold) from profiles present in several European countries (i.e. Germany, France, Poland, Austria, Switzerland, Denmark, Slovakia and Czech Republic), and only three profiles clustering with samples from the Arctic (i.e. Saint Lawrence Island (Alaska) and Svalbard archipelago (Norway); figure 2b; electronic supplementary material, table S2).

When the NMDS ordination analysis was performed, the representation of the genetic distances between profiles in two dimensions had the lowest stress (0.05; indicating the best representation), with the three clusters identified being linked with the location (i.e. provinces) (figure 3a) and the mtDNA haplotypes (figure 3b).

(c) Genetic diversity of EmsB profiles in Alberta was
lower than would be expected in endemic areas

The sample completeness profiles increased with diversity order ($q \ge 0$), implying that there was undetected diversity
<figure><img src="image_7.png" /><figcaption>**Figure 2.** Cluster analysis of EmsB profiles of *Echinococcus multilocularis* from Western Canada and correspondence with the most similar profiles globally. (a) Dendrogram constructed by hierarchical clustering analysis based on EmsB genotypic data from *E. multilocularis* worms collected in Alberta (AB), Saskatchewan (SK) and British Columbia (BC) from 2012 to 2017. Two samples of *E. granulosus sensu stricto* (G1) were used as an outgroup. Approximately unbiased *p*-values (red numbers on nodes, in per cent) were calculated with multi-scale bootstrap resampling (1000 bootstrap replicates). A genetic distance threshold of 0.08 (red line) was used to identify unique profiles. Worms for the same host with indistinguishable profiles were pooled and the number of asterisks represents the number of worms for each profile. In total, 16 unique profiles were identified and grouped in three clusters based on unbiased *p*-values (sig. α = 0.05). (b) Dendrogram constructed with EmsB profiles from Europe, Asia and North America (from the EWET database), and the 16 EmsB profiles identified in our study in Alberta (orange and blue for North American and European clades, respectively). In total, 1275 samples from 17 countries were compared with the profiles obtained from our study; thereafter, five profiles with the lowest genetic distance to each of our profiles (P1–P16) were used to build the dendrogram.</figcaption></figure>

(mostly less abundant EmsB profiles) in all subregions (except for AB-CW), and nearly all abundant EmsB profiles had been detected. The areas with the lowest sample coverage ($q = 1$) were AB-NE and AB-CE (89%) and the area with the highest coverage was AB-CW (100%). When Alberta was analysed as a single area, the sample completeness

profile was 100% for all diversity orders (electronic supplementary material, table S3).

The size-based rarefaction and extrapolation analysis revealed that for each subarea the sampling curve stabilized for orders $q = 1$ and 2, but not for $q = 0$ (electronic supplementary material, figure S1), implying that the asymptotic
<figure><img src="image_8.png" /><figcaption>**Figure 3.** NMDS ordination plot of EmsB profiles of *Echinococcus multilocularis* collected in Western Canada from 2012 to 2017. Genetic distance was calculated based on the Euclidean distance. Coloured symbols represent the location (a) and mtDNA haplotype (b); two-dimensional stress = 0.05.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 4.** *Echinococcus multilocularis* samples collected in Alberta from 2012 to 2017, and coverage-based rarefaction and extrapolation sampling curves per subregion in Alberta, Canada. (a) Subregions in Alberta: North-West (AB-NW), North-East (AB-NE), Central-West (AB-CW), Central-West (AB-CE), South-West (AB-SW) and South-East (AB-SE). The dots represent the geographical locations of the parasite samples analysed. (b) Non-asymptotic coverage-based rarefaction (solid lines) and extrapolation (dashed lines) curves for orders q = 0, 1 and 2, (i.e. Chao1 richness, Shannon and inverse Simpson diversity) for each subregion in Alberta. Confidence intervals (95%) are represented by shading. Solid dots denote the observed data points. (AB-NW was excluded from the analysis due to the small sample size).</figcaption></figure>

diversity estimates could be used to infer true diversities (Shannon and Simpson), but not to estimate true genetic richness. When comparing the asymptotic estimates for *q* = 1 and 2 (electronic supplementary material, table S3), the highest diversity indices were found in AB-SW (5.9, 4.5) and the lowest in AB-CW (2.5, 2.4) and AB-SE (3.1, 2.2), with similar values for the other subregions.

Since data were insufficient to infer the true richness, non-asymptotic coverage-based rarefaction and extrapolation curves were calculated (figure 4) and measures were computed up to a standardized coverage value of $C_{max}$ = 94.5% (electronic supplementary material, table S3). For this $C_{max}$, the corresponding highest richness estimate ($q=0$) was found in AB-CE (8.2) and the lowest in AB-CW (2.7). For evenness profiles, all values of *q* showed the highest evenness of EmsB profiles in AB-CW and the lowest in AB-SE.

Analysis of *β*-diversity was performed by constructing a dendrogram based on the Bray–Curtis distance to compare genetic diversity between areas. There was a higher similarity between the central and northern regions, and a second branch grouping for the southern regions (electronic supplementary material, figure S2), with low dissimilarity among all regions (0.27 to 0.53). When comparing non-asymptotic diversity indices from subregions in Alberta with diversity reported in historical endemic areas from Europe (recalculated based on $C_{max}$ = 94.5%), the range of values (min. and max.) for all *q* orders was higher in European countries than in most subregions in Alberta. Moreover, when Alberta was analysed as a single area, estimates for orders *q*=1 and *q*=2 were lower than in all areas in Europe (except for South Germany; electronic supplementary material, table S3). The same pattern was observed for the evenness profiles, with
<figure><img src="image_10.png" /><figcaption>**Figure 5.** Spatial distribution of the EmsB profiles of *Echinococcus multilocularis* collected in Western Canada from 2012 to 2017. (a,b) The distributions of P13 and P15, the two main EmsB profiles in Alberta; (c) the distribution of profiles P4 and P5, present only in Alberta (AB) and Saskatchewan (SK); (d) the distribution of the profiles P6 and P11, present exclusively in BC and AB. (e) The distribution of EmsB profiles per WMU in Alberta. The sizes of the pie charts in (e) represent the numbers of samples collected.</figcaption></figure>

Alberta being the most uneven area (0.2–0.3), which indicates the presence of a dominant EmsB profile and low diversity.

### (d) Spatial distribution of EmsB profiles in Western Canada

Profiles were spatially clustered, with profiles P13 and P15 found exclusively in Alberta (figure 5), and those from BC (P6 and P11) mostly localized only there, except for P11, which was found also in northern AB, close to the BC border. Similarly, in southern SK, only P4 and P5 were detected, with both profiles also present in southern AB. No profiles from BC were detected in SK, or *vice versa*. To assess the hypothesis of genetic isolation by geographical distance, we used a Mantel test comparing a genetic distance matrix between EmsB profiles and a matrix of the geographical distances between samples, based on Euclidean pairwise distances [52]. There was no significant correlation between genetic and geographical distances when testing samples only from AB ($r = 0.02, p = 0.30$). However, the Mantel test yielded a significant correlation ($r = 0.4, p < 0.001$) when including samples from BC and SK in the analysis. The Mantel correlogram indicated a significant positive correlation (more similar samples) within 50 to 550 km, but a negative correlation (more dissimilar samples) when > 600 km apart (electronic supplementary material, figure S3a).

In a dbRDA using the two spatial variables, latitude and longitude, 21.1% of the variance of the genetic matrix was explained by geographical distance (for samples from all three provinces). Still, only the longitude—the axis dividing provinces—was significant ($F = 30.36, p < 0.001$) (electronic supplementary material, figure S3b).

## 4. Discussion

In this study, we investigated the origin of the European-like strains of *E. multilocularis* in Western Canada, their spread

dynamics and potential spatio-temporal scenarios of their process of invasion in North America. We assessed the genetic diversity of the parasite in this region, comparing EmsB profiles and mtDNA sequences; most genetic variants from Western Canada were closely related to the original European clade. Furthermore, the diversity of the European genetic variants found in Western Canada was low compared to that expected for a long-established endemic strain, supporting the hypothesis of relatively recent introductions and discounting the hypothesis of North America as an ancestral endemic focus for this strain. Moreover, the spatial analysis indicated genetic discontinuity, only evident over large geographical scales, suggesting various introduction events in Western Canada.

When using mtDNA as a genetic marker, most haplotypes detected were closely related to the European clade, with few point mutations from European isolates [18,24–26]. However, the number of European samples with sequences available for the three mtDNA genes was limited. Therefore, although not likely, it cannot be ruled out that haplotypes from Western Canada may also be present in European countries but not yet detected. The Bayesian phylogenetic tree implied different evolutionary origins and pathways of the haplotypes found in Western Canada. The ECA, EAB and ESK were grouped in the same branch, whereas BC1 was in a different branch and showed a higher divergence with the closest haplotype SK1. This division was also evident in their geographical distribution, as BC1 was detected in British Columbia, but not in Saskatchewan, whereas ECA and EAB were only in Alberta. These results suggest different spatio-temporal scenarios, with multiple introductions of the European strains to Western Canada, possibly occurring at different times in the past centuries, causing worm populations to be more or less isolated. Indeed, various sources of invasion have been proposed, including the translocation of domestic dogs from European endemic areas [12], introduced European red foxes imported for sport hunting [54], and/or introduction of the parasite via translocation of intermediate hosts with international shipping [9].
During the process of invasion, the likelihood of success-
ful establishment will strongly depend on the propagule
pressure (i.e. the number of introduction events and the
number of infective stages released). Therefore, the repeated
release of a large number of individuals in multiple locations
can facilitate the long-term establishment and increase the
rate of spread of non-native populations [55]. However, mul-
tiple introductions and high genetic variation do not seem to
be indispensable for a successful invasion [56]. For a parasite
like *E. multilocularis*, reproductive traits (i.e. hermaphrodit-
ism, asexual multiplication in intermediate hosts, self- and
cross-fertilization and development of thousands of sexually
mature adult worms in each definitive host) could have also
been determinant in the first stages of invasion, boosting local
spread, but producing in turn geographical isolation of the
different genotypes introduced.

The EmsB microsatellite enabled detection of more profiles (16 in total) in comparison to only seven mtDNA haplotypes, helping to identify local genetic variation and trace the invasion process of the parasite. When compared to a world database, most of the profiles we found (13 of 16) clustered within the European group, having a low genetic distance to European isolates (ranging from 0.05 to 0.1). Moreover, half of the profiles were genetically indistinguishable from known European ones, supporting the hypothesis of a recent introduction of these European genetic variants. By contrast, the remaining three EmsB profiles (P1-P3), which were genotyped as North American haplotype N2 based on mtDNA, clustered within the Arctic clade and were genetically more distant from the Alaskan and Svalbard isolates (0.21 to 0.42), suggesting a more prolonged genetic isolation, consistent with the greater divergence between the haplotypes N1 and N2 [20,57].

Complete agreement between the two markers, in correlation with the geographical location, is not guaranteed. Umhang et al. [22] described samples from Poland clustering within the European clade, based on mtDNA, whereas they were genotyped as Asian-type profiles using EmsB, suggesting some degree of introgression between the two strains. Regardless, these two markers facilitated differentiation between North American and European genetic variants, with no evidence of introgression between them. All EmsB profiles clustered within the European clade were identified as European-like haplotypes based on mtDNA, whereas EmsB profiles identified as North American strains were in a different cluster. Moreover, most mtDNA haplotypes were represented by more than one EmsB profile from the same cluster, with samples from Saskatchewan clustering apart from samples from British Columbia, as observed with mtDNA haplotypes.

In Europe, the expansion history of the parasite from the
historical endemic core area into peripheral regions has been
governed by a mainland-island system of transmission, in
which the ancestral focus in central Europe served as a 'main-
land' supplying the peripheral areas (islands), perhaps due to
dispersals of fox populations. From a genetic perspective, this
resulted in invasion events with introduction of only a few
genetic profiles, resulting in low genetic diversity in the colo-
nized region compared to greater diversity in the ancestral
endemic foci [33,50,58]. In our study, the genetic diversity
indices recalculated for Switzerland (endemic area), based
on C<sub>max</sub> = 94.5%, were as high as 15 (q = 1) and 11.5 (q = 2),
with an estimated richness of 23 profiles. In comparison, in
the province of Alberta (with 661 848 km<sup>2</sup> and almost twice
as big as Germany), only 15 EmsB profiles were found,

with genetic diversity indices between 2.3 to 5.1 (q = 1) and
2.1 to 4.1 (q = 2) in the five subareas, indicating low genetic
diversity. This is inconsistent with the hypothesis of an unde-
tected historical endemic area of the European-like clade.
Moreover, the genetic indices recalculated for Switzerland
were more than three times higher than those calculated for
all of Alberta.

In this Canadian province, the highest genetic richness was found in the central east region, and the highest genetic diversity based on the abundance of each EmsB profile was found in the southwestern region. However, the degree of differentiation between *E. multilocularis* populations from all five subregions was low, with widespread distribution of two profiles (P13 and P15) representing 77% of the samples from this province. Moreover, these two profiles might be considered as part of one population, due to: (i) the low genetic distance between them (0.14), (ii) their overlap in geographical distribution and (iii) the relatively high frequency of co-occurrence in the same host (19 hosts). These results were consistent with previous studies using mtDNA showing a higher prevalence and wider geographical distribution of the haplotype ECA, the causative agent of the most recent human cases of AE in Alberta [25,26]. Therefore, the presence of a single predominant genetic variant in Alberta supported the hypothesis of a single invasion event responsible for the initial establishment of a small number of *E. multilocularis* individuals in that province.

The subsequent spread across the province was likely aided
by red foxes, but mostly by the most abundant host, the coyote,
which has a larger home range (greater than 100 km²) and
higher dispersal distances (up to 300 km) [59–62]. In a previous
study comparing coyotes and red foxes from Alberta infected
with the ECA haplotype, the difference in intensity of infection
(worm burden) between both hosts was significantly higher in
coyotes than foxes [26], which could be related to a lack of co-
evolved resistance in coyotes to the European strains. There-
fore, the coyote, being a naïve host, might be the primary
source of environmental contamination with eggs of these
strains. In this study, we did not perform an analysis of intra-
host diversity to compare red foxes and coyotes due to the
low sample size of red foxes. However, we observed a higher
prevalence of the profile P13 (identified as the ECA haplotype
based on mtDNA) in both coyotes and red foxes, and a low
prevalence of North American profiles, which is consistent
with previous studies using deep amplicon sequencing [26].
Further analyses of differences in the prevalence and intensity
of infection of European and North American genetic variants
in coyotes compared to red foxes will be important to under-
stand the role of the two hosts in the transmission and
spread of the European strains and the intra-host competition
between North American and European strains.

To estimate the expansion rate of the parasite, Takumi et al.
[63] used a mathematical model predicting a spreading rate of
2.7 km per year, from an endemic to a non-endemic area in
the Netherlands. By contrast, the presence of highly vagile
hosts, like coyotes in Western Canada, could have facilitated
and accelerated the spread of the European strain. Thus, it is
likely that the process of colonizing the entire province of
Alberta was completed in a shorter time frame, expanding
from the south to the north, as evidenced by the highest differ-
entiation between the Southern and the Northern regions of
Alberta, and the distribution of the least common profiles in dis-
tant regions. This pattern, with one single cluster of profiles
being the most prevalent, and low genetic diversity, was similar
to that found in European regions initially not considered endemic. For example, in a survey conducted in Poland, only one profile was the most prevalent and was distributed across the country [33]. Similarly, in Denmark and Sweden, one profile represented 68.4 and 55.5% of worms, respectively [44].

In a previous study using EmsB to assess the genetic diversity of the parasite at a continental scale in Europe [50], only 5% of the genetic variability was explained by geographical distance. In our study, there was no significant correlation when analysing samples from Alberta alone. However, there was a significant correlation when including samples from British Columbia and Saskatchewan, with more dissimilar samples that were far apart (greater than 600 km). In the distance-based redundancy analysis, 21.1% of the genetic variability was explained by geographical distance, with longitude being the only significant spatial variable. Our results could indicate independent invasion events in each province and isolation by distance, with the exchange of some genetic variants between neighbouring regions, thanks to the dispersal movement of definitive hosts. Furthermore, these invasion events were likely influenced by geographical barriers, e.g. the Rocky Mountains between Alberta and British Columbia.

Despite the smaller sample size in Saskatchewan and British Columbia, our results were consistent with previous studies using mtDNA characterization in these areas. In British Columbia, Gesy et al. [18] reported only one European-like haplotype (BC1) in the immediate area of Quesnel, BC. Likewise, in a study in Saskatchewan, only one European-like haplotype (SK1) was detected in the central area of the province. Yet, seven haplotypes that belonged to the North American clade were detected in the southern region [24]. However, in our study in Alberta, the abundance of the North American strain was surprisingly low, compared to the European strains. Nonetheless, the genetic variability of the North American strain (i.e. number of disparate profiles in a few hosts) was high (three EmsB profiles in eight worms from two hosts), which is consistent with what would be expected from an ancestral endemic strain. The difference in the prevalence, genetic diversity and distribution between the North American and the European strains observed in Alberta, using nuclear and mitochondrial markers, was relevant, suggesting some degree of competitive interactions, with a population expansion of the European over the North American strain [26].

This was the first large-scale study using both nuclear and mitochondrial markers to assess the genetic diversity of *E. multilocularis* in North America. The use of EmsB improved characterization of parasite diversity, even at a fine geographical scale, whereas the use of mtDNA aided in unravelling the evolutionary process of the invasion of the European strain. As evidenced in this and previous studies, the implementation of genetic population structure analysis is a powerful tool to trace the origins and history expansion of parasites and their hosts in the areas invaded and for the assessment of public health risks. For example, a genetic study of raccoons (*Procyon lotor*) in the Netherlands (where they were previously absent) and their GI nematode parasite

(*Baylisascaris procyonis*) showed that most of the Dutch raccoons and their roundworms were introduced through ex-captive individuals, which ultimately aided in the development of control measures of these invasive populations [64]. Recently, a genomic analysis allowed the evaluation of the process of colonization of the Americas by *Schistosoma mansoni*, a blood fluke that infects humans and that was introduced into the Americas from Africa during the Transatlantic slave trade [65]. In this study, no evidence of population bottlenecks was observed, suggesting that *S. mansoni* parasites were pre-adapted to the Americas and able to establish with relative ease.

Further work on the study of genetic diversity of *E. multilocularis* including wild and domestic hosts from British Columbia and Saskatchewan, as well as from eastern Canadian provinces where the parasite has been recently detected, and in central US states, will be pivotal to understanding the current distribution and expansion trend of the European strains in continental North America, including the role of the main definitive hosts (coyotes and foxes) in harbouring and spreading these strains.

**Data accessibility.** All data generated or analysed during this study are included in this article (and its supplementary information files). The dataset of EmsB profiles is available in the Figshare repository: (https://doi.org/10.6084/m9 figshare.16818607) [66].

The data are provided in the electronic supplementary material [67].

**Authors' contributions.** M.A.S.: conceptualization, formal analysis, investigation, methodology, visualization, writing—original draft; G.U.: data curation, formal analysis, investigation, methodology, visualization and writing—review and editing; C.K.: investigation, methodology and writing—review and editing; D.M.G.: data curation, investigation and writing—review and editing; K.E.R.: funding acquisition, supervision and writing—review and editing; M.M.: resources, supervision and writing—review and editing; J.G.: writing—review and editing; A.M.: conceptualization, funding acquisition, investigation, methodology, supervision, visualization and writing—original draft.

All authors gave final approval for publication and agreed to be held accountable for the work performed therein.

**Conflict of interest declaration.** We declare we have no competing interests.

**Funding.** MITACS (Mathematics of Information Technology and Complex Systems), Accelerate Program, grant no. 10018836 (to A.M.); ELANCO Canada Limited, grant no. 10017067 (to A.M.); ACA (Alberta Conservation Association), grant no. 10010809 (to A.M.); Natural Sciences & Engineering Research Council of Canada (NSERC), Discovery Grant, Grant no.10013904, URGC SEED Grant; grant no.10007230 (to A.M.); NSERC, Discovery Grant, grant no. RGPIN-2018-03913 (to K.E.R.); City of Calgary, grant no.10007719 (to A.M.).

**Acknowledgements.** We gratefully acknowledge the ATA (Alberta Trappers' Association), the trappers who provided animal carcasses, Vanessa Bastid (Nancy Laboratory for Rabies and Wildlife) who helped with the molecular analysis, Dimitri Giunchi (University of Pisa) who helped with the Bayesian analysis, Maria Cavedon and Jessica Theoret (University of Calgary) who helped with spatial analysis, John Patrick Kastelic (University of Calgary) for revising the final draft of the paper, and all undergraduate students that contributed to laboratory work.

## References

1. Sax DF, Gaines SD. 2003 Species diversity: from global decreases to local increases. *Trends Ecol.*

*Evol.* **18**, 561–566. (doi:10.1016/S0169-5347(03)00224-6)

2. Young HS, Parker IM, Gilbert GS, Sofia Guerra A, Nunn CL. 2017 Introduced species, disease ecology,
and biodiversity-disease relationships. *Trends Ecol. Evol.* **32**, 41–54. (doi:10.1016/j.tree.2016.09.008)

3. Hatcher MJ, Dick JTA, Dunn AM. 2012 Disease emergence and invasions. *Funct. Ecol.* **26**, 1275–1287. (doi:10.1111/j.1365-2435.2012.02031.x)

4. Ogden NH et al. 2019 Emerging infectious diseases and biological invasions: a call for a One Health collaboration in science and management. *R. Soc. Open Sci.* **6**, 181577. (doi:10.1098/rsos.181577)

5. Thompson R.CA, Deplazes P, Lymbery AJ. 2017 *Advances in parasitology: Echinococcus and Echinococcosis*, part A, 1st edn. New York, NY: Academic Press.

6. Eckert J et al. 2001 *WHO/OIE manual on echinococcosis in humans and animals: a public health problem of global concern*. Paris, France: World Organization for Animal Health.

7. Thompson RCA. 2017 Biology and systematics of *Echinococcus*. *Adv. Parasitol.* **95**, 65–109. (doi:10.1016/bs.apar.2016.07.001)

8. Torgerson PR, Robertson LJ, Enemark HL, Foehr J, Van Der Giessen JWB, Kapel CMO, Klun I, Trevisan C. 2020 Source attribution of human echinococcosis: a systematic review and meta-analysis. *PLoS Negl. Trop. Dis* **14**, e0008382. (doi:10.1371/journal.pntd.0008382)

9. Davidson RK, Romig T, Jenkins E, Tryland M, Robertson LJ. 2012 The impact of globalisation on the distribution of *Echinococcus multilocularis*. *Trends Parasitol.* **28**, 239–247. (doi:10.1016/j.pt.2012.03.004)

10. Eckert J, Conraths FJ, Tackmann K. 2000 Echinococcosis: an emerging or re-emerging zoonosis? *Int. J. Parasitol.* **30**, 1283–1294. (doi:10.1016/S0020-7519(00)00130-2)

11. Giraudoux P et al. 2003 Interactions between landscape changes and host communities can regulate *Echinococcus multilocularis* transmission. *Parasitology* **127**, S121–S131. (doi:10.1017/S0031182003003512)

12. Jenkins EJ, Peregrine AS, Hill JE, Somers C, Gesy K, Barnes B, Gottstein B, Polley L. 2012 Detection of European strain of *Echinococcus multilocularis* in North America. *Emerg. Infect. Dis.* **18**, 1010–1012. (doi:10.3201/eid1806.111420)

13. Gottstein B, Stojkovic M, Vuitton DA, Millon L, Marcinkute A, Deplazes P. 2015 Threat of alveolar echinococcosis to public health: a challenge for Europe. *Trends Parasitol.* **31**, 407–412. (doi:10.1016/j.pt.2015.06.001)

14. Massolo A et al. 2018 Unexpected *Echinococcus multilocularis* infections in shepherd dogs and wolves in south-western Italian Alps: a new endemic area? *Int. J. Parasitol. Parasites Wildl.* **7**, 309–316. (doi:10.1016/j.ijppaw.2018.08.001)

15. Kotwa J, Isaksson M, Jardine C, Douglas Campbell G, Berke O, Pearl DL, Mercer NJ, Osterman-Lind E, Peregrine A. 2019 *Echinococcus multilocularis* infection, Southern Ontario, Canada. *Emerg. Infect. Dis.* **25**, 265–272. (doi:10.3201/eid2502.180299)

16. Deplazes P, Hagglin D, Gloor S, Romig T. 2004 Wilderness in the city: the urbanization of *Echinococcus multilocularis*. *Trends Parasitol.* **20**, 77–84. (doi:10.1016/j.pt.2003.11.011)

17. Liccioli S, Giraudoux P, Deplazes P, Massolo A. 2015 Wilderness in the 'city' revisited: different urbes shape transmission of *Echinococcus multilocularis* by altering predator and prey communities. *Trends Parasitol.* **31**, 297–305. (doi:10.1016/j.pt.2015.04.007)

18. Gesy K, Hill JE, Schwantje H, Liccioli S. 2013 Establishment of a European-type strain of *Echinococcus multilocularis* in Canadian wildlife. *Parasitology* **140**, 1133–1137. (doi:10.1017/S0031182013000607)

19. Henttonen H, Fuglei E, Gower CN, Haukisalmi V, Ims RA, Niemimaa J, Yoccoz NG. 2001 *Echinococcus multilocularis* on Svalbard: introduction of an intermediate host has enabled the local life-cycle. *Parasitology* **123**(Pt.6), 547–552. (doi:10.1017/S0031182001008800)

20. Nakao M, Xiao N, Okamoto M, Yanagida T, Sako Y, Ito A. 2009 Geographic pattern of genetic variation in the fox tapeworm *Echinococcus multilocularis*. *Parasitol. Int.* **58**, 384–389. (doi:10.1016/j.parint.2009.07.010)

21. Konyaev SV et al. 2013 Genetic diversity of *Echinococcus spp.* in Russia. *Parasitology* **140**, 1637–1647. (doi:10.1017/S0031182013001340)

22. Umhang G, Knapp J, Wassermann M, Bastid V, Peytavin de Garam C, Boué F, Cencek T, Romig T, Karamon J. 2021 Asian admixture in European *Echinococcus multilocularis* populations: new data from Poland comparing EmsB microsatellite analyses and mitochondrial sequencing. *Front. Vet. Sci.* **7**, 620722. (doi:10.3389/fvets.2020.620722)

23. Karamon J, Stojlecki K, Samorek-Pierog M, Bilskazajac E, Rozycki M, Chmurzynska E, Sroka J, Zdybel J, Cencek T. 2017 Genetic diversity of *Echinococcus multilocularis* in red foxes in Poland: the first report of a haplotype of probable Asian origin. *Folia Parasitol. (Praha)* **64**, 7. (doi:10.14411/fp.2017.007)

24. Gesy KM, Jenkins EJ. 2015 Introduced and native haplotypes of *Echinococcus multilocularis* in wildlife in Saskatchewan, Canada. *J. Wildl. Dis.* **51**, 743–748. (doi:10.7589/2014-08-214)

25. Massolo A et al. 2019 European *Echinococcus multilocularis* identified in patients in Canada. *N. Engl. J. Med.* **381**, 384–385. (doi:10.1056/NEJMCM184975)

26. Santa MA, Rezansoff AM, Chen R, Gilleard JS, Musiani M, Ruckstuhl KE, Massolo A. 2021 Deep amplicon sequencing highlights low intra-host genetic variability of *Echinococcus multilocularis* and high prevalence of the European-type haplotypes in coyotes and red foxes in Alberta, Canada. *PLoS Negl. Trop. Dis.* **15**, e0009428. (doi:10.1371/journal.pntd.0009428)

27. Houston S et al. 2021 Epidemiological and clinical characteristics of alveolar Echinococcosis: an emerging infectious disease in Alberta, Canada. *Am. J. Trop. Med. Hyg.* **104**, 1863–1869. (doi:10.4269/ajtmh.20-1577)

28. Schurer JM, Tsybina P, Gesy KM, Kolapo TU, Skinner S, Hill JE, Jenkins EJ. 2021 Molecular evidence for local acquisition of human alveolar Echinococcosis in

Saskatchewan, Canada. *J. Infect. Dis.* **223**, 1015–1018. (doi:10.1093/infdis/jiaa473)

29. Polish LB et al. 2022 European haplotype of *Echinococcus multilocularis* in the United States. *N. Engl. J. Med.* **387**, 1902–1904. (doi:10.1056/nejmc2210000)

30. Polish LB, Pritt B, Barth TFE, Gottstein B, O'Connell EM, Gibson PC. 2020 First European haplotype of *Echinococcus multilocularis* identified in the United States: an emerging disease? *Clin. Infect. Dis.* **72**, 1117–1123. (doi:10.1093/cid/ciaa245)

31. Umhang G, Knapp J, Hormaz V, Raoul F, Boué F. 2014 Using the genetics of *Echinococcus multilocularis* to trace the history of expansion from an endemic area. *Infect. Genet. Evol.* **22**, 142–149. (doi:10.1016/j.meegid.2014.01.018)

32. Knapp J, Bart JM, Maillard S, Gottstein B, Piarroux R. 2010 The genomic *Echinococcus microsatellite* EmsB sequences: from a molecular marker to the epidemiological tool. *Parasitology* **137**, 439–449. (doi:10.1017/S0031182009991612)

33. Umhang G, Karamon J, Hormaz V, Knapp J, Cencek T, Boué F. 2017 A step forward in the understanding of the presence and expansion of *Echinococcus multilocularis* in Eastern Europe using microsatellite EmsB genotyping in Poland. *Infect. Genet. Evol.* **54**, 176–182. (doi:10.1016/j.meegid.2017.07.004)

34. Knapp J, Guislain MH, Bart JM, Raoul F, Gottstein B, Giraudoux P, Piarroux R. 2008 Genetic diversity of *Echinococcus multilocularis* on a local scale. *Infect. Genet. Evol.* **8**, 367–373. (doi:10.1016/j.meegid.2008.02.010)

35. Santa MA, Pastran SA, Klein C, Duignan P, Ruckstuhl K, Romig T, Massolo A. 2018 Detecting co-infections of *Echinococcus multilocularis* and *Echinococcus canadensis* in coyotes and red foxes in Alberta, Canada using real-time PCR. *Int. J. Parasitol. Parasites Wildl.* **7**, 111–115. (doi:10.1016/j.ijppaw.2018.03.001)

36. Gesy K, Pawlik M, Kapronczai L, Wagner B, Elkin B, Schwantje H, Jenkins E. 2013 An improved method for the extraction and quantification of adult *Echinococcus* from wildlife definitive hosts. *Parasitol. Res.* **112**, 2075–2078. (doi:10.1007/s00436-013-3371-x)

37. Paradis E. 2010 pegas: an R package for population genetics with an integrated-modular approach. *Bioinformatics* **26**, 419–420. (doi:10.1093/bioinformatics/btp696)

38. Ronquist F et al. 2012 MrBayes 3.2: efficient Bayesian phylogenetic inference and model choice across a large model space. *Syst. Biol.* **61**, 539–542. (doi:10.1093/sysbio/sys029)

39. Schliep KP. 2011 phangorn: phylogenetic analysis in R. *Bioinformatics* **27**, 592–593. (doi:10.1093/bioinformatics/btq706)

40. Knapp J et al. 2017 EWET: Data collection and interface for the genetic analysis of *Echinococcus multilocularis* based on EmsB microsatellite. *PLoS ONE* **12**, e0183849. (doi:10.1371/journal.pone.0183849)

41. Legendre P, Legendre L. 2012 *Numerical ecology*. Oxford, UK: Elsevier.
42. Shimodaira H. 2004 Approximately unbiased tests of regions using multistep-multiscale bootstrap resampling. *Ann. Stat.* **32**, 2616–2641. (doi:10.1214/009053604000000823)

43. Knapp J, Bart JM, Glowatzki ML, Ito A, Gerard S, Maillard S, Piarroux R, Gottstein B. 2007 Assessment of use of microsatellite polymorphism analysis for improving spatial distribution tracking of *Echinococcus multilocularis*. *J. Clin. Microbiol.* **45**, 2943–2950. (doi:10.1128/JCM.02107-06)

44. Knapp J, Umhang G, Wahlström H, Al-Sabi M.NS, Ågren EO, Enemark HL. 2019 Genetic diversity of *Echinococcus multilocularis* in red foxes from two Scandinavian countries: Denmark and Sweden. *Food Waterborne Parasitol.* **14**, e00045. (doi:10.1016/j.fawpar.2019.e00045)

45. Kruskal JB. 1964 Nonmetric multidimensional scaling: a numerical method. *Psychometrika* **29**, 115–129. (doi:10.1007/bf02289694)

46. Suzuki R, Shimodaira H. 2006 Pvclust: an R package for assessing the uncertainty in hierarchical clustering. *Bioinformatics* **22**, 1540–1542. (doi:10.1093/bioinformatics/btl117)

47. Chao A et al. 2020 Quantifying sample completeness and comparing diversities among assemblages. *Ecol. Res.* **35**, 292–314. (doi:10.1111/1440-1703.12102)

48. Chao A, Gotelli NJ, Hsieh TC, Sander EL, Ma KH, Colwell RK, Ellison AM. 2014 Rarefaction and extrapolation with Hill numbers: a framework for sampling and estimation in species diversity studies. *Ecol. Monogr.* **84**, 45–67. (doi:10.1890/13-0133.1)

49. Chao A, Jost L. 2012 Coverage-based rarefaction and extrapolation: standardizing samples by completeness rather than size. *Ecology* **93**, 2533–2547. (doi:10.1890/11-1952.1)

50. Knapp J et al. 2009 Genetic diversity of the Cestode *Echinococcus multilocularis* in red foxes at a Continental Scale in Europe. *PLoS Negl. Trop. Dis.* **3**, e452. (doi:10.1371/journal.pntd.0000452)

51. Oksanen J et al. 2013 vegan:Community ecology package. *R package version*, 2.0-2. Software http://CRAN.R-project.org/package=vegan.

52. Mantel N. 1967 The detection of disease clustering and a generalized regression approach. *Cancer Res.* **27**, 209–220.

53. Legendre P, Anderson MJ. 1999 Distance-based redundancy analysis: testing multispecies responses in multifactorial ecological experiments. *Ecol. Monogr.* **69**, 1–24. (doi:10.1890/0012-9615(1999)069[0001:DBRATM]2.0.CO;2)

54. Kamler JF, Ballard WB. 2002 A review of native and nonnative red foxes in North America. *Wildl. Soc. Bull.* (1973–2006) **30**, 370–379. (https://www.jstor.org/stable/3784493)

55. Lockwood JL, Cassey P, Blackburn T. 2005 The role of propagule pressure in explaining species invasions. *Trends Ecol. Evol.* **20**, 223–228. (doi:10.1016/j.tree.2005.02.004)

56. Dlugosch KM, Anderson SR, Braasch J, Cang FA, Gillette HD. 2015 The devil is in the details: genetic variation in introduced populations and its contributions to invasion. *Mol. Ecol.* **24**, 2095–2111. (doi:10.1111/mec.13183)

57. Davidson RK, Lavikainen A, Konyaev S, Schurer J, Miller AL, Oksanen A, Skirnisson K, Jenkins E. 2016 *Echinococcus* across the north: current knowledge, future challenges. *Food Waterborne Parasitol.* **4**, 39–53. (doi:10.1016/j.fawpar.2016.08.001)

58. Umhang G et al. 2021 Unravelling the genetic diversity and relatedness of *Echinococcus multilocularis* isolates in Eurasia using the EmsB microsatellite nuclear marker. *Infect. Genet. Evol.* **92**, 104863. (doi:10.1016/j.meegid.2021.104863)

59. Gosselink TE, Van Deelen TR, Warner RE, Joselyn MG. 2003 Temporal habitat partitioning and spatial use of coyotes and red foxes in East-Central Illinois. *J. Wildl. Manage.* **67**, 90–103. (doi:10.2307/3803065)

60. Dlugosch KM, Parker IM. 2008 Founding events in species invasions: genetic variation, adaptive evolution, and the role of multiple introductions. *Mol. Ecol.* **17**, 431–449. (doi:10.1111/j.1365-294X.2007.03538.x)

61. Gese EM, Rongstad OJ, Mytton WR. 1988 Home range and habitat use of coyotes in Southeastern Colorado. *J. Wildl. Manage.* **52**, 640–646. (doi:10.2307/3800923)

62. Harrison DJ. 1992 Dispersal characteristics of juvenile coyotes in Maine. *J. Wildl. Manage.* **56**, 128–138. (doi:10.2307/3808800)

63. Takumi K, de Vries A, Chu ML, Mulder J, Tenuis P, van der Giessen J. 2008 Evidence for an increasing presence of *Echinococcus multilocularis* in foxes in the Netherlands. *Int. J. Parasitol.* **38**, 571–578. (doi:10.1016/j.ijpara.2007.09.014)

64. Maas M et al. 2022 Population genetics, invasion pathways and public health risks of the raccoon and its roundworm *Baylisascaris procyonis* in northwestern Europe. *Transbound. Emerg. Dis.* **69**, 2191–2200. (doi:10.1111/tbed.14218)

65. Platt RN et al. 2022 Genomic analysis of a parasite invasion: colonization of the Americas by the blood fluke *Schistosoma mansoni*. *Mol. Ecol.* **31**, 2242–2263. (doi:10.1111/mec.16395)

66. Santa MA, Umhang G, Klein C, Grant DM, Ruckstuhl KE, Musiani M, Gilleard JS, Massolo A. 2023 It's a small world for parasites: evidence supporting the North American invasion of European *Echinococcus multilocularis*. Figshare. (doi:10.6084/m9 figshare.16818607)

67. Santa MA, Umhang G, Klein C, Grant DM, Ruckstuhl KE, Musiani M, Gilleard JS, Massolo A. 2023 It's a small world for parasites: evidence supporting the North American invasion of European *Echinococcus multilocularis*. Figshare. (doi:10.6084/m9 figshare.c.6440232)