The genome of E. oligarthrus was assembled from a combination of two paired-end libraries sequenced in the Illumina MiSeq platform. The genome assembly involved several steps. First, a preliminary de novo assembly using SPAdes 3.6 (Bankevich et al., 2012; Safonova et al., 2014) was performed. The assembly sequences were screened against the NCBI nt database, using Nucleotide–Nucleotide BLAST v2.6.0+ (available at ftp://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/nt.gz) in megablast mode, with an e-value cutoff of 1e−25 and a culling limit of 2 and using DIAMOND tblastx against SwissProt (Buchfink et al., 2015). Raw, paired-end Illumina reads were mapped against the assembly using Bowtie2 (Langmead and Salzberg, 2012). The output was converted to a BAM file using Samtools (Li and Durbin, 2009). Blobtools v1.1 (Laetsch and Blaxter, 2017) was used to create taxon-annotated GC-coverage plots for E. oligarthrus genome assembly and to identify the target sequences and target reads using as input the Nucleotide–Nucleotide and DIAMOND BLAST (Buchfink et al., 2015) and the raw read mapping results. Sequences that did not match the Platyhelminthes taxon as a top BLAST hit at the phylum level were filtered out. After removing contaminants, sequences that did match as a top BLAST hit at the phylum level and whose base coverage was >4× were used to recover the target reads. Before the assembly, we used Jellyfish to create a k-mer histogram and the convergence was analyzed with genomescope in order to evaluate whether after cleaning with blobtools the remaining reads could be used to obtain an acceptable assembly. The target reads extracted from the bam files were used to re-assemble the E. oligarthrus genome using SPAdes 3.12 (Bankevich et al., 2012; Safonova et al., 2014). The scaffolds of E. oligarthus were obtained with Chromosomer (Tamazian et al., 2016) and using E. multilocularis as the reference genome [WormBase ParaSite, Version WBPS12 (WS267)]. Redundant unplaced and unlocalized sequences were discarded from the final assembly. The fragment sequence was considered unplaced if two or more alignments located on different reference sequences or unlocalized if two or more alignments located on the same reference sequences according to Tamazian et al. (2016). The contigs that did not map to the reference genome were included in the final assembly. The standard quality metrics of the assembly such as N50, the total number of contigs and the total length of the assembly were evaluated using QUAST (Gurevich et al., 2013). Putative non-target contigs shorter than 500 bp in length were removed from the final assembly. The completeness of the gene space was validated using BUSCO2 (Simão et al., 2015) and eukaryote CEGGs database. Furthermore, the core of cestodes genes [genes contained in all cestode species according to Maldonado et al. (2017)] was screened using BLAST. Coverage and depth coverage were calculated with custom scripts. Depth coverage refers to the number of times that the same region or position in the reference genome is represented by the assembled genome. Coverage refers to the percentage of the total length of the reference genome that is represented by the assembled genome.

The gene annotation of E. oligarthrus was performed by transferring the gene annotations from the E. multilocularis genome [WormBase ParaSite, Version WBPS12 (WS267)] using a manual and own scripting approach. CDS and protein sequences of E. multilocularis and the scaffolds of E. oligarthrus were used for this purpose. First, the gene allocated regions were identified using BLAST with an e-value cutoff of 1e−10 and one best target hit. The CDS fragments of E. oligarthrus were extracted using bedtools (Quinlan and Hall, 2010; Quinlan, 2014) and the transcripts were re-assembled using Chromosomer (Tamazian et al., 2016). GeneWise (Birney et al., 2004) was used to find the correct frameshift of the CDS and to obtain the final datasets of CDS and proteins. The coordinates of the annotations were assessed using Exonerate (Slater et al., 2005) and added to a final GFF file. The performance of gene annotation and basic statistics for E. oligarthrus gene models, including the average intron/exon lengths and the number of introns, were calculated using Eval (Keibler and Brent, 2003). The functional gene annotation was performed using InterProScan-5.7-48.0 (Quevillon et al., 2005) and InterPro2GO databases were used to assign Gene Ontology (GO) terms (Ashburner et al., 2000). Gene models were subjected to BLAST search (Stemmer et al., 2013) against UniprotDB and Blast2GO (https://www.blast2go.com/) was used to define the final annotations and GO terms stats. The GO terms were analyzed using GO.db database implementing R, Bioconductor version: Release (3.8) with custom scripts. Proteins studied in this work were searched using BLAST (Stemmer et al., 2013) against UniProtKB/Swiss-Prot databases. Protein domains were screened against PFAM, and Prosite databases using PFAM_scan (Finn, 2006) or HMMscan 3.0. Global pairwise alignments were performed using Needle software (Rice et al., 2000) in order to identify the orthologous genes related to host–parasite interactions between E. oligarthrus and E. multilocularis (Brehm and Koziol, 2017). Identity and coverage stats were calculated and only hits with more than 50% of coverage were selected.

Molecular markers of the Echinococcus nuclear genome were downloaded from GenBank (Kinkar et al., 2017). For each locus, homologous regions were extracted from complete Echinococcus genomes [WormBase ParaSite, Version WBPS12 (WS267)] and the E. oligarthrus genome obtained in this work using custom scripts. Homologous genomic sequences from Taenia solium were employed as outgroup. The DNA sequences were aligned with ClustalX (v2.0.12) and multiple alignments were edited with BioEdit (v7.1.3). Phylogenetic analyses of concatenated nuclear markers were performed using the maximum likelihood method and Tamura-Nei model implemented by MEGAX software (Kumar et al., 2018). The bootstrap consensus trees inferred from 1000 replicates were retained. The percentage of replicate trees in which the associated taxa clustered together in more than 50% in the bootstrap test are shown for each node. The concatenated data resulted in 1552 bp length for seven taxa. A second phylogeny was calculated using MrBayes 3.1.2 (Ronquist et al., 2012). The evolutionary model was set on generalized time reversible substitution model with gamma-distributed rate variation across the sites and a proportion of invariable sites (GTR +G + I model) with at least 200 samples from the posterior probability distribution, and diagnostics calculated every 1000 generations. We use an exponential prior on the branch length with mean = 0.1 substitutions/site.

Variant calling was performed as described in the “Whole-Genome SNP analysis” section. Genome-wide SNPs were used to perform phylogeny analysis as follows: first, heterozygous SNPs were removed and only the homozygous SNPs with a depth coverage >20× and strictly covered in all the Echinococcus species were selected to correct for complete lineage sorting. Afterwards, the homozygous SNP loci were concatenated and the resulting sequence alignment was used to create a phylogenetic tree by implementing the Maximum Likelihood and Bayesian method as above. PartitionFinder 2 (Lanfear et al., 2016) was used to select the best-fit partitioning schemes and models of evolution for the phylogenetic analysis. The selected model was implemented to perform phylogeny analysis.

Genome sequences of Echinococcus species were ordered in chromosomes according to the last version of the E. multilocularis genome [WormBase ParaSite, Version WBPS12 (WS267)] using Chromosomer (Tamazian et al., 2016) and were used as mapping templates for the variant calling analyses. Variant calling was performed as follows: first, all of the raw reads from E. canadensis G7, E. granulosus G1, E. multilocularis, and E. oligarthrus libraries were processed and filtered by quality. Then, the reads were first mapped against their own reference genomes and then against the genomes of the other three Echinococcus species using bowtie2 (Langmead and Salzberg, 2012). Mapping statistics were calculated with Bamtools (Barnett et al., 2011), and duplications were marked and discarded using picard-tools v-2.18 (http://broadinstitute.github.io/picard/). The variant calling was performed using bcftools (Narasimhan et al., 2016; Danecek and McCarthy, 2017) and GATK (McKenna et al., 2010) using the following parameters: variation frequency was set >40% with a depth coverage of at least 70% of the total mean coverage; the base quality of both reference site and variation site was set to >30. Insertion and deletions (indels) were filtered out with VCFtools (Danecek et al., 2011), and SNPs with less than 10 bp far from indels were removed to avoid false-positive SNP. In order to annotate the heterozygous and homozygous polymorphic sites, the reads were first mapped against their own reference genomes and then against the other target genomes. Heterozygous sites were retained only if both forward and reverse reads mapped against the reference and alternative allele at a given nucleotide position whose depth coverage was at least 70% of the total mean coverage supporting that position. Homozygous polymorphic sites were annotated if the forward and reverse reads mapped onto the alternative allele with at least 70% of the total mean coverage supporting that position and if there were no reads supporting the reference allele. Homozygous and heterozygous variant sites were registered for all of the species. The Transition/transversion ratios were calculated using VCFtools (Danecek et al., 2011), and the annotation and classification of SNPs based on the effect of annotated genes were carried out with SnpEff v4.0 (Cingolani et al., 2012). Graphics were built using R software. (https://www.r-project.org/).

The E. oligarthrus genome was assembled from two paired-end libraries sequenced in Illumina MiSeq. High-quality genomic DNA was purified from two cysts isolated from its natural host, D. azarae, which was naturally infected and found in Iguazú National Park, Misiones province, Argentina. Microscopic and macroscopic findings indicated that it was a neotropical Echinococcus species. The species and the genotype were confirmed by polymerase chain reaction (PCR) amplification of cytochrome oxidase 1 (cox1) followed by direct sequencing. The E. oligarthrus group 1 cox1 sequence was determined (Arrabal et al., 2017). The E. oligarthrus genome assembly was performed using a de novo assembly strategy, and the best assembly was chosen based on its quality metrics. Due to the type of the metacestode development that consists of close interaction with the host’s tissue and in spite of having done a careful extraction of the parasite material, we sequenced a high proportion of host’s DNA. Indeed, preliminary analysis of raw reads revealed that the sequencing yield of parasite DNA was ∼22%, whereas the remaining ∼78% of the sequences was from presumably D. azarae, whose genome has not been sequenced yet. Therefore, we performed a screening of potential contaminants before obtaining the final assembly. Since the genome of D. azarae is unknown and non-related sequences are available, we made a thorough identification of sequences truly derived from the target genome using Blob tools software (Laetsch and Blaxter, 2017). This software was specially designed for DNA contamination assessment and is particularly useful for organisms of parasitic origin. Here, we performed a taxonomic selection of the target sequences using the Platyhelminthes taxon (Taxonomy ID: 6157) and created a taxon-annotated GC-coverage plot that allowed us to identify Platyhelminthes contigs and discard all the non-target sequences ( Supplementary Figure 1 ). After removing the contaminants, sequences that did match as a top BLAST hit at the phylum level and whose base coverage was >4× and %GC ∼ 41% (typical GC content of Echinococcus genus) were used to recover the target reads and re-assemble the E. oligarthrus genome (see Materials and methods for more details). The final assembly contained 74513 contigs that were further used to obtain the genome scaffolds. Redundant unplaced and unlocalized sequences (19.7 Mb) were removed from the contigs to obtain the final scaffolds assembly (see Materials and methods for more details). The contigs that did not map to the reference genome (did not allocate in scaffolds) were included in the final assembly. The final E. oligarthrus genome assembly was composed of 3764 sequences comprising a total of 86.22 Mb with an average GC content of 41% and showed an N50 ∼ 10 Mb and ∼20× in depth coverage. The assembly metrics for contigs and scaffolds are shown in Table 1 and Supplementary Table S1 . The nuclear genome of E. oligarthrus showed 76.3% coverage and ∼90% identity with the genome of E. multilocularis whose genome size is ∼115 Mb [WormBase ParaSite, Version WBPS12 (WS267)]. The percentage average of nucleotide identity among Echinococcus chromosomes ranged from 88.7% to 92.7% ( Supplementary Table S2 ). The assembly also included the E. oligarthrus mitochondrial genome composed of one scaffold that was obtained from the contigs Eoli_00665 and Eoli_02629. The scaffold length was 13,893 bp with 98% coverage and 96% nucleotide identity to the E. oligarthrus mitochondrial reference genome (GenBank accession number AB208545). As expected, the nucleotide identity of the mitochondrial genomes was lower in comparison with other Echinococcus species ( Supplementary Figure 2 ). To assess the completeness of the genome assembly, we evaluated the gene space using “Benchmarking Universal Single-Copy Orthologs” (BUSCO) (Simão et al., 2015), which measures the genome completeness based on evolutionarily informed expectations of gene content. In this analysis, we identified 77.6% (235 of the 303 core genes) that are expected to be present in all metazoans, including 125 complete and duplicated, 110 fragmented, and 68 missing orthologs. The fragmented nature of the assembly may have prevented many genes from meeting the stringent matching criteria implemented by BUSCO. Indeed, the BLAST results suggest that most of the core genes are identifiable in the genome, even though many genes are present as fragments within the assembly. Also, due to the fragmentation level of the genome assembly, a hybrid and ab initio gene prediction is likely impracticable. Since all the available gene annotation transfer tools tried here showed low efficiency, the gene annotation was performed by manual gene annotation transfer as described in the Materials and methods section. First, we used BLAST to search for the best suitable Echinococcus species to be used as a template. In this regard, we found that E. multilocularis was the most suitable species because of the integrity of the assembly and the higher identity with the genome of E. oligarthrus ( Supplementary Table S2 ). The final set of E. oligarthrus genes comprised 8753 genes coding for proteins (4494 genes with coverage > 50% and 4259 genes with coverage < 50%) ( Supplementary Table S3 ). The whole set of genes represented the 85.0% of the total genes expected to be found in species of the genus Echinococcus. The predicted proteins were also screened against the conserved core of genes that are present in all the cestodes species according to Maldonado et al. (2017) using BLAST. In this regard, 4872 out of 5203 genes were found for E. oligarthrus using e values of <1e−12, which comprised the 93.6% of the total conserved core of genes in cestodes, indicating that the genome contains useful molecular data ( Supplementary Table S4 ).

GO terms were assigned to the 40% of E. oligarthrus proteins. In relation to the Molecular Function GO terms frequency, the two main categories found were “binding (GO:0005488)” and “catalytic activity (GO:0003824)”, which is in accordance with the GO terms frequency observed in other cestode genomes (Hahn et al., 2014) and Echinococcus species (Maldonado et al., 2017). For the Biological Process GO term, the highest frequencies observed were the categories “cellular process (GO:0009987)” and “metabolic process (GO:0008152)”, also in accordance with the GO terms frequency observed in the related organism ( Supplementary Figure 3 ).

From the total gene-set, we searched for genes that have already been described as having a role in host–parasite interactions in Echinococcus (Brehm and Koziol, 2017). In this regard, we found 56 genes involved in host–parasite interactions whose coverages and identities to E. multilocularis orthologous genes were >50% ( Supplementary Table S5 ). Particularly, two genes encoding for the epidermal growth factor (EGF) tyrosine kinase receptors (Eoli_000075800 and Eoli_000617300) with a high percentage of identity (59.8% and 52.6%) and coverage (65.3% and 51.4%) were identified. Also, we found a gene encoding for fibroblast growth factor (FGF) receptor tyrosine kinase (Eoli_000833200) whose percentage of identity and coverage values were 76% and 100%, respectively. Regarding nuclear receptor hormones, we found nine orthologous genes with high identity and coverage (83.8% and 51.8% on average, respectively) including the cestode-specific nuclear hormone receptor Eoli_000937000. Non-kinase receptors were also found, comprising four genes, including Frizzled G protein-coupled receptors (GPCR) (Eoli_000682100) with high identity and coverage (88.6% and 87.6%, respectively). Finally, we identified the amino acid transporters (DAACS family), lipid binding proteins (FABPs), and antigens (AgB and Eg95), specific and conserved proteins in cestodes/Echinococcus.

The nucleotide variation in the genomes of different Echinococcus species was assessed through variant calling analysis. Genetic variants were identified among the genomes of E. oligarthrus, E. multilocularis, E. granulosus G1, and E. canadensis G7. Here, we focused on the study of single-nucleotide polymorphisms (SNPs). In order to perform this analysis, NGS raw reads were mapped against a reference genome composed of chromosomes and contigs. For all the analyses, the reads were first mapped against their own reference genomes and then against the genome of the corresponding analyzed species. Homozygous and heterozygous variant sites were identified and were marked in both the reference and the alternative allele (see Materials and Methods for more details). First, we evaluated the intraspecific variation of the Echinococcus genomes. As described in Maldonado et al. (2017). E. granulosus G1 genome exhibited the highest number of intraspecific variant sites (74,796 SNPs, 0.65 SNPs/kb), followed by E. oligarthrus (23,301 SNP, 0.23 SNP/kb), E. canadensis G7 (10,791 SNPs, 0.095 SNPs/kb), and E. multilocularis (1,287 SNPs, 0.011 SNPs/kb). With regard to the genetic diversity observed among the Echinococcus species, we observed that the SNP distribution was similar to the distribution described in our previous research (Maldonado et al., 2017). E. canadensis G7 and E. granulosus G1 showed a higher number of SNPs (842,322, Ts/Tv = 2.97) between each other than between E. canadensis G7 and E. multilocularis (314,176, Ts/Tv = 2.97), and in comparison to E. granulosus G1 and E. multilocularis (272,138, Ts/Tv = 2.94). Furthermore, the pair E. canadensis G7 and E. granulosus G1 also showed almost 10 orders of magnitude more SNPs than between E. oligarthrus and E. multilocularis (38,911, Ts/Tv = 3.25). The number of SNPs between E. oligarthrus and E. granulosus (126,472 Ts/Tv = 2.99) was similar to the number of SNPs between E. oligarthrus and E. canadensis G7 (171,135, Ts/Tv = 3.00). We also identified homozygous and heterozygous variant sites for the four Echinococcus species, for both the reference and the alternative allele in each case. The number of homozygous SNPs was 313,992 for E. canadensis G7 and E. multilocularis, 266,180 for E. granulosus G1 and E. multilocularis, 38,762 for E. oligarthrus and E. multilocularis, 830,768 for E. canadensis G7 and E. granulosus G1, 125,147 for E. oligarthrus and E. granulosus G1, and 170,049 for E. oligarthrus and E. canadensis G7 ( Table 2 ). Moreover, the SNP density was assessed in terms of the number of SNPs per 1-Mb length of the Echinococcus chromosomes. For all the Echinococcus genomes, the highest SNP density (SNPs/Mb of the chromosome) was found for chromosome 1. The chromosomes 2, 3, and 4 showed a slightly lower SNP density than chromosome 1 and chromosomes 5, 6, 7, 8, and 9 showed the lowest SNP density ( Figure 1 ). However, the intraspecific SNP density distribution for E. granulosus G1 was higher in chromosomes 1, 5, and 9. We also evaluated the SNP distribution of E. oligarthrus in the coding and non-coding genomic regions (intergenic, exons, and introns) of E. multilocularis. In this regard, the 38,911 SNPs distributed with a higher rate in exons and introns rather than in the intergenic regions. However, those distributed in the coding regions exhibited a higher rate of synonymous changes (67.3%) rather than missense changes (32.6%) ( Supplementary Figure 4 ).

In order to evaluate the contribution of the SNPs to the genetic diversity among the different Echinococcus species, we performed phylogenetic analyses by implementing three different approaches: the first used the whole genome variant sites for the four Echinococcus species (only the SNPs), which consisted of analyzing 244,246 sites in each genome arising a total of 2.08 Mb; the second approach used genome regions for the four Echinococcus species whose depth coverage was >20× and that strictly contained variant sites supported by more than 20 reads, which involved the analysis of 40,179,279 sites in each genome arising in a total of ∼162 Mb. For the third approach, we used only coding regions for the four Echinococcus species whose depth coverage was >20× and that strictly contained variant sites supported by more than 20 reads. This involved the analysis of 42,200 sites in each genome arising in a total of 168.8 kb. In all the cases, only the homozygous SNPs were used to perform phylogenetic analyses. The sites and sequences retained here were concatenated and the resulting alignment was used to create the phylogenetic trees by implementing the maximum likelihood and Bayesian methods (see Materials and methods for more details).

The construction of the phylogenetic tree implementing the Bayesian method and using whole-genome SNPs showed a topology that demonstrated a higher genetic distance between E. canadensis G7 and E. granulosus G1 in comparison with the common node from which E. multilocularis and E. oligarthrus diverge equidistantly at a very short distance ( Figure 2A ). E. multilocularis and E. oligarthrus exhibited a low genetic diversity between each other. PartitionFinder 2 (Lanfear et al., 2016) was used to select the best-fit partitioning schemes and models of evolution for the phylogenetic analysis. Transversional substitution model (TVMef) was the best-fitted model. Phylogenetic analysis using TVMef was consistent with the previous result ( Figure 2B ). The topologies of the trees were also consistent with the topology observed using the Maximum Likelihood method ( Supplementary Figure 5A ). Phylogenetic analysis using genomic regions with depth coverage >20× sites and coding regions with depth coverage >20× also showed the same topology, except for the branches that exhibited different lengths ( Supplementary Figures 5B, C ). The total length of the coding sequences used here was ∼42.2 kb and the genes sampled and used for this purpose are listed in Supplementary Table S6 .

In addition, we registered how many polymorphic loci are shared among the Echinococcus species containing the same polymorphism. In previous research, we determined that the number of shared loci was higher when we used E. multilocularis as the reference genome, demonstrating to be the basal species and discarding E.granulosus G1 and E. canadensis G7 as possible candidates (Maldonado et al., 2017). In order to evaluate whether E. oligarthrus occupies a basal position in the genus, we incorporated E. oligarthrus and analyzed the number of shared loci using both E. multilocularis and E. oligarthrus as reference genomes and compared the results to each other. Due to the different genome quality assemblies and coverage, we normalized the number of SNPs in shared loci using the effective length of the sampled regions under the assumption that all the genome regions are equally subjected to mutation. In this regard, we found that 12,282 loci shared the same nucleotide change among E. canadensis G7, E. granulosus G1, and E. oligarthrus with respect to E. multilocularis (E. multilocularis used as reference). Moreover, the number of shared loci between E. oligarthrus and E. granulosus G1, between E. oligarthrus and E. canadensis G7, and between E. canadensis G7 and E. granulosus G1 was 4115, 8919, and 92,326, respectively. Most of the loci were unique for E. granulosus G1 (226,350) and E. canadensis G7 (335,728) rather than for E. oligarthrus (40,386). On the other hand, the number of shared loci containing the same nucleotide change among E. canadensis G7, E. granulosus G1, and E. multilocularis was 12,277 with respect to E. oligarthrus (E. oligarthrus used as reference). The number of shared loci between E. multilocularis and E. granulosus G1, between E. multilocularis and E. canadensis G7, and between E. canadensis G7 and E. granulosus G1 was 8675, 8875, and 54,792 respectively. Furthermore, and similar to the above results, the highest numbers of unique loci were for E. granulosus G1 (133,417) and E. canadensis G7 (242,815), rather than for E. multilocularis (60,455) ( Figure 3 ).

In order to gain accuracy and supporting evidence for our previous results, we reconstructed the Echinococcus phylogeny with nuclear molecular markers previously described (Kinkar et al., 2017). Nuclear molecular marker sequences of all the available Echinococcus species were concatenated and then were aligned and analyzed using both the Bayesian and maximum likelihood methods. This analysis involved the study of 1552 nucleotides that allowed to root the trees and reinforced the basal position of E. oligarthrus ( Supplementary Figure 5D ). The phylogenetic topology obtained was consistent with previous phylogenetic analyses of (Nakao et al., 2013b) that have placed E. oligarthrus (and E. vogeli) in a basal position for the genus Echinococcus. Similar phylogenetic topologies were obtained with the two methods employed here ( Supplementary Figure 5E ).