Skip to main content

Sequencing of the complete mitochondrial genome of a fish-parasitic flatworm Paratetraonchoides inermis (Platyhelminthes: Monogenea): tRNA gene arrangement reshuffling and implications for phylogeny



Paratetraonchoides inermis (Monogenea: Tetraonchoididae) is a flatworm parasitising the gills of uranoscopid fishes. Its morphological characteristics are ambiguous, and molecular data have never been used to study its phylogenetic relationships, which makes its taxonomic classification controversial. Also, several decades of unsuccessful attempts to resolve the relationships within the Monogenea present a strong indication that morphological datasets may not be robust enough to be used to infer evolutionary histories. As the use of molecular data is currently severely limited by their scarcity, we have sequenced and characterized the complete mitochondrial (mt) genome of P. inermis. To investigate its phylogenetic position, we performed phylogenetic analyses using Bayesian inference and maximum likelihood approaches using concatenated amino acid sequences of all 12 protein-coding genes on a dataset containing all available monogenean mt genomes.


The circular mt genome of P. inermis (14,654 bp) contains the standard 36 genes: 22 tRNAs, two rRNAs, 12 protein-encoding genes (PCGs; Atp8 is missing) and a major non-coding region (mNCR). All genes are transcribed from the same strand. The A + T content of the whole genome (82.6%), as well as its elements, is the highest reported among the monogeneans thus far. Three tRNA-like cloverleaf structures were found in mNCR. Several results of the phylogenomic analysis are in disagreement with previously proposed relationships: instead of being closely related to the Gyrodactylidea, Tetraonchidea exhibit a phylogenetic affinity with the Dactylogyridea + Capsalidea clade; and the order Capsalidea is neither basal within the subclass Monopisthocotylea, nor groups with the Gyrodactylidea, but instead forms a sister clade with the Dactylogyridea. The mt genome of P. inermis exhibits a unique gene order, with an extensive reorganization of tRNAs. Monogenea exhibit exceptional gene order plasticity within the Neodermata.


This study shows that gene order within monopisthocotylid mt genomes is evolving at uneven rates, which creates misleading evolutionary signals. Furthermore, our results indicate that all previous attempts to resolve the evolutionary history of the Monogenea may have produced at least partially erroneous relationships. This further corroborates the necessity to generate more molecular data for this group of parasitic animals.


Classification and phylogeny of the Monogenea have been debated for decades without reaching a consensus. Therefore, several classification systems, based on different morphological characters, coexist (e.g. [1,2,3,4,5,6,7]). These often use inconsistent terminology; e.g. Monogenea vs Monogenoidea, Monopisthocotylea vs Polyonchoinea, etc. (for further details see Table 1 in Mollaret et al. [8]). To simplify comparison with previous studies and discussion of our results, we follow the naming system of Boeger & Kritsky [4], with four exceptions: we use order rank for Tetraonchidea (which was split into two parts by Boeger & Kritsky: Tetraonchidae and Amphibdellatidae were assigned to Dactylogyridea, and the remaining Bothitrematidae and Tetraonchoididae were assigned to Gyrodactylidea), Monogenea instead of Monogenoidea, Monopisthocotylea instead of Polyonchoinea, and Polyopisthocotylea instead of Oligonchoinea.

Table 1 Annotated mitochondrial genome of Paratetraonchoides inermis

Paratetraonchoides inermis Bychowsky, Gussev & Nagibina, 1965 (Monogenea: Tetraonchoididae) is a parasitic flatworm usually found on gills of fishes belonging to the family Uranoscopidae. It was originally assigned to the Tetraonchoididae (Tetraonchidea) based on several morphological characteristics: a considerably extended ribbon-shaped body, absence of eyes and middle hooks, and a digestive system typical for this family [9]. However, as the parasite uses 16 marginal gyrodactylid-type (hinged) edge hooks to attach itself to the gills of its host, later it was reassigned to Gyrodactylidea by Boeger & Kritsky [4]. Concerning relationships among the orders Dactylogyridea, Tetraonchidea, Capsalidea and Gyrodactylidea (DTCG group henceforth), Bychowsky [1] argued that the ancestors of the Tetraonchidea were morphologically closer to the Gyrodactylidea than to the Dactylogyridea (note that capsalids were assigned to the Dactylogyridea in this classification). In contrast, in taxonomic classification presented by Lebedev [2] includes the Tetraonchidea, together with the Dactylogyridea and Capsalidea, into the superorder Dactylogyria, whereas the Gyrodactylidea was elevated to the superorder level (Gyrodactylia). This classification was partly supported by Justine [5], based on spermatology data: sperm pattern four was found in the orders Tetraonchidea and Dactylogyridea, whereas sperm patterns 2a and 2c were found in the Capsalidea and Gyrodactylidea, respectively. The fact that these traits produce incongruent results, and that several decades of research have failed to produce a consensus on the phylogeny of DTCG group, presents a strong indication that morphological datasets may not provide a sufficiently robust method to establish a comprehensive phylogenetic hypothesis for DTCG group and the Monogenea.

With rapid advances in sequencing techniques, accompanied by exponentially increasing the availability of molecular data, molecular phylogenetics is becoming the tool of choice for resolving the evolutionary relationships within this group of animals [10, 11]. Despite this, the availability of molecular data for DTCG remains limited and unbalanced (some orders are underrepresented) to infer their phylogenetic relationships with high resolution. For example, currently (July 2017), only 12 molecular records are available in the GenBank database for the Tetraonchidea. The only molecular marker used so far to study the phylogenetics of DTCG was 18S ribosomal RNA [11]. This marker produced a topology somewhat similar to Justine’s results [5]: (Dactylogyridea + Tetraonchidea) + (Monocotylidea + (Capsalidea + Gyrodactylidea)). However, as single-gene markers may not carry a sufficient amount of information to provide high phylogenetic resolution, future studies shall probably increasingly rely on much more powerful phylogenomic approaches [12, 13].

Owing to the abundance of mitochondria in animal tissues, maternal inheritance, the absence of introns, the small size of genomes in metazoans, and an increasingly large set of available orthologous sequences, metazoan mitochondrial (mt) genomes have become a popular tool in population genetics [14], phylogenetics [15, 16] and diagnostics [12]. Comparisons of the arrangements of genes have garnered much scientific attention, as they carry a wealth of complementary resources with potential applications in molecular systematics [17, 18]. Although gene order within neodermatan (i.e. parasitic flatworms) mt genomes is assumed to be remarkably conserved [16, 19, 20], several exceptions have been reported, including African/Indian schistosomes [19, 21], polyopisthocotylids [16, 22, 23] and a single monopisthocotylid [24].

Here, we sequenced, annotated and characterized the first tetraonchid mt genome sequence, belonging to P. inermis. We compared it structurally to all available (July 2017) neodermatan mt genomes and used all available monogenean mt genome sequences to reconstruct the phylogenetic relationships within the entire class.


Specimen collection and DNA extraction

Monogeneans were collected on 27th September 2016 from the gills of Ichthyscopus lebeck (Bloch & Schneider, 1801) obtained from Dong-He market in Zhoushan, Zhejiang Province, China (29°56′–40°62′N; 122°18′–12°30′E). Paratetraonchoides inermis was identified morphologically by the hard parts of the haptor (dorsal and ventral connective bars, marginal hooks) and reproductive organs (male copulatory organ and vaginal armament) [9] under a stereomicroscope and a light microscope. The parasites were preserved in 99% ethanol and stored at 4 °C. Total genomic DNA was extracted from about 120 entire parasites using the TIANamp Micro DNA Kit (Tiangen Biotech, Beijing, China), according to the manufacturer’s recommended protocol, and stored at -20 °C.

Amplification and sequencing

Conducted as described before [25], with minor modifications: partial sequences of nad5, nad1, cox1, cox3, cox2 and rrnS genes were amplified by PCR using six degenerate primer pairs (Additional file 1). Based on these newly sequenced fragments, we designed specific primers for amplification and sequencing of the whole mitogenome (Additional file 1: Table S1). PCR was performed in a 20 μl reaction mixture, containing 7.4 μl dd H2O, 10 μl 2× PCR buffer (Mg2+, dNTP plus, Takara, Dalian, China), 0.6 μl of each primer, 0.4 μl rTaq polymerase (250 U, Takara), and 1 μl of DNA template. Amplification was conducted under the following conditions: initial denaturation at 98 °C for 2 min, followed by 40 cycles at 98 °C for 10 s, 48–60 °C for 15 s, 68 °C for 1 min/kb, and a final extension at 68 °C for 10 min. PCR products were sequenced bi-directionally at Sangon Company (Shanghai, China) using the primer-walking strategy.

Sequence annotation and analyses

After BLASTn analysis [26], the mitochondrial genomic sequence was assembled manually in a stepwise manner. The mt genome was aligned against the mt genomic sequences of other published monogeneans using MAFFT 7.149 [27] to determine approximate gene boundaries. The annotation was further fine-tuned using Geneious [28] adopting one capsalid mt genome, Neobenedenia melleni (MacCallum, 1927) (JQ038228) as the reference, and finally recorded in a Word document. Protein-coding genes (PCGs) were found by searching for ORFs (employing genetic code 9, echinoderm mitochondrial) and checking nucleotide alignments against the reference genome in Geneious. All tRNAs were identified using ARWEN [29], DOGMA [30], and MITOS [31] web servers. The two rRNAs, rrnL and rrnS, were also preliminarily found using MITOS, and their precise boundaries determined by alignment with closely related orthologs in Geneious. The NCBI submission file and tables with statistics for mt genomes were created using a home-made GUI-based program, MitoTool [32]. A nucleotide composition table was then used to make the broken line graph of skewness and A + T content in ggplot2 [33]. Codon usage and relative synonymous codon usage (RSCU) for 12 protein-encoding genes (PCGs) of seven analyzed monopisthocotylids were initially computed with MEGA 5 [34], then further sorted using custom-made Python scripts [35], and finally imported into ggplot2 to draw the RSCU figure. Rearrangement events in studied mt genomes and pairwise comparisons of gene orders of seven monogeneans were analyzed with CREx web tool [36] using the common interval measurement.

Phylogenetic analyses

Phylogenetic analyses were conducted using amino acid sequences of PCGs of the newly sequenced mt genome (P. inermis) and all 17 monogenean mt genomes available in the GenBank (Additional file 2: Table S2). Two species of the order Tricladida, Crenobia alpina (Dana, 1766) (KP208776) and Obama sp. MAP-2014 (NC_026978), were used as outgroups, as suggested in our previous study [25]. A fasta file with nucleotide sequences for all 12 PCGs was extracted from GenBank files and translated into amino acid sequences (employing genetic code 9, echinoderm mitochondrial) using MitoTool. All genes were aligned in batches with MAFFT, integrated into another GUI-based program written by us, BioSuite [37]. BioSuite was also used to concatenate these alignments and generate phylip and nexus format files. Phylogenetic analyses were conducted using maximum likelihood (ML) and Bayesian inference (BI) methods. Selection of the most appropriate evolutionary model for the dataset was carried out using ProtTest [38]. Based on the Akaike information criterion, MtArt + I + G + F was chosen as the optimal model for ML analysis, and LG + I + G + F for the BI analysis. ML analysis was performed in RaxML GUI [39] using a ML + rapid bootstrap (BP) algorithm with 1000 replicates. BI analysis was performed in MrBayes 3.2.6 [40] with default settings, and 5 × 106 metropolis-coupled MCMC generations. Finally, phylograms and gene orders were visualized and annotated by iTOL [41] with the help of several dataset files generated by MitoTool.

Results and discussion

Genome organization and base composition

The circular mitochondrial genome of P. inermis is 14,654 bp in size (GenBank accession number KY856918). The mitogenome is comprised of 12 protein-encoding genes (PCGs), 22 tRNA genes, two rRNA genes, and a major non-coding region (mNCR) (Fig. 1). The genome lacks the Atp8 gene, which is common for flatworms [42]. All genes are transcribed from the same strand. Six overlaps and 22 intergenic regions were found in the genome (Table 1). Noteworthy, the A + T content of the whole genome (82.6%), concatenated PCGs (81.9%), concatenated rRNA genes (81.1%), concatenated tRNA genes (82.6%), and even individual elements (three codon positions, two rRNAs, all 12 individual PCGs) of the genome, are the highest reported among the monogenean mitogenomes characterized so far (Additional file 2: Table S2; Fig. 2). Nucleotide skewness of the mt genome (as well as its elements) did not differ from other monogeneans (Fig. 2).

Fig. 1
figure 1

Visual representation of the circular mitochondrial genome of Paratetraonchoides inermis. Protein-coding genes (12) are red, tRNAs (22) are yellow, rRNAs (2) are green, and non-coding regions are grey

Fig. 2
figure 2

A + T content and skewness of individual elements and the complete genome. Species are coloured according to their taxonomic placement at the order level

Protein-coding genes and codon usage

Coalesced PCGs were 9996 bp in size, with a notably high A + T content of 81.9%. This was also reflected in individual PCGs: from 75.1% (cox1) to 87.9% (nad2 and nad3) (Table 2). Apart from cox1 (which used GTG), ATG was the initial codon for all other PCGs. Among the terminal codons, 11 out of 12 were TAA, whereas nad4L used TAG. No abbreviated stop codons (T--) were found (Table 1).

Table 2 Nucleotide composition and skewness of different elements of the studied mitochondrial genome

Codon usage, RSCU, and codon family proportion (corresponding to the amino acids usage) were investigated among seven monopisthocotylid representatives (Fig. 3). Except for Tetrancistrum nebulosi (Young, 1967), the most abundant codon families were Leu2, Ile, and Phe. This is comparable to Lepidoptera [43] and Nemertea [44]. Noteworthy, the studied mt genome exhibited a strong preference for the A + T-rich members of these four codon families (>10%, Phe, Ile, Leu2 and Asn in Fig. 3), whereas three codons mainly composed of G + C (CGC, GCG and CUG) were not found at all. This all corresponds well to the exceptionally high A + T bias of this mt genome. Overall, A + T-rich codons were favored over synonymous codons with lower A + T content among all seven considered monopisthocotylids (Fig. 3). This A + T preference is notably exemplified by the Leu2 family (TTA and TTG), where the TTA codon accounted for 86.92 ± 4.64%.

Fig. 3
figure 3

Relative synonymous codon usage (RSCU) of seven monopisthocotylid representatives. Codon families are labelled on the x-axis. Values on the top of the bars denote amino acid usage

Transfer RNA genes

All 22 standard tRNAs were found (Table 1), and most of them exhibited the conventional cloverleaf structure. Exceptions were trnS1 (AGN) and trnC, which lacked DHU arms. The unorthodox trnS1 (AGN) is commonplace among all sequenced monogeneans [12, 15, 16, 22,23,24, 45,46,47,48,49,50,51], and possibly even all flatworms [42]. The unpaired DHU-arm in tRNA cys was also reported in most monogeneans apart from Pseudochauhanea macrorchis, M. sebastis, Polylabris halichoeres (Wang & Yang, 1998), T. nebulosi and N. melleni [16, 22, 23, 45, 52]. A slight preference for the A nucleotide (AT skewness, 0.013) was found in concatenated tRNAs of the P. inermis mt genome, which is an exclusive feature among the analyzed monogeneans (Fig. 2), all of which exhibit a preference for the T nucleotide.

Non-coding regions

The major non-coding region (mNCR), 1201 bp in size and located between nad5 and cox3, had a slightly higher A + T content (86%) than other parts of the genome (Table 2). Within the mNCR, there were two minor repetitive regions, both consisting of two repeats, 19 and 16 bp in size. Three tRNA-like cloverleaf structures were found in the mNCR (Additional file 3: Figure S1), among which trnS1-like and trnL1-like sequences contained modified standard anti-codons (ACT and AAG respectively), whereas trnS2-like had a standard TAA anticodon. Average sequence similarity values between the three tRNA-like pseudo-genes and the corresponding functional monogenean tRNA homologs were low (41.45 ± 4.61% trnS1-like, 37.14 ± 3.84% trnL1-like, and 40.32 ± 6.01% trnS2-like), which indicates that they may not be functional. Such tRNA-like sequences were also observed in mNCRs of many lepidopteran insects [53, 54]. These three tRNA-like genes could be a remnant of the tandem-duplication-random-loss (TDRL) process, and the associated heightened rates of substitutions and indels in duplicated genes. A similar hypothesis was put forward by Cameron [55] with regard to the presence of tRNA-like sequences in the mNCRs of many lepidopteran insects [53, 54]. However, due to the limited data we have at disposal, functionality and presence of such tRNA-like sequences in other closely related species of the Tetraonchidea and other monogeneans remain speculative.


Both phylogenetic analysis methods (BI and ML) produced phylograms with concordant branch topologies and high statistical support: all bootstrap support values were ≥ 88, and all Bayesian posterior probabilities were 1.0. Since both phylograms had the same topology, only the latter is shown (Fig. 4). Tree topology indicates the existence of two major clades: subclass Monopisthocotylea (Gyrodactylidea, Capsalidea, Tetraonchidea and Dactylogyridea) and subclass Polyopisthocotylea (Mazocraeidea). The Monopisthocotylea clade was further sub-divided into two clades, (Tetraonchidea + (Dactylogyridea + Capsalidea)) and Gyrodactylidea, both robustly supported (BP/BPP = 88/1 and 100/1, respectively).

Fig. 4
figure 4

Phylogeny of the five orders inferred using concatenated amino acid sequences of 12 protein-coding genes. Scale-bar represents the estimated number of substitutions per site. Star symbol indicates that both methods produced the maximum statistical support value (BP = 100, BPP = 1.0), elsewhere both values are shown above the node as ML/BI. The number and distribution of hooks, number of anchors and spermatozoon patterns for the five orders are given to the right of the figure. The 14/12 + 2 in the first column refer to 14/12 marginal hooks +2 central hooks

Based on the topology obtained in our phylogenetic analysis, the order Tetraonchidea appears to be much more closely related to the Dactylogyridea + Capsalidea clade than to the order Gyrodactylidea. This result is not fully congruent with any of the previously proposed classifications [1,2,3,4,5, 56]. As the entire order Tetraonchidea was represented by a single species of the Tetraonchoididae (P. inermis) in our study, this topology should be interpreted with some caution. This relationship appears to be supported by spermiogenetic and spermatozoal ultrastructural characters [5] (Tetraonchidea and Dactylogyridea possess sperm pattern 4, exhibiting one axoneme; Fig. 4), but as the Gyrodactylidea and Capsalidea both possess two axonemes (patterns 2a and 2c), we can conclude that sperm morphology and mitochondrial phylogenomics produce incongruent signals. Our results also appear to reject the validity of taxonomic grouping of the Dactylogyridea, Tetraonchidea and Gyrodactylidea by the possession of 16 gyrodactylid-type marginal hooks, proposed by Bychowsky [1]. Similarly, the Tetraonchoididae were reassigned to the Gyrodactylidea by Boeger & Kritsky [4] mainly by hinged ‘gyrodactylid’ hooks, which is also incongruent with our results. Our results also do not support the basal phylogenetic position of the Capsalidea proposed by Boeger & Kritsky [4]. Finally, the results are also in disagreement with phylogenetic classifications based on molecular data: 18S ribosomal RNA sequences produced a topology in which capsalids were phylogenetically closer to gyrodactylids than dactylogyrids [11].

Although our phylogenetic framework failed to reach a consensus with any of the previous studies, either those based on morphological or on molecular data, it provides important new insights into the evolutionary history of the four monogenean orders, the Gyrodactylidea, Dactylogyridea, Capsalidea and Tetraonchidea. Morphological traits are often believed to exhibit a high frequency of homoplasy, especially in (parasitic) microscopic animals, whether as a consequence of subjective, or merely simplistic, definitions of a character state (artifact) [57], or of a convergent evolution caused by similar selection pressures on different taxonomic groups [58]. The existence of numerous incompatible phylogenetic hypotheses regarding DTCG group and the entire class Monogenea [1,2,3,4,5,6,7, 56] presents excellent proof that wrong conclusions are often reached when poorly-chosen or numerically insufficient morphological characters are invoked. Although molecular phylogenetics is a promising tool to address this issue [10, 11], future studies should rely on molecular markers that carry a sufficient amount of information to provide high phylogenetic resolution [12, 13].

Gene order

Paratetraonchoides inermis exhibits an extensive reorganization of tRNAs in comparison to all other sequenced monogenean mt genomes (Fig. 5). However, disregarding the tRNA genes, the order of PCGs and rRNA genes within its mt genome conforms to the common neodermatan pattern [20, 59]. A gene order similarity matrix (Additional file 4: Table S3) based on all 36 genes also indicates that P. inermis was the most dissimilar among the compared monopisthocotylids, even including the unique order-possessing A. forficulatus [24]. The transformational pathway from P. inermis to the most similar gene arrangement, found in T. nebulosi and Paragyrodactylus variegatus (You, King, Ye & Cone, 2014), required two coupled transposition events, as well as three coupled long-range rearrangement operations, of which two were a TDRL, and one was a transposition (Additional file 5: Figure S2). Disregarding the two “non-standard” mt genomes, P. inermis and A. forficulatus, the remaining monopisthocotylids exhibit a remarkably conserved gene order [20, 23].

Fig. 5
figure 5

The 21 unique gene orders in neodermatan mitochondrial genomes filtered from 107 species. Representatives and the corresponding taxonomic category at the class/subclass level are shown on the left; star symbol denotes that the gene order is shared by Monogenea and Cestoda. Pattern types used here to classify gene orders are indicated on the right. Also, see Additional file 6: Figure S3)

For a better comparison of gene order among neodermatans, we used MitoTool and iTOL to extract and visualize all available sequences for species of the Cestoda and Trematoda (Additional file 6: Figure S3) and filter out all non-unique gene orders. This resulted in a set of 21 unique gene orders: 11 for the Monogenea, four for the Cestoda (also see our recent publication [60]) and seven for the Trematoda (note that the eleventh gene order in Fig. 5 was shared by the Cestoda and Monogenea). This indicates an exceptional plasticity in the mitochondrial gene order in the Monogenea, as they merely represent 16.8% (18/107) of all available neodermatan mt genomes, but account for 52.4% (11/21) of all unique gene orders. In general, gene order within neodermatan mt genomes is relatively conserved: all of the Cestoda [60] and a majority of the Trematoda mt genomes exhibited only minor variations in the tRNA order (Fig. 5; pattern 1a). Exceptions were only the African/Indian schistosomes (pattern 2, with interchange of two gene blocks: nad5-cox3-cytb-nad4L-nad4 and atp6-nad2, and interchange of nad1 and nad3; Fig. 5) and the monogenean subtaxon Polyopisthocotylea (pattern 3). This was also observed by Webster et al. [19], but they did not have the latest two sequenced monopisthocotylid mitogenomes (P. inermis and A. forficulatus) at their disposal. The majority of Monopisthocotylidea species indeed do exhibit the most common gene order pattern 1a, but these two mt genomes both exhibit extensively altered gene orders: A. forficulatus exhibits pattern 4, with a transposal of an entire gene block [24]; and the newly sequenced P. inermis is the sole representative of the pattern 1b, exhibiting an extensive reshuffling of tRNA genes (Fig. 5).

We hypothesise that the most common pattern (1a) might be the primitive gene order (plesiomorphy) from which patterns 1b, 2, 3 and 4 were derived. This hypothesis is in agreement with previous studies [19, 21, 61], which considered the gene order of Asian schistosomatid mt genomes as the ancestral state, and the gene order of African/Indian schistosomatids a derived trait. Two non-standard gene patterns (4 and 1b) contradict the two hypotheses regarding the putative monogenean plesiomorphic gene arrangement atp6-nad2-trnV-trnA-trnD-nad1-trnN-trnP-trnI [23], and the discriminating markers between Monopisthocotylea and Polyopisthocotylea: rrnL-trnC-rrnS and trnN-trnP-trnI-trnK [48] (Table 3).

Table 3 Occurrence of gene blocks in the five proposed gene order patterns

Such rare gene arrangements are believed to be a promising tool for molecular systematics and phylogenetic reconstruction because mitochondrial gene order reversal events are very rare, and unique orders rarely occur independently in separate lineages [18]. However, the latest two sequenced mt genomes (P. inermis and A. forficulatus) show that monopisthocotylids do not possess a synapomorphic gene order, and instead suggest that gene order within this group may be evolving at uneven rates. This can create misleading evolutionary signals, which was observed before in some taxonomic groups [62,63,64,65]. Thus, while taking into consideration that our insight is curbed by the sparsity of available mt genomes, this finding provides a strong note of caution to researchers wishing to use gene order rearrangements as a tool for neodermatan phylogeny.


Despite the limited availability of molecular data, our analysis provides three findings particularly worth noting. Firstly, there is no support for the sister-group relationship between the Gyrodactylidea and Tetraonchidea [1], nor for the allocation of the family Tetraonchoididae to Gyrodactylidea [4]. Instead, the Tetraonchidea exhibits a phylogenetic affinity with the Dactylogyridea + Capsalidea clade, which indirectly supports Lebedev’s traditional classification [2]. Secondly, the order Capsalidea is neither basal within the subclass Monopisthocotylea [4], nor forms a sister group with the Gyrodactylidea [10, 11]. Instead, it forms a sister clade with the Dactylogyridea, which is in full agreement with the two latest mitochondrial phylogenomic studies [24, 25] and lends further support to the traditional classifications by Bychowsky [1] and Lebedev [2]. Thirdly, the mitogenome of P. inermis provides several interesting findings from the genomic perspective as well: the unprecedentedly high A + T content of the entire genome and its elements, three tRNA-like sequences found in the mNCR, and a unique gene order. The latter indicates that gene order within monopisthocotylids may be evolving at uneven rates, thus creating misleading evolutionary signals. Heightened AT bias can confound phylogenetic inference [66] and the inclusion of only a handful of representatives for three orders (one for Tetraonchidea, two for Dactylogyridea and three for Capsalidea) in our analyses severely limits the phylogenetic resolution. Therefore, we are currently not able to generate a comprehensive phylogenetic hypothesis for the high-level phylogeny of Monopisthocotylea subclass, nor to conduct tests rigorous enough to be able to reject/accept with confidence the hypotheses put forward by the previous studies. Denser sampling and use of strategies alleviating potential compositional biases are needed to evaluate our phylogenetic results and resolve the phylogeny of monogeneans. Our work offsets the scarcity of molecular data for the order Tetraonchidea to some extent, providing a base both for the future fragmentary dataset studies (morphology data and single gene sequence-based molecular markers), as well as the future mitochondrial phylogenomics studies.



Bayesian inference




Bayesian posterior probability


Dactylogyridea, Tetraonchidea, Capsalidea and Gyrodactylidea


Maximum likelihood


Major non-coding region


Protein-encoding genes


Relative synonymous codon usage




  1. Bychowsky BE. [Monogenetic trematodes. Their systematics and phylogeny.] Moscow-Leningrad: Izdatel’stvo Akademiya Nauk SSSR. 1957: 509 pp (In Russian).

  2. Lebedev BI. Monogenea in the light of new evidence and their position among platyhelminths. Angew Parasitol. 1988;29(3):149–67.

    CAS  PubMed  Google Scholar 

  3. Malmberg G. On the ontogeny of the haptor and the evolution of the Monogenea. Syst Parasitol. 1990;17(1):1–65.

    Article  Google Scholar 

  4. Boeger WA, Kritsky DC. Phylogeny and a revised classification of the Monogenoidea Bychowsky, 1937 (Platyhelminthes). Syst Parasitol. 1993;26(1):1–32.

    Article  Google Scholar 

  5. Justine J-L. Cladistic study in the Monogenea (Platyhelminthes), based upon a parsimony analysis of spermiogenetic and spermatozoal ultrastructural characters. Int J Parasitol. 1991;21(7):821–38.

    Article  Google Scholar 

  6. Yamaguti S. Systema helminthum. Vol. IV. Monogenea and Aspidocotylea. New York: Interscience Publishers; 1963.

    Google Scholar 

  7. Llewellyn J. Taxonomy, genetics and evolution of parasites. J Parasitol. 1970;56:493–504.

    Google Scholar 

  8. Mollaret I, Jamieson BG, Justine J-L. Phylogeny of the Monopisthocotylea and Polyopisthocotylea (Platyhelminthes) inferred from 28S rDNA sequences. Int J Parasitol. 2000;30(2):171–85.

    Article  CAS  PubMed  Google Scholar 

  9. Bykhovsky BE, Gussev AV, Nagibina LF. [Monogenetic trematodes of the family Tetraonchoididae Bychowsky, 1951.] Tr Zool Inst Akad Nauk SSSR. 1965;35:140–66 (In Russian).

  10. Olson P, Littlewood D. Phylogenetics of the Monogenea - evidence from a medley of molecules. Int J Parasitol. 2002;32(3):233–44.

    Article  CAS  PubMed  Google Scholar 

  11. Šimková A, Plaisance L, Matějusová I, Morand S, Verneau O. Phylogenetic relationships of the Dactylogyridae Bychowsky, 1933 (Monogenea: Dactylogyridea): the need for the systematic revision of the Ancyrocephalinae Bychowsky, 1937. Syst Parasitol. 2003;54:1–11.

    Article  PubMed  Google Scholar 

  12. Huyse T, Buchmann K, Littlewood DTJ. The mitochondrial genome of Gyrodactylus derjavinoides (Platyhelminthes: Monogenea) - a mitogenomic approach for Gyrodactylus species and strain identification. Gene. 2008;417(1–2):27–34.

    Article  CAS  PubMed  Google Scholar 

  13. Delsuc F, Tsagkogeorga G, Lartillot N, Philippe H. Additional molecular support for the new chordate phylogeny. Genesis. 2008;46(11):592–604.

    Article  PubMed  Google Scholar 

  14. Shao R, Barker S. Mitochondrial genomes of parasitic arthropods: implications for studies of population genetics and evolution. Parasitology. 2007;134:153–67.

    Article  CAS  PubMed  Google Scholar 

  15. Perkins EM, Donnellan SC, Bertozzi T, Whittington ID. Closing the mitochondrial circle on paraphyly of the Monogenea (Platyhelminthes) infers evolution in the diet of parasitic flatworms. Int J Parasitol. 2010;40(11):1237–45.

    Article  CAS  PubMed  Google Scholar 

  16. Park JK, Kim KH, Kang S, Kim W, Eom KS, Littlewood DTJ. A common origin of complex life cycles in parasitic flatworms: evidence from the complete mitochondrial genome of Microcotyle sebastis (Monogenea: Platyhelminthes). BMC Evol Biol. 2007;7:11.

    Article  PubMed  PubMed Central  Google Scholar 

  17. Rokas A, Holland PWH. Rare genomic changes as a tool for phylogenetics. Trends Ecol Evol. 2000;15(11):454–9.

    Article  CAS  PubMed  Google Scholar 

  18. Boore JL, Fuerstenberg SI. Beyond linear sequence comparisons: the use of genome-level characters for phylogenetic reconstruction. Philos Trans R Soc Lond Ser B Biol Sci. 2008;363:1445–51.

    Article  Google Scholar 

  19. Webster BL, Littlewood DTJ. Mitochondrial gene order change in Schistosoma (Platyhelminthes: Digenea: Schistosomatidae). Int J Parasitol. 2012;42(3):313–21.

    Article  CAS  PubMed  Google Scholar 

  20. Wey-Fabrizius AR, Podsiadlowski L, Herlyn H, Hankeln T. Platyzoan mitochondrial genomes. Mol Phylogenet Evol. 2013;69(2):365–75.

    Article  PubMed  Google Scholar 

  21. Littlewood DTJ, Lockyer AE, Webster BL, Johnston DA, Le TH. The complete mitochondrial genomes of Schistosoma haematobium and Schistosoma spindale and the evolutionary history of mitochondrial genome changes among parasitic flatworms. Mol Phylogenet Evol. 2006;39(2):452–67.

    Article  CAS  PubMed  Google Scholar 

  22. Zhang J, Wu X, Xie M, Xu X, Li A. The mitochondrial genome of Polylabris halichoeres (Monogenea: Microcotylidae). Mitochondrial DNA. 2011;22:3–5.

    Article  CAS  PubMed  Google Scholar 

  23. Zhang J, Wu X, Xie M, Li A. The complete mitochondrial genome of Pseudochauhanea macrorchis (Monogenea: Chauhaneidae) revealed a highly repetitive region and a gene rearrangement hot spot in Polyopisthocotylea. Mol Biol Rep. 2012;39(8):8115–25.

    Article  CAS  PubMed  Google Scholar 

  24. Bachmann L, Fromm B, de Azambuja LP, Boeger WA. The mitochondrial genome of the egg-laying flatworm Aglaiogyrodactylus forficulatus (Platyhelminthes: Monogenoidea). Parasit Vectors. 2016;9(1):1.

    Article  Google Scholar 

  25. Zhang D, Zou H, Wu SG, Li M, Jakovlic I, Zhang J, Chen R, et al. Sequencing, characterization and phylogenomics of the complete mitochondrial genome of Dactylogyrus lamellatus (Monogenea: Dactylogyridae). J Helminthol. 2017:1–12.

  26. Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. Basic local alignment search tool. J Mol Biol. 1990;215(3):403–10.

    Article  CAS  PubMed  Google Scholar 

  27. Katoh K, Standley DM. MAFFT multiple sequence alignment software version 7: improvements in performance and usability. Mol Biol Evol. 2013;30(4):772–80.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. Kearse M, Moir R, Wilson A, Stones-Havas S, Cheung M, Sturrock S, Buxton S, et al. Geneious Basic: an integrated and extendable desktop software platform for the organization and analysis of sequence data. Bioinformatics. 2012;28(12):1647–9.

    Article  PubMed  PubMed Central  Google Scholar 

  29. Laslett D, Canback B. ARWEN: a program to detect tRNA genes in metazoan mitochondrial nucleotide sequences. Bioinformatics. 2008;24(2):172–5.

    Article  CAS  PubMed  Google Scholar 

  30. Wyman SK, Jansen RK, Boore JL. Automatic annotation of organellar genomes with DOGMA. Bioinformatics. 2004;20(17):3252–5.

    Article  CAS  PubMed  Google Scholar 

  31. Bernt M, Donath A, Juhling F, Externbrink F, Florentz C, Fritzsch G, Putz J, et al. MITOS: improved de novo metazoan mitochondrial genome annotation. Mol Phylogenet Evol. 2013;69(2):313–9.

    Article  PubMed  Google Scholar 

  32. Zhang D. MitoTool software. 2016. Accessed 23 Oct 2016.

    Google Scholar 

  33. Hadley W. ggplot2: Elegant graphics for data analysis. J Stat Softw. 2010;35(1):65–88.

    Google Scholar 

  34. Tamura K, Peterson D, Peterson N, Stecher G, Nei M, Kumar S. MEGA5: Molecular Evolutionary Genetics Analysis using maximum likelihood, evolutionary distance, and maximum parsimony methods. Mol Biol Evol. 2011;28(10):2731–9.

  35. Zhang D. rscuStack package. 2017. Accessed 17 June 2017.

  36. Bernt M, Merkle D, Ramsch K, Fritzsch G, Perseke M, Bernhard D, Schlegel M, et al. CREx: inferring genomic rearrangements based on common intervals. Bioinformatics. 2007;23(21):2957–8.

    Article  CAS  PubMed  Google Scholar 

  37. Zhang D. BioSuite software. 2016. Accessed 9 Oct 2016.

    Google Scholar 

  38. Abascal F, Zardoya R, Posada D. ProtTest: selection of best-fit models of protein evolution. Bioinformatics. 2005;21(9):2104–5.

    Article  CAS  PubMed  Google Scholar 

  39. Silvestro D, Michalak I. raxmlGUI: a graphical front-end for RAxML. Org Divers Evol. 2011;12(4):335–7.

    Article  Google Scholar 

  40. Ronquist F, Teslenko M, van der Mark P, Ayres DL, Darling A, Höhna S, Larget B, et al. MrBayes 3.2: efficient Bayesian phylogenetic inference and model choice across a large model space. Syst Biol. 2012;61(3):539–42.

    Article  PubMed  PubMed Central  Google Scholar 

  41. Letunic I, Bork P. Interactive tree of life (iTOL) v3: an online tool for the display and annotation of phylogenetic and other trees. Nucleic Acids Res. 2016;44:242–5.

    Article  Google Scholar 

  42. Le TH BD, DP MM. Mitochondrial genomes of parasitic flatworms. Trends Parasitol. 2002;18(5):206–13.

    Article  PubMed  Google Scholar 

  43. Salvato P, Simonato M, Battisti A, Negrisolo E. The complete mitochondrial genome of the bag-shelter moth Ochrogaster lunifer (Lepidoptera, Notodontidae). BMC Genomics. 2008;9:331.

    Article  PubMed  PubMed Central  Google Scholar 

  44. Chen HX, Sun SC, Sundberg P, Ren WC, Norenburg JL. A comparative study of nemertean complete mitochondrial genomes, including two new ones for Nectonemertes cf. mirabilis and Zygeupolia rubens, may elucidate the fundamental pattern for the phylum Nemertea. BMC Genomics. 2012;13:139.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  45. Zhang J, Wu X, Li Y, Zhao M, Xie M, Li A. The complete mitochondrial genome of Neobenedenia melleni (Platyhelminthes: Monogenea): mitochondrial gene content, arrangement and composition compared with two Benedenia species. Mol Biol Rep. 2014;41(10):6583–9.

    Article  CAS  PubMed  Google Scholar 

  46. Huyse T, Plaisance L, Webster BL, Mo TA, Bakke TA, Bachmann L, Littlewood DTJ. The mitochondrial genome of Gyrodactylus salaris (Platyhelminthes: Monogenea), a pathogen of Atlantic salmon (Salmo salar). Parasitology. 2007;134(5):739–47.

    Article  CAS  PubMed  Google Scholar 

  47. Plaisance L, Huyse T, Littlewood DTJ, Bakke TA, Bachmann L. The complete mitochondrial DNA sequence of the monogenean Gyrodactylus thymalli (Platyhelminthes: Monogenea), a parasite of grayling (Thymallus thymallus). Mol Biochem Parasitol. 2007;154(2):190–4.

    Article  CAS  PubMed  Google Scholar 

  48. Ye F, King SD, Cone DK, You P. The mitochondrial genome of Paragyrodactylus variegatus (Platyhelminthes: Monogenea): differences in major non-coding region and gene order compared to Gyrodactylus. Parasit Vectors. 2014;7:377.

    Article  PubMed  PubMed Central  Google Scholar 

  49. Ye F, Easy RH, King SD, Cone DK, You P. Comparative analyses within Gyrodactylus (Platyhelminthes: Monogenea) mitochondrial genomes and conserved polymerase chain reaction primers for gyrodactylid mitochondrial DNA. J Fish Dis. 2016;40(4):541–55.

    Article  PubMed  Google Scholar 

  50. Zou H, Zhang D, Li WX, Zhou S, Wu SG, Wang GT. The complete mitochondrial genome of Gyrodactylus gurleyi (Platyhelminthes: Monogenea). Mitochondrial DNA Part B. 2016;1(1):383–5.

    Article  Google Scholar 

  51. Zhang D, Zou H, Zhou S, Wu SG, Li WX, Wang GT. The complete mitochondrial genome of Gyrodactylus kobayashii (Platyhelminthes: Monogenea). Mitochondrial DNA Part B. 2016;1(1):154–5.

    Article  Google Scholar 

  52. Zhang J, Wu X, Li Y, Xie M, Li A. The complete mitochondrial genome of Tetrancistrum nebulosi (Monogenea: Ancyrocephalidae). Mitochondrial DNA. 2014;27(1):1–2.

    Google Scholar 

  53. Kim MJ, Kang AR, Jeong HC, Kim KG, Kim I. Reconstructing intraordinal relationships in Lepidoptera using mitochondrial genome data with the description of two newly sequenced lycaenids, Spindasis takanonis and Protantigius superans (Lepidoptera: Lycaenidae). Mol Phylogenet Evol. 2011;61(2):436–45.

    Article  PubMed  Google Scholar 

  54. Park JS, Kim MJ, Jeong SY, Kim SS, Kim I. Complete mitochondrial genomes of two gelechioids, Mesophleps albilinella and Dichomeris ustalella (Lepidoptera: Gelechiidae), with a description of gene rearrangement in Lepidoptera. Curr Genet. 2016;62(4):809–26.

    Article  CAS  PubMed  Google Scholar 

  55. Cameron SL. Insect mitochondrial genomics: implications for evolution and phylogeny. Annu Rev Entomol. 2014;59:95–117.

    Article  CAS  PubMed  Google Scholar 

  56. Lambert A. Oncomiracidiums et phylogénèse des Monogènes (Plathelminthes), 2ème partie: Structures argyrophiles des oncomiracidiums et phylogénèse des Monogenea. Ann Parasitol Hum Comp. 1980;55:281–325.

    Article  CAS  PubMed  Google Scholar 

  57. Perkins EM, Donnellan SC, Bertozzi T, Chisholm LA, Whittington ID. Looks can deceive: molecular phylogeny of a family of flatworm ectoparasites (Monogenea: Capsalidae) does not reflect current morphological classification. Mol Phylogenet Evol. 2009;52(3):705–14.

    Article  CAS  PubMed  Google Scholar 

  58. Poulin R, Morand S. The diversity of parasites. Q Rev Biol. 2000;75(3):277–93.

    Article  CAS  PubMed  Google Scholar 

  59. Egger B, Bachmann L, Fromm B. Atp8 is in the ground pattern of flatworm mitochondrial genomes. BMC Genomics. 2017;18(1):414.

    Article  PubMed  PubMed Central  Google Scholar 

  60. Li WX, Zhang D, Boyce K, Xi BW, Zou H, Wu SG, Li M, et al. The complete mitochondrial DNA of three monozoic tapeworms in the Caryophyllidea: a mitogenomic perspective on the phylogeny of eucestodes. Parasit Vectors. 2017;10(1):314.

    Article  PubMed  PubMed Central  Google Scholar 

  61. Wang Y, Wang CR, Zhao GH, Gao JF, Li MW, Zhu XQ. The complete mitochondrial genome of Orientobilharzia turkestanicum supports its affinity with African Schistosoma spp. Infect Genet Evol. 2011;11(8):1964–70.

    Article  CAS  PubMed  Google Scholar 

  62. Babbucci M, Basso A, Scupola A, Patarnello T, Negrisolo E. Is it an ant or a butterfly? Convergent evolution in the mitochondrial gene order of Hymenoptera and Lepidoptera. Genome Biol Evol. 2014;6(12):3326–43.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  63. Dowton M, Austin AD. Evolutionary dynamics of a mitochondrial rearrangement “hot spot” in the Hymenoptera. Mol Biol Evol. 1999;16(2):298–309.

    Article  CAS  PubMed  Google Scholar 

  64. San MD. A hotspot of gene order rearrangement by tandem duplication and random loss in the vertebrate mitochondrial genome. Mol Biol Evol. 2005;23(1):227–34.

    Article  Google Scholar 

  65. Wang J-G, Zhang D, Jakovlić I, Wang W-M. Sequencing of the complete mitochondrial genomes of eight freshwater snail species exposes pervasive paraphyly within the Viviparidae family (Caenogastropoda). PLoS One. 2017; doi:10.1371/journal.pone.0181699.

  66. Talavera G, Vila R. What is the phylogenetic signal limit from mitogenomes? The reconciliation between mitochondrial and nuclear data in the Insecta class phylogeny. BMC Evol Biol. 2011;11(1):315.

    Article  PubMed  PubMed Central  Google Scholar 

Download references


The authors would like to thank Dr. Cong J. Hua for the assistance in sampling the monogeneans. We would also like to thank the editor and the two anonymous reviewers for the time they have invested into reviewing our manuscript.


This work was funded by the Earmarked Fund for China Agriculture Research System (CARS-45-15), the National Natural Science Foundation of China (31572658), and the Major Scientific and Technological Innovation Project of Hubei Province (2015ABA045).

Availability of data and materials

All data are fully available without restriction. The datasets supporting the conclusions of this article are included within the article and its additional files. The mitogenomic sequence is available from the GenBank under the accession number KY856918.

Author information

Authors and Affiliations



WXL and DZ designed the study. DZ, HZ, JZ and RC conducted the experiments. DZ, WXL and IJ conducted the data analysis. DZ and IJ wrote the paper. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Wen X. Li.

Ethics declarations

Ethics approval and consent to participate

All applicable international, national, and institutional guidelines for the care and use of animals were followed.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interest.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional files

Additional file 1: Table S1.

Primers used to amplify and sequence the mitochondrial genome of Paratetraonchoides inermis. NCR is non-coding region. (DOCX 16 kb)

Additional file 2: Table S2.

The list of monogenean species and outgroups used for comparative mitogenomic and phylogenetic analyses. (XLSX 10 kb)

Additional file 3: Figure S1.

Secondary structures of the three tRNA-like sequences found in the major noncoding region. (PDF 216 kb)

Additional file 4: Table S3.

Pairwise comparison of 21 unique mitochondrial DNA gene orders among 107 neodermatan species. The analysis is based on the order of all 36 genes: 12 protein-coding, two rRNA and 22 tRNA genes. Taxa with identical gene order were removed and representative randomly chosen. Scores indicate the similarity between gene orders: the higher the score, the more similar the gene order, where “1254” represents an identical gene order. (XLSX 13 kb)

Additional file 5: Figure S2.

Rearrangement pathway from Paratetraonchoides inermis to Tetrancistrum nebulosi. (PDF 294 kb)

Additional file 6: Figure S3.

Gene orders of 107 neodermatan mitochondrial genomes. Identical gene orders are indicated (and numbered) by a vertical black line on the right. To facilitate the comparison and distance calculation in CREx program, we have re-annotated the mt genomes for which trnS1 (AGN) , S2 (UCN) , L1 (CUN) and L2 (UUR) were ambiguously annotated with the help of ARWEN and MitoTool programs. Mt. genomes for which it was impossible to produce a reliable and consistent annotation were not used in the analysis. In detail: no tRNAs could be predicted by ARWEN for positions 11,564 to 11,635 in Orthocoelium streptocoelium NC_028071; 11,604 to 11,667 of Metorchis orientalis NC_028008; 11,606 to 11,665 of Fasciolopsis buski NC_030528; 11,717 to 11,785 of Fischoederius elongatus NC_028001; and 7304 to 7362 of Gastrothylax crumenifer NC_027833; rrnS was absent from Paragonimus westermani NC_002354, nad2 was absent from Fasciolopsis buski KX449331, nad3 was absent from Artyfechinostomum sufrartyfex KX943545; a duplicated trnS1 (AGN) was found in Metagonimus yokogawai NC_023249, and a duplicated trnC was found in Schistosoma mansoni NC_002545 and all the Schistosoma japonicum isolates. (PDF 4935 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Zhang, D., Zou, H., Wu, S.G. et al. Sequencing of the complete mitochondrial genome of a fish-parasitic flatworm Paratetraonchoides inermis (Platyhelminthes: Monogenea): tRNA gene arrangement reshuffling and implications for phylogeny. Parasites Vectors 10, 462 (2017).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: