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

Background 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. Results 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. Conclusions 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. Electronic supplementary material The online version of this article (10.1186/s13071-017-2404-1) contains supplementary material, which is available to authorized users.


Background
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.
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 ribbonshaped 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°3 0′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.

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. Proteincoding 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 × 10 6 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) Table S2; Fig. 2). Nucleotide skewness of the mt genome (as well as its elements) did not differ from other monogeneans (Fig. 2).

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).
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%.

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.

Phylogeny
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). 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   [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 "nonstandard" mt genomes, P. inermis and A. forficulatus,  the remaining monopisthocotylids exhibit a remarkably conserved gene order [20,23]. 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).
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.

Conclusions
Despite the limited availability of molecular data, our analysis provides three findings particularly worth cox1 G T rrnL rrnS cox2 H M cox3 C L1 K Y S2 nad6 L2 R nad5 E cytb nad4L nad4 Q F atp6 nad2 V A D nad1 N P I nad3 S1 W cox1 rrnL C rrnS cox2 E nad6 Y L1 S2 L2 R nad5 G cox3 H cytb nad4L nad4 T F Q M atp6 nad2 V A D nad1 N P I K nad3 S1 W cox1 G T rrnL rrnS cox2 M H cox3 C K nad6 Y L1 S2 L2 R nad5 E cytb nad4L nad4 Q F atp6 nad2 V A D nad1 N P I nad3 S1 W cox1 T rrnL C rrnS cox2 nad6 Y L1 S2 L2 R nad5 G cox3 E H cytb nad4L nad4 Q F M atp6 nad2 A D nad1 N P I K nad3 V W S1 cox1 T rrnL C rrnS cox2 E nad6 Y L1 S2 R nad5 L2 G cox3 H cytb nad4L nad4 Q F M atp6 nad2 V A D nad1 N P I K nad3 S1 W cox1 T rrnL C rrnS cox2 nad6 Y L1 S1 L2 R nad5 E G cox3 H cytb nad4L nad4 Q F M atp6 nad2 V A D nad1 N P I K nad3 S2 W cox1 T rrnL C rrnS L1 S2 L2 cox2 E nad6 Y R nad5 G cox3 H cytb nad4L nad4 Q F M atp6 nad2 V A D nad1 N P I K nad3 S1 W cox1 T rrnL C rrnS cox2 nad6 L1 Y S2 L2 R nad5 G cox3 E H cytb nad4L nad4 Q F M atp6 nad2 A D nad1 N P I K nad3 W V S1 cox1 T rrnL C rrnS cox2 nad6 Y L1 S2 L2 R nad5 E G cox3 H cytb nad4L nad4 Q F M atp6 nad2 V A D nad1 N P I K nad3 S1 W cox1 T rrnL C rrnS cox2 E nad6 Y L1 S2 L2 R nad5 G cox3 H cytb nad4L nad4 F Q M atp6 nad2 V A D nad1 N P I K nad3 S1 W cox1 T rrnL C rrnS cox2 E nad6 Y S2 L1 L2 R nad5 G cox3 H cytb nad4L nad4 Q F M atp6 nad2 V A D nad1 N P I K nad3 S1 W cox1 T rrnL C rrnS cox2 E nad6 L1 L2 Y S2 R nad5 G cox3 H cytb nad4L nad4 Q F M atp6 nad2 V A D nad1 N P I K nad3 S1 W cox1 T rrnL C rrnS cox2 nad6 L1 S2 L2 R nad5 G Y cox3 H cytb nad4L nad4 Q F M atp6 nad2 V A D nad1 N E P I K nad3 S1 W cox1 T rrnL C rrnS cox2 E nad6 Y L1 S2 L2 R nad5 G cox3 H cytb nad4L nad4 Q F M atp6 nad2 V A D nad1 N P I K nad3 S1 W cox1 T rrnL C rrnS cox2 nad6 Y L1 S2 L2 R nad5 G E cox3 H cytb nad4L nad4 Q F M atp6 nad2 V A D nad1 N P I K nad3 S1 W cox1 T rrnL rrnS cox2 E nad6 L2 Y L1 V G S2 C R nad5 cox3 H cytb nad4L nad4 Q F M atp6 nad2 A D nad1 I K nad3 W N P S1 cox1 T rrnL C rrnS cox2 E nad6 Y L1 S2 L2 R nad5 cox3 G H cytb nad4L nad4 Q F M atp6 nad2 V A D nad1 N P I K nad3 S1 W cox1 G T rrnL rrnS cox2 M H cox3 K Y R C L1 nad6 S2 L2 nad5 E cytb nad4L nad4 Q F atp6 nad2 V A D nad1 N P I nad3 S1 W cox1 T rrnL C rrnS cox2 E nad6 Y S2 V D A nad1 N P I K nad3 L1 R L2 nad5 G cox3 H cytb nad4L nad4 Q M F atp6 nad2 S1 W cox1 T rrnL C rrnS cox2 E nad6 Y L1 Q M S2 L2 R nad5 G cox3 H cytb nad4L nad4 F atp6 nad2 V A D nad1 N P I K nad3 S1 W cox1 T rrnL C rrnS cox2 nad6 Y L1 S2 N I F atp6 nad2 A L2 R nad5 G cox3 E H cytb nad4L nad4 Q K nad3 D nad1 V P M W S1  Fig. 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) noting. Firstly, there is no support for the sistergroup 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.

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) Abbreviations BI: Bayesian inference; BP: Bootstrap; BPP: Bayesian posterior probability; DTCG: Dactylogyridea, Tetraonchidea, Capsalidea and Gyrodactylidea; ML: Maximum likelihood; mNCR: Major non-coding region; PCGs: Proteinencoding genes; RSCU: Relative synonymous codon usage; TDRL: Tandemduplication-random-loss Table 3 Occurrence of gene blocks in the five proposed gene order patterns Gene block/ Pattern cox1-T-rrnL rrnL-C-rrnS cox2-E-nad6 atp6-nad2-V-A-D-nad1 nad1-N-P-I-K-nad3 indicates that a gene block occurs in the corresponding pattern, whereas "-" indicates that the a gene block is broken up. "P" indicates that the gene block is present in all species exhibiting the pattern 1a, except for trematodes. "Y" with a superscript letter indicates the existence of minor exceptions in the gene block, wherein "Y a " denotes Benedenia seriolae; "Y b " denotes Trichobilharzia regenti, Schistosoma mekongi and Schistosoma japonicum; and "Y d " denotes Brachycladium goliath