INTRODUCTION
Echinococcosis or hydatidosis is a serious parasitic zoonosis which affects various internal organs of terrestrial mammals including humans, livestock and wildlife. Pathogenic organisms of the genus Echinococcus (Cestoda: Taeniidae) require an obligatory mammalian host. Canids and felids are listed as definitive hosts for adult tapeworms and many ungulates, lagomorphs and rodents as intermediate hosts for metacestode larvae. Humans are involved as accidental intermediate hosts by ingesting eggs derived from definitive hosts. In metacestode infections, affected organs show three types of macropathological lesion known as unilocular (unicystic), polycystic and alveolar hydatids. The fox tapeworm E. multilocularis in the Holarctic region is responsible for alveolar hydatid, while polycystic hydatid is caused by Neotropical E. vogeli (Rausch and Bernstein, Reference Rausch and Bernstein1972; D'Alessandro and Rausch, Reference D'Alessandro and Rausch2008). Metacestodes of Neotropical E. oligarthrus and Tibetan E. shiquicus develop into unicystic hydatids in rodents and lagomorphs, respectively (Xiao et al. Reference Xiao, Qiu, Nakao, Li, Yang, Chen, Schantz, Craig and Ito2005; D'Alessandro and Rausch, Reference D'Alessandro and Rausch2008). However, in large mammals including humans, it is difficult to distinguish morphologically the causative agents of unilocular hydatids.
Historically, the type species E. granulosus had been considered as an agent of human echinococcosis causing both unilocular and alveolar hydatids up until the 1950s (Rausch, Reference Rausch and Thompson1986; Tappe and Frosch, Reference Tappe and Frosch2007). The involvement of E. multilocularis as a human pathogen was confirmed by Vogel (Reference Vogel1957), based on the morphological observation of adult tapeworms originated from the alveolar lesions of human cases. Even after this evaluation, the taxonomy of unilocular hydatids and their adult worms was still controversial because of the presence of morphologically similar species and subspecies differentiated based mainly on host specificity (Sweatman and Williams, Reference Sweatman and Williams1963; Williams and Sweatman, Reference Sweatman and Williams1963; Verster, Reference Verster1965). Since unsuitable taxa were invalidated or synonymized into E. granulosus, only four morphospecies of E. granulosus, E. multilocularis, E. oligarthrus and E. vogeli were recognized as valid until recently (Rausch and Nelson, Reference Rausch and Nelson1963; Rausch, Reference Rausch1967; Rausch and Bernstein, Reference Rausch and Bernstein1972; Eckert et al. Reference Eckert, Gemmell, Meslin and Pawłowski2001). The taxonomic simplification led to the recognition of ‘strains’ of E. granulosus, i.e. epidemiologically distinct entities using particular ungulates as intermediate hosts (Thompson et al. Reference Thompson, Lymbery and Constantine1995; Eckert et al. Reference Eckert, Gemmell, Meslin and Pawłowski2001).
The pioneering works of molecular taxonomy of Echinococcus (Bowles et al. Reference Bowles, Blair and McManus1992, Reference Bowles, Blair and McManus1994; Bowles and McManus, Reference Bowles and McManus1993) certified the four morphospecies and defined the genotypes G1 to G8 of E. granulosus according to the variations of partial mitochondrial DNA (mtDNA) sequences of genes for cytochrome c oxidase subunit 1 (cox1) and NADH dehydrogenase subunit 1 (nad1). These genotypes corresponded to the biological strains as follows: G1 (sheep strain), G2 (Tasmanian sheep strain), G3 (buffalo strain), G4 (horse strain), G5 (cattle strain), G6 (camel strain), G7 (pig strain) and G8 (cervid strain). However, subsequent phylogenetic study using the data sets of mtDNA (cox1 and nad1) and nuclear DNA (ribosomal internal transcribed spacer 1) suggested that E. granulosus genotypes are paraphyletic (Bowles et al. Reference Bowles, Blair and McManus1995). Nevertheless, the genotypes and strains were continuously used for epidemiological and ecological purposes. Subsequently, the genotypes of E. granulosus G9 (unidentified strain) and G10 (Fennoscandian cervid strain) were discovered (Scott et al. Reference Scott, Stefaniak, Pawlowski and McManus1997; Lavikainen et al. Reference Lavikainen, Lehtinen, Meri, Hirvelä-Koski and Meri2003), but the status of G9 is still uncertain (Kedra et al. Reference Kedra, Swiderski, Tkach, Dubinsky, Pawlowski, Stefaniak and Pawlowski1999). Recent field surveys denied the endemicity of G8 and G10 (Thompson et al. Reference Thompson, Boxell, Ralston, Constantine, Hobbs, Shury and Olson2006; Moks et al. Reference Moks, Jõgisalu, Valdmann and Saarma2008; Himsworth et al. Reference Himsworth, Jenkins, Hill, Nsungu, Ndao, Thompson, Covacin, Ash, Wagner, McConnell, Leighton and Skinner2010; Bryan et al. Reference Bryan, Darimont, Hill, Paquet, Thompson, Wagner and Smits2012). Both genotypes occur sympatrically in the boreal forests of the Holarctic region.
Mitochondrial genome (mt genome) sequencing of E. granulosus G1 and G4 genotypes (Le et al. Reference Le, Pearson, Blair, Dai, Zhang and McManus2002) and E. multilocularis (Nakao et al. Reference Nakao, Yokoyama, Sako, Fukunaga and Ito2002) accelerated the taxonomic revisions of E. granulosus genotypes. The first revision elevated E. granulosus equinus (Williams and Sweatman, Reference Sweatman and Williams1963) to the species rank and resurrected E. ortleppi (Thompson and McManus, Reference Thompson and McManus2002). The genotypes G1, G2 and G3 (sheep and buffalo strains) were, therefore, treated as E. granulosus sensu stricto (s.s.), G4 (horse strain) as E. equinus and G5 (cattle strain) as E. ortleppi. The second revision, based on mt genome sequences (Nakao et al. Reference Nakao, McManus, Schantz, Craig and Ito2007), completely supported the first revision, and E. granulosus canadensis (Cameron, Reference Cameron1960; Sweatman and Williams, Reference Sweatman and Williams1963) was further elevated to the species rank. The genotypes G6, G7 and G8 (camel, pig and cervid strains) were unified into E. canadensis because of their close genetic relationships. Although the mt genome sequence of G10 was still undetermined, this genotype was included in E. canadensis by considering mtDNA genealogy published previously (Lavikainen et al. Reference Lavikainen, Lehtinen, Meri, Hirvelä-Koski and Meri2003). The second revision also showed that E. shiquicus, a newly found species in the Tibetan Plateau (Xiao et al. Reference Xiao, Qiu, Nakao, Li, Yang, Chen, Schantz, Craig and Ito2005), is sister to E. multilocularis.
These revisions have widely been accepted, but the species status of E. canadensis is still controversial (Cruz-Reyes et al. Reference Cruz-Reyes, Constantine, Boxell, Hobbs and Thompson2007; Thompson, Reference Thompson2008) because biologically different strains were unified. Based mainly on the specimens from North America, Sweatman and Williams (Reference Sweatman and Williams1963) split the cervid strain into two subspecies, E. granulosus borealis and E. granulosus canadensis, occurring, respectively, in the wolf–moose sylvatic cycle and wolf/dog-reindeer semi-domestic cycle. The first was thought to be indigenous to boreal forests of North America, and the latter to be introduced from Norway to arctic north-western Canada (Sweatman and Williams, Reference Sweatman and Williams1963). The genotypes G8 and G10 were speculated to correspond to E. granulosus borealis and E. granulosus canadensis, respectively (Lavikainen, Reference Lavikainen, Lehtinen, Meri, Hirvelä-Koski and Meri2003), but the intermediate host specificity and geographic distribution of the genotypes overlap remarkably (Lavikainen et al. Reference Lavikainen, Lehtinen, Laaksonen, Agren, Oksanen and Meri2006; Thompson et al. Reference Thompson, Boxell, Ralston, Constantine, Hobbs, Shury and Olson2006; Moks et al. Reference Moks, Jõgisalu, Valdmann and Saarma2008). Comparative pathology of human cystic echinococcosis in boreal and arctic zones suggests malignant and benign disease types probably due to G8 and G10, respectively (Wilson et al. Reference Wilson, Diddams and Rausch1968; McManus et al. Reference McManus, Zhang, Castrodale, Le, Pearson and Blair2002). Rausch (Reference Rausch2003) denied the subspecies because of their sympatric occurrence and treated them as the northern biotype of E. granulosus. By contrast, the camel and pig strains (the genotypes G6 and G7) are maintained primarily by dog-mediated domestic life cycles from tropical to temperate zones, indicating that G6 and G7 are ecologically and geographically segregated from G8 and G10. Therefore, a taxonomic opinion that G6 and G7 should be treated as a single species, E. intermedius, has been presented (Cruz-Reyes et al. Reference Cruz-Reyes, Constantine, Boxell, Hobbs and Thompson2007; Thompson, Reference Thompson2008). The resurrection of this old specific name seems to be inappropriate (Nakao et al. Reference Nakao, Yanagida, Okamoto, Knapp, Nkouawa, Sako and Ito2010b), but the idea would be supported if G6 and G7 were sister to G8 and G10 (Fig. 1). However, molecular surveys using short mtDNA sequences suggested that G8 and G10 are paraphyletic and G10 is more closely related to G6 and G7 (Lavikainen et al. Reference Lavikainen, Lehtinen, Laaksonen, Agren, Oksanen and Meri2006; Moks et al. Reference Moks, Jõgisalu, Valdmann and Saarma2008).
Schema showing possible relationships among E. canadensis genotypes/strains. Closed circles in cladograms show possible ancestors susceptible to cervids.
The previously published phylogeny based on mt genome sequences (Nakao et al. Reference Nakao, McManus, Schantz, Craig and Ito2007) lacked two important taxa of E. felidis and E. canadensis G10. The former African species had been included in E. granulosus sensu lato as a subspecies or strain, but its validity was confirmed recently (Hüttner et al. Reference Hüttner, Nakao, Wassermann, Siefert, Boomker, Dinkel, Sako, Mackenstedt, Romig and Ito2008). In this study, the mitochondrial phylogeny is revised by adding these taxa, particularly to examine interrelationships among E. canadensis genotypes. Based on complete and partial sequences of the mitochondrial cox1 gene, interrelationships among E. canadensis genotypes were further examined by both phylogenetic tree and parsimony network analyses using many geographic isolates from various localities. These statistical analyses provided important information to resolve the species status of E. canadensis.
MATERIALS AND METHODS
Taxonomic units and their DNA sequences
Echinococcus taxa used for phylogeny reconstruction were listed in Table 1. Most of their mt genomes have already been published (Le et al. Reference Le, Pearson, Blair, Dai, Zhang and McManus2002; Nakao et al. Reference Nakao, Yokoyama, Sako, Fukunaga and Ito2002, Reference Nakao, McManus, Schantz, Craig and Ito2007). In this study, mt genomes were newly sequenced using eggs of E. felidis isolated from lion faeces in Uganda (Hüttner et al. Reference Hüttner, Nakao, Wassermann, Siefert, Boomker, Dinkel, Sako, Mackenstedt, Romig and Ito2008) and protoscoleces of E. canadensis G10 collected from a lung cyst of a moose in Finland. Additionally, larval isolates of E. granulosus s.s. and E. equinus of previous studies (Nakao et al. Reference Nakao, Li, Han, Ma, Xiao, Qiu, Wang, Yanagida, Mamuti, Wen, Moro, Giraudoux, Craig and Ito2010a; Knapp et al. Reference Knapp, Nakao, Yanagida, Okamoto, Saarma, Lavikainen and Ito2011) were further used for mt genome sequencing. Versteria mustelae (mt genome accession number AB732957) served as an outgroup for inferring the phylogeny of Echinococcus (Nakao et al. Reference Nakao, Lavikainen, Iwaki, Haukisalmi, Konyaev, Oku, Okamoto and Ito2013).
Mitochondrial genomes of Echinococcus spp. analysed in this study
a Numbers having the prefix ‘NC_’, which means complete genomic molecule, are reference sequences (RefSeq) in the GenBank database.
A phylogenetic tree comparison among geographic isolates of E. canadensis genotypes was made by using complete cox1 sequences of the genotypes G6, G7, G8 and G10 from 7 countries. A total of 30 isolates were used for this comparison (Table 2). Among them, the complete sequences were newly determined in 9 isolates of G6 from Russia, Peru and Ethiopia, 2 isolates of G7 from Peru, 3 isolates of G8 from Russia and 12 isolates of G10 from Finland, Russia and USA. The Russian isolates originated from regions of European Russia (Arkhangelsk), Asian Russia (Altai Krai and Altai) and Far East Russia (Sakha). The complete cox1 sequences of E. ortleppi (NC_011122) and E. multilocularis (NC_000928) were added for this comparison, and the latter served as an outgroup.
Isolates used for phylogenetic analyses of E. canadensis genotypes
a The mitochondrial haplotypes were determined based on partial (308 bp) and complete (1608 bp) nucleotide sequences of cox1. n.d., not determined.
b L, larva; A, adult.
Partial cox1 sequences of E. canadensis genotypes were also compared by a haplotype network analysis using 250 geographic isolates from 20 countries (Table 2). For this analysis, 220 cox1 sequences, published in peer-reviewed scientific journals, were selected from DNA databases (Bowles et al. Reference Bowles, Blair and McManus1992; Kamenetzky et al. Reference Kamenetzky, Canova, Guarnera and Rosenzvit2000; Lavikainen et al. Reference Lavikainen, Lehtinen, Meri, Hirvelä-Koski and Meri2003; Obwaller et al. Reference Obwaller, Schneider, Walochnik, Gollackner, Deutz, Janitschke, Aspoek and Auer2004; Bart et al. Reference Bart, Morariu, Knapp, Ilie, Pitulescu, Anghel, Cosoroaba and Piarroux2006; Thompson et al. Reference Thompson, Boxell, Ralston, Constantine, Hobbs, Shury and Olson2006; Varcasia et al. Reference Varcasia, Canu, Kogkos, Pipia, Scala, Garippa and Seimenis2007; Moks et al. Reference Moks, Jõgisalu, Valdmann and Saarma2008; Moro et al. Reference Moro, Nakao, Ito, Schantz, Cavero and Cabrera2009; Sharbatkhori et al. Reference Sharbatkhori, Mirhendi, Jex, Pangasa, Campbell, Kia, Eshraghian, Harandi and Gasser2009; Abushhewa et al. Reference Abushhewa, Abushhiwa, Nolan, Jex, Campbell, Jabbar and Gasser2010; Jabbar et al. Reference Jabbar, Narankhajid, Nolan, Jex, Campbell and Gasser2011; Sharifiyazdi et al. Reference Sharifiyazdi, Oryan, Ahmadnia and Valinezhad2011; Hailemariam et al. Reference Hailemariam, Nakao, Menkir, Lavikainen, Yanagida, Okamoto and Ito2012). The other 30 sequences were our original ones determined in this study or published previously (Nakao et al. Reference Nakao, McManus, Schantz, Craig and Ito2007; Konyaev et al. Reference Konyaev, Yanagida, Ingovatova, Shoikhet, Nakao, Sako, Bondarev and Ito2012). The cox1 sequence of E. ortleppi (NC_011122) was also included in the network analysis.
PCR amplification and sequencing
As reported previously (Lavikainen et al. Reference Lavikainen, Lehtinen, Laaksonen, Agren, Oksanen and Meri2006; Knapp et al. Reference Knapp, Nakao, Yanagida, Okamoto, Saarma, Lavikainen and Ito2011), genomic DNA was extracted from parasites and used as a template for polymerase chain reaction (PCR) to amplify mtDNA fragments. The PCR was performed in a final volume of 50 μl containing 1 μl of genomic DNA (∼100 ng), 0·25 μ m of each primer, 2·5 mm of each dNTP, 0·5 unit of ExTaq DNA polymerase (TaKaRa BIO INC., Otsu, Shiga, Japan) and the manufacturer-supplied reaction buffer. A three-step thermal cycle (94 °C for 30 sec, 55 °C for 30 sec and 72 °C for 1–3 min) was repeated 35 times to amplify mtDNA ranging in length from 0·3 to 3 kb. Direct sequencing of PCR products was carried out using a BigDye terminator cycle sequencing kit and ABI genetic analyzer 3500 (Applied Biosystems, Foster City, CA, USA). Each of PCR primers was used as a sequencing primer. Long PCR products more than 1·5 kb were read by primer walking.
The whole mt genome sequences of E. felidis, E. canadensis G10, E. granulosus s.s. and E. equinus were determined by a method based on sequence islands (Nakao et al. Reference Nakao, Sako and Ito2003). In brief, initial amplifications of mtDNA were carried out using primers designed from conserved regions, and the resultant sequences allowed us to design new primers whereby remaining unknown regions were amplified. All sequences covering complete mtDNA were compiled into a total sequence.
The complete cox1 sequences of geographic isolates of E. canadensis genotypes were determined as reported previously (Hüttner et al. Reference Hüttner, Nakao, Wassermann, Siefert, Boomker, Dinkel, Sako, Mackenstedt, Romig and Ito2008; Konyaev et al. Reference Konyaev, Yanagida, Ingovatova, Shoikhet, Nakao, Sako, Bondarev and Ito2012).
Phylogeny reconstruction by mt genome sequences
Nucleotide sequences of 12 protein-coding genes were translated to amino acid sequences using the echinoderm mitochondrial genetic code (Nakao et al. Reference Nakao, Sako, Yokoyama, Fukunaga and Ito2000; Telford et al. Reference Telford, Herniou, Russell and Littlewood2000). The deduced amino acid sequences were aligned by the T-coffee multiple sequence aligner (Taly et al. Reference Taly, Magis, Bussotti, Chang, Di Tommaso, Erb, Espinosa Carrasco, Kemena and Notredame2011). The software CodonAlign 2.0 (Hall, Reference Hall2004) subsequently created nucleotide alignments under the guidance of the amino acid alignments. Sequences of 2 ribosomal RNA (rRNA) and 22 transfer RNA (tRNA) genes were aligned gene-by-gene using the T-coffee with the optional program R-coffee especially designed for RNA. These alignments were concatenated to make operational data sets as follows: the set #1 of nucleotide sequences for all genes (13 369 sites), the set #2 of nucleotide sequences for 12 protein-coding genes (10 098 sites) and the set #3 of deduced amino acid sequences for 12 proteins (3366 sites).
Maximum likelihood (ML) phylogenies were reconstructed by the program PhyML 3.0 (Guindon et al. Reference Guindon, Dufayard, Lefort, Anisimova, Hordijk and Gascuel2010) using each of the data sets with its appropriate substitution model. The model GTR+G+I was applied for the sets #1 and #2 and the model JTT+G+I for the set #3. The initial phylogenetic tree was estimated by the BIONJ algorithm (Gascuel, Reference Gascuel1997), and the robustness of ML trees were tested by bootstrapping with 500 replicates. Phylogenies were also inferred by the Bayesian program MrBayes 3.2.1 (Ronquist et al. Reference Ronquist, Teslenko, van der Mark, Ayres, Darling, Höhna, Larget, Liu, Suchard and Huelsenbeck2012) using the same data sets and substitution models as those of ML analyses. The Markov chain Monte Carlo (MCMC) analysis was run for 1 million generations and sampled every 100 generations to estimate the posterior probabilities of phylogenetic trees. The run produced 10 000 trees, of which initial 1000 trees were treated as burn-in. The ML and Bayesian inferences were conducted at least twice for each data set to verify the consistency of results. The graphical viewer FigTree 1.3 (http://tree.bio.ed.ac.uk) was used to draw phylogenetic trees. All the trees were rooted with the outgroup taxon V. mustelae.
Pairwise divergence values between E. canadensis genotypes
Using mt genome data of E. canadensis genotypes, the software MEGA 5 (Tamura et al. Reference Tamura, Peterson, Peterson, Stecher, Nei and Kumar2011) was employed to calculate pairwise divergence values of cytochrome b (cob) and cox1 genes. The values of the nucleotide sequences were calculated under the Kimura 2-parameter model (Kimura, Reference Kimura1980) with a gamma setting of 0·5. The values of the deduced amino acid sequences were computed using JTT model with a gamma setting of 0·5. As a genetic yardstick approach, pairwise divergence values of cob and cox1 were computed between intraspecific variants of E. multilocularis (Nakao et al. Reference Nakao, Xiao, Okamoto, Yanagida, Sako and Ito2009) and between closely related sister species of Taenia saginata and T. asiatica (Jeon et al. Reference Jeon, Lee, Kim, Hwang and Eom2005, Reference Jeon, Kim and Eom2007). Furthermore, pairwise comparisons were made between E. ortleppi and E. canadensis genotypes.
Phylogenetic tree and parsimony network analyses
Complete or partial nucleotide sequences of cox1 gene from geographic samples of E. canadensis genotypes were aligned manually. The software TCS 1.2 (Clement et al. Reference Clement, Posada and Crandall2000) was employed to identify individual haplotypes from the alignments (Table 2). Phylogenetic trees of E. canadensis genotypes were inferred by ML and Bayesian analyses using the nucleotide data set of complete cox1 sequences (1608 sites). The parameter settings of the analyses were the same as described above, and the resultant trees were rooted by the outgroup taxon E. multilocularis. A parsimony network of mtDNA haplotypes among various geographic samples of E. canadensis genotypes was illustrated by the TCS 1.2 using the nucleotide data set of partial cox1 sequences (308 sites). The network estimation was run at 95% connection limit. Haplotype frequencies in the network were not used for this study.
RESULTS
Mt genome-based phylogeny of Echinococcus
Phylogenetic trees of Echinococcus were inferred by ML and Bayesian methods using three data sets of mitochondrial genomes. Both statistical methods led to the same result in every data set. The data set #1 of nucleotide sequences for all genes produced a robust phylogeny (Fig. 2A), which was almost identical to a previously published one (Nakao et al. Reference Nakao, McManus, Schantz, Craig and Ito2007). This newly revised tree depicted sister species relationships between E. felidis and E. granulosus and between E. multilocularis and E. shiquicus. Concerning E. canadensis genotypes, G8 and G10 were paraphyletic and the latter was sister to G6 and G7. The same tree topology was obtained from the data set #2 of protein-coding nucleotide sequences (Fig. 2B) and from the data set #3 of deduced amino acid sequences (Fig. 2C). However, branch lengths and node reliabilities of the resultant three trees were different depending on the data sets used.
ML and Bayesian phylogentic trees of Echinococcus inferred from mitochondrial genome data sets. Phylograms were depicted only by ML method, because the resultant topology was the same as that of Bayesian method. Scale bars represent the estimated number of substitutions per site. An outgroup was omitted from the trees. Values of each node are ML bootstrap percentages. Bayesian posterior probability percentages are shown in parenthesis. (A) The phylogram made from the data set of all genes (13 369 nucleotide sites). (B) The phylogram made from the data set of 12 protein-coding genes (10 098 nucleotide sites). (C) The phylogram made from the data set of 12 proteins (3366 amino acid sites).
The mt genome sequences of E. equinus and E. granulosus s.s. determined in this study appeared to contain few intraspecific variations when compared with previously published sequences of other isolates (Le et al. 2002). The present mitochondrial trees including these additional taxa demonstrated a level of sequence variation between E. canadensis G6 and G7 to be almost equivalent to that between the isolates of E. equinus or between the isolates of E. granulosus s.s. (Fig. 2).
A genetic yardstick approach for E. canadensis
By using the complete sequences of cob and cox1 genes, pairwise divergence values were compared in E. canadensis genotypes (Table 3). Characteristics of the nucleotide and amino acid values are quite similar irrespective of the genes examined. The values of concatenated nucleotide sequences (cob+cox1, 2676 sites) among G6, G7 and G10 ranged narrowly from 0·004 to 0·022, showing that G6 and G7 are closely related to G10. The comparisons revealed that G8 is distinct from the other genotypes. In the four genotypes, a maximum divergence value of the concatenated nucleotide sequences reached to 0·04 between G7 and G8.
Pairwise divergence values of mitochondrial cob and cox1 genes between E. canadensis genotypes and between closely related taeniid taxa
a A reference mtDNA sequence of Echinococcus ortleppi (NC_011122) was compared with E. canadensis genotypes.
b A reference mtDNA sequence of E. multilocularis (NC_000928) was compared with that of the most distantly related variant (the mitochondrial haplotype O1) in the Inner Mongolia (Nakao et al. Reference Nakao, Xiao, Okamoto, Yanagida, Sako and Ito2009).
c Reference mtDNA sequences of T. saginata (NC_009938) and T. asiatica (NC_004826) were compared with each other.
As shown in Table 3, the divergence values of the nucleotide and amino acid sequences among G6, G7 and G10 ranged at intraspecific levels when compared with those between geographic variants of E. multilocularis. However, the values between G8 and the other genotypes are slightly lower than those between closely related sister species of T. saginata and T. asiatica, suggesting that G8 is a borderline case in differentiating into species. On the other hand, the values between E. ortleppi and E. canadensis genotypes are higher than or equal to those between T. saginata and T. asiatica, supporting the validity of E. ortleppi.
Relationships among geographic isolates of E. canadensis
To confirm the paraphyly of G8 and G10, ML and Bayesian phylogenetic trees were inferred using complete cox1 sequences of E. canadensis genotypes from various localities (Table 2 and Fig. 3). Both statistical methods resulted in the same tree topology, showing that haplotypes of G10 from Finland, Russia and USA are sister to those of G6 and G7 from Poland, Peru, Kazakhstan, Russia and Ethiopia. Haplotypes of G8 from Russia and USA are sister to a clade including all haplotypes of G6, G7 and G10. Although the number of samples examined were insufficient for firm conclusions, haplotypes of G10 in the Palearctic and Nearctic regions were divided into two groups depending on their locality. A similar tendency was observed in haplotypes of G8.
Phylogenetic relationships among geographic isolates of E. canadensis genotypes inferred by ML and Bayesian methods using the data set of complete mitochondrial cox1 sequences (1608 nucleotide sites). Phylograms were depicted only by ML method, because the resultant topology was the same as that of Bayesian method. A scale bar represents the estimated number of substitutions per site. An outgroup was omitted from the trees.
A haplotype network of E. canadensis genotypes was also made using partial cox1 sequences from worldwide localities (Table 2 and Fig. 4). The sequences from 20 countries were assembled into 17 haplotypes. In this analysis, G6 and G7 were indistinguishable because of low-resolution due to a very short sequence of the gene examined. Seven haplotypes of G6/G7 clustered in conjunction with six haplotypes of G10. The network showed that three haplotypes of G8 were distantly related to the other haplotypes. In particular, haplotypes of G6/G7 did not have any obvious correlation with locality, suggesting a possibility of anthropogenic movement of domestic animals worldwide.
Mitochondrial haplotype network of E. canadensis genotypes. The haplotypes of the genotypes G6/7 (6/7a–6/7 g), G8 (8a – 8c) and G10 (10a – 10f) were based on the short sequences of cox1 obtained from this study and published papers (Table 2).
DISCUSSION
A recent taxonomic revision defined that the family Taeniidae consists of the genera Hydatigera, Taenia, Versteria and Echinococcus, and showed that Versteria is sister to Echinococcus (Nakao et al. Reference Nakao, Lavikainen, Iwaki, Haukisalmi, Konyaev, Oku, Okamoto and Ito2013). A phylogenetic study, based on nuclear DNA sequences, demonstrated the close genetic relationships among members of the genus Echinococcus, indicating that the genus is a young group in which speciation and global radiation occurred rapidly (Knapp et al. Reference Knapp, Nakao, Yanagida, Okamoto, Saarma, Lavikainen and Ito2011). In this nuclear DNA phylogeny, a Neotropical clade consisting of E. oligarthrus and E. vogeli is sister to all other remaining species. The basal position of the Neotropical species is also observed in a phylogeny based on mitochondrial genomes (Nakao et al. Reference Nakao, McManus, Schantz, Craig and Ito2007). A mitochondrial phylogeny revised in this study strongly supported the previous one. Notable points of this revision are the confirmation of sister species relationship between E. granulosus and E. felidis and the clarification of relationships among E. canadensis genotypes. Although the relatedness of the genotypes has remained obscure in the nuclear DNA phylogeny (Knapp et al. Reference Knapp, Nakao, Yanagida, Okamoto, Saarma, Lavikainen and Ito2011), the present study could characterize it because of the rapid evolution of mtDNA as compared with nuclear DNA.
Previous studies using partial sequences of mtDNA (Hüttner et al. Reference Hüttner, Nakao, Wassermann, Siefert, Boomker, Dinkel, Sako, Mackenstedt, Romig and Ito2008) and nuclear DNA (Knapp et al. Reference Knapp, Nakao, Yanagida, Okamoto, Saarma, Lavikainen and Ito2011) suggested that E. felidis is sister to E. granulosus. A mt genome-based phylogeny shown in this study firmly confirmed this sister relationship. Discovery of sister species is important to help understand their evolutionary process. Both E. felidis and E. granulosus are sympatric in Africa, and they are assumed to have a common Asian ancestor (Hüttner et al. Reference Hüttner, Nakao, Wassermann, Siefert, Boomker, Dinkel, Sako, Mackenstedt, Romig and Ito2008). Recent phylogeographic studies suggesting the Middle East origin of E. granulosus (Casulli et al. Reference Casulli, Interisano, Sreter, Chitimia, Kirkova, La Rosa and Pozio2012; Yanagida et al. Reference Yanagida, Mohammadzadeh, Kamhawi, Nakao, Sadjjadi, Hijjawi, Abdel-Hafez, Sako, Okamoto and Ito2012) support this hypothesis. The lion tapeworm, E. felidis, probably invaded Africa together with lions in the late Pliocene, considering the Asian origin of modern Felidae (Johnson et al. Reference Johnson, Eizirik, Pecon-Slattery, Murphy, Antunes, Teeling and O'Brien2006). The present mitochondrial phylogeny also confirms a sister species relationship between E. shiquicus and E. multilocularis.
The utilization of mammalian hosts by Echinococcus species is probably specialized through their evolutionary process in which host-switches occurred. Each of the species depends upon characteristic definitive and intermediate hosts to perpetuate its own life cycle, except for E. canadensis genotypes. The present study strongly suggests that a clade consisting of E. ortleppi and E. canadensis genotypes is the youngest within the genus. In particular, the short branches of E. canadensis genotypes indicate two possibilities: (1) that they are a single species but in the course of speciation, or (2) that they are recently diverged species. As far as we know, T. saginata and T. asiatica could be treated as the closest sister species pair within taeniid parasites. Recent genetic studies revealed that reproductive isolation or hybridization barriers are incomplete between T. saginata and T. asiatica (Okamoto et al. Reference Okamoto, Nakao, Blair, Anantaphruti, Waikagul and Ito2010; Yamane et al. Reference Yamane, Suzuki, Tachi, Li, Chen, Nakao, Nkouawa, Yanagida, Sako, Ito, Sato and Okamoto2012). The maximum divergence values of mitochondrial cob and cox1 genes between E. canadensis genotypes are at slightly lower level than those between T. saginata and T. asiatica, suggesting that all of the genotypes belong to a single species, or that G8 could be separated from G6, G7 and G10 at species level. However, such a genetic yardstick approach can not serve as a universal proxy for species delimitation.
The phylogenetic or cladistic consideration of E. canadensis genotypes is another way to evaluate their species status. As suggested previously (Lavikainen et al. Reference Lavikainen, Lehtinen, Laaksonen, Agren, Oksanen and Meri2006; Moks et al. Reference Moks, Jõgisalu, Valdmann and Saarma2008), the present mt genome-based phylogeny clearly showed that G8 and G10 belonging to the cervid strain are paraphyletic. The tree topology strongly suggests that G6 and G7 (the camel and pig strains) are descendants of a common ancestor, which utilized cervids as intermediate hosts. As summarized in Fig. 1, the close genetic relatedness of G10 to G6 and G7 confirms an invalidity of treating the cervid strain (G8 and G10) as an independent species from the camel and pig strains (G6 and G7). Similarly, the treatment of G6 and G7 as an independent species seems to be inadequate. If a species were to be proposed for G6 and G7, additional species should be prepared for each of G8 and G10. Such a meticulous revision is unrealistic among the closely related genotypes. The strain differentiation of E. canadensis genotypes, dependent on the susceptibilities to intermediate hosts, is very important for epidemiological purposes. Nevertheless, a taxonomic revision based on the strains is unacceptable in view of molecular cladistics. By disregarding a concept of the strains, G6, G7 and G10 might be assembled as a separate species from another species for G8. Ecological, morphological or genetic criteria to differentiate G6, G7 and G10 from G8 are prerequisites for assigning species rank, but presently there are no effective criteria. Since G8 and G10 are distributed sympatrically (Thompson et al. Reference Thompson, Boxell, Ralston, Constantine, Hobbs, Shury and Olson2006; Moks et al. Reference Moks, Jõgisalu, Valdmann and Saarma2008; Himsworth et al. Reference Himsworth, Jenkins, Hill, Nsungu, Ndao, Thompson, Covacin, Ash, Wagner, McConnell, Leighton and Skinner2010; Bryan et al. Reference Bryan, Darimont, Hill, Paquet, Thompson, Wagner and Smits2012), splitting into two subspecies is also inadequate for them.
Various geographic isolates of E. canadensis genotypes were also analysed in this study in order to clarify their interrelationships. The paraphyly of G8 and G10 was highly supported by a phylogenetic tree analysis using complete cox1 sequences. Moreover, a haplotype network analysis based on the short cox1 sequences of many geographic isolates confirmed the close relatedness of G10 to G6 and G7. Haplotypes defined in this study did not have strict correlations with geographic distribution, particularly in haplotypes belonging to G6 and G7. However, an assemblage of interesting haplotypes was found in the Far East Russia. A Russian isolate derived from a reindeer was identified as the G6 haplotype EcRUS1 (Table 2 and Fig. 3). As far as we know, this is the first case of G6 infection in cervid, implying that G6 evolved from an ancestor susceptible to cervid (Fig. 1). More interestingly, this haplotype was found in the Sakha (Yakutia) Republic, together with the G8 haplotype EcRUS3 (collected from moose and wolf) and the G10 haplotypes EcRUS4 (moose), EcRUS5 (wolf), EcRUS6 (reindeer) and EcRUS7 (human). The sympatric distribution of G6, G8 and G10 in Yakutia is quite important in studying their species status and evolutionary process in wildlife.
The specific name of E. canadensis seems to be the most suitable for handling the closely related genotypes. However, in this article we withhold our judgement on their species status. The final decision will be made after completion of population genetic studies using nuclear DNA markers (Knapp et al. Reference Knapp, Nakao, Yanagida, Okamoto, Saarma, Lavikainen and Ito2011).
ACKNOWLEDGEMENTS
We are grateful to K. Beckmen and A. Oksanen for providing metacestodes of E. canadensis and to M. Wassermann and T. Romig for providing eggs of E. felidis. Thanks are also due to A. Saveljev for his help in sample collection.
FINANCIAL SUPPORT
This study was supported by a Grant-in-Aid for Scientific Research (No. 22590376) from Japan Society for the Promotion of Science (JSPS) to M. N., by RFBR research project (Nos. 13-04-92107 and 13-04-10140) to S. K. and by Grants-in-Aid for Scientific Research (Nos. 21256003 and 24256002) from JSPS, JSPS-Asia/Africa Scientific Platform Fund (2006–2011) and the Special Coordination Fund for Promoting Science and Technology from the Ministry of Education, Culture, Sports, Science and Technology in Japan (MEXT) (2010–2012) to A. I.
