To the best of our knowledge this is the first application of the NGS approach to completely characterise mt genomes and nuclear rRNA operons of a trematode in a single step. Given the fact that within parasitic flatworms (Neodermata) there have been just two monogenean species of Gyrodactylus sequenced using NGS platform (Illumina HiSeq 2000^TM) to develop an automated in silico mt genome assembly approach [18], this is the first use of a NGS platform with limited output to easily characterise multiple-copy genomic loci within parasitic flatworms. We combined a total of four indexed TruSeq Nano library samples (only two samples published here) prepared from untreated DNA extractions on a single Illumina MiSeq run to assess the capability of a simple and straightforward approach to reliably characterise the genomic loci of interest. Of the two loci, nuclear rRNA operons received about five- to nine-fold greater coverage than the mt genomes, whose minimal values in turn were 61 and 53 reads in D. spathaceum and D. pseudospathaceum, respectively. We consider this a sufficiently high coverage that allows for combining two- or three-fold more flatworm specimens on a single MiSeq run using the same library preparation strategy as employed here in order to ensure that characterisation of complete novel mt genomes and nuclear rRNA operons is as cost-effective as possible while receiving sufficient coverage in future studies.

Our study made available reference mt genomes of closely related Diplostomum spp. thus providing a rich resource for future approaches to species delimitation that, in turn, will enable the exploration of molecular epidemiology within a large group of widespread fish pathogens. Recent application of cox1 ‘barcoding’ proved essential in discovering vast previously unrecognised species diversity within the genus [8, 13, 14]. The currently used primers for amplification of the cox1 ‘barcode’ region of [12] work reasonably well and seem to allow for differentiation and identification of a number of Diplostomum spp. However, problematic low divergence has already been demonstrated within three cryptic species complexes, “D. baeri”, “D. huronense” and “D. mergi” (see [8, 14]). We predict that expansion of the exploration of genetic diversity across host populations will inevitably reveal cryptic species thus making species delineation crucial for defining the host ranges and geographical distribution of Diplostomum spp. and the successful development of population genetic studies. Our study profited from the availability of DNA from the adult stages of Diplostomum spp. and taxonomic expertise. Given that most of the recently molecularly delimited species-level lineages of Diplostomum cannot be formally described (awaiting the collection of the adult stages), we see an increasing need to employ the NGS approach presented here to characterise any morphologically determinable Diplostomum specimens available to allow for building a solid baseline for molecular barcoding on which ecological and epidemiological studies could be based in future. Moreover, it is already exciting to visualise an idea of using the outputs of NGS, most easily the mt genome sequences, to test their power in addressing the phylogenetic interrelationships of Diplostomum spp. when analysed at the amino acid level.

Although the currently utilised cox1 ‘barcode’ fragment seems useful in assigning individual isolates of Diplostomum to both described species and novel molecularly defined lineages, its use for species delineation has a number of disadvantages. First, single molecular markers for robust molecular systematic estimates are insufficient, and multi-gene approaches are preferred; the recently applied coalescent-based approaches [14] and the molecular resources provided here hold significant promise for species delimitation within Diplostomum. Secondly, phylogenetic estimates based on the cox1 ‘barcode’ fragment [8, 14] generally lack nodal support at the internal nodes that define the interrelationships of individual lineages and therefore the phylogenetic utility of this region is limited. Finally, and most importantly, we have shown that the cox1 gene in fact represents the most conserved protein-coding region of the mt genome of Diplostomum spp. (Fig. 1), as already shown for other parasitic flatworms (e.g. [35, 36]) and that the ‘barcode’ region is even more conserved. The comparative mt genome assessment allowed us to identify genes with the greatest interspecific variation that are promising for the development of molecular diagnostic markers for species recognition [35]. Certainly the best candidates are the nad4 and nad5 genes that are both relatively long and the least conserved, even when comparing such closely related taxa as D. spathaceum and D. pseudospathaceum, and thus have the potential to add phylogenetic signal and possibly resolve current taxonomic problems. The use of these new markers in large-scale population studies on Diplostomum spp. will offer the advantage of simultaneous species delimitation and assessment of intraspecific genetic variation and thus boost molecular epidemiology studies based on the most accessible life-cycle stages, i.e. the larval forms in the snail and fish populations.

The accelerated development of methods for next-generation biodiversity assessment such as environmental DNA (eDNA) and metabarcoding [37] offer significant promise for large-scale spatial studies related to disease ecology. The first study using water samples to assess the presence of the pathogenic parasite of frogs Ribeiroia ondatrae in wetland habitats indicates the high potential for detecting free-living larval stages of macroparasite infectious agents in the aquatic environment [38]. Our study provides a resource for the design of primer pairs targeting very short DNA sequences that would allow the use of NGS tools to identify unambiguously Diplostomum spp. in both bulk samples (natural assemblages of metacercariae in fish eyes) and freshwater samples containing DNA of the dispersive free-living larval stages (miracidia, cercariae).

The novel data on mt genomes and nuclear rRNA operons of Diplostomum spp. allowed the first genome wide estimation of the phylogenetic relationships of the order Diplostomida, one of the two fundamental lineages from which extant digeneans have diversified [39]. Surprisingly, the currently available mitogenomic data for the Digenea analysed at both the amino acid and nucleotide levels invariably recovered, albeit with much weaker statistical support, the Diplostomidae as a sister lineage of the order Plagiorchiida rather than as a basal lineage of the Diplostomida, the latter inferred by both the benchmark phylogeny of the Digenea of [39] and the current rDNA-based phylogenetic analysis for the set of taxa used in mt genome-based analysis (Fig. 2). This is perhaps the most striking finding given the significance and depth of the basal dichotomy in the Digenea between the Plagiorchiida and Diplostomida on both the molecular [39] and life history levels [40]. If the mt-based phylogeny reflects actual organismal phylogeny, with Diplostomidae as an early diverging lineage within the Plagiorchiida, the seemingly synapomorphic life history characteristics of Diplostomida would need to be reconsidered. For example, the active penetration of a vertebrate host by cercaria in the Diplostomida (with the exception of Brachylaimoidea), or the utilization of tetrapod vertebrates as definitive hosts (with the exception of the Aporocotylidae), would be viewed as plesiomorphies for the Digenea. However, addressing such questions would currently be premature since the mt genomes of Diplostomum spp. could be compared solely with those of taxa representing the most derived lineages of both the Diplostomida and the Plagiorchiida. Additional mt genome data, especially for the earlier diverging lineages of the Diplostomida (Brachylaimoidea, Clinostomidae, Aporocotylidae and Spirorchidae) and Plagiorchiida (e.g. Bivesiculoidea, Transversotrematoidea) would allow testing of this hypothesis.

In light of the conflict between nuclear ribosomal and mitochondrial protein estimates of phylogeny, it may be less surprising that mt genome organisation of Diplostomum spp. also closely matches that of the nearly perfectly conserved mt genomes of representatives of the Plagiorchiida, rather than that of any representative of the Schistosomatidae, the only members of the order Diplostomida with mt genomes characterised prior to our study. Compared with the mt gene order in all species of the Plagiorchiida, there was a single minor alteration in the mt genomes of Diplostomum spp., i.e. a reciprocally switched position of immediately neighbouring genes trnP and trnN (Fig. 3). This tRNA gene position switch likely represents an autapomorphy for Diplostomum, because trnP is always found immediately downstream to trnN throughout remaining trematodes as well as cestodes, the sister lineage of the trematodes (e.g. [34]). Species of the plagiorchiidan genera Dicrocoelium, Paragonimus and Paramphistomum then display a further common switch of neighbouring trnE and trnG gene positions; however, this minor reorganisation does not seem to reflect actual cladogenetic events and might have arisen several times during digenean evolution. Contrary to the members of the Plagiorchiida, species of the Schistosomatidae display a significantly altered gene order, including a few tRNA transpositions (as seen in Trichobilharzia regenti and two Asian schistosomes) and two larger mt genome region rearrangement events within the more derived African schistosomes (see [17, 41] for details). Further mt genome data for a wider range of taxa, especially for representatives of the remaining families of the Diplostomida (see above), would provide important insights into the patterns and possible mechanisms of mt genome evolution in early divergent digeneans.