Sequencing, characterization and phylogenomics of the complete mitochondrial genome of Dactylogyrus lamellatus (Monogenea: Dactylogyridae)

Abstract

Despite the worldwide distribution and pathogenicity of monogenean parasites belonging to the largest helminth genus, Dactylogyrus, there are no complete Dactylogyrinae (subfamily) mitogenomes published to date. In order to fill this knowledge gap, we have sequenced and characterized the complete mitogenome of Dactylogyrus lamellatus, a common parasite on the gills of grass carp (Ctenopharyngodon idella). The circular mitogenome is 15,187 bp in size, containing the standard 22 tRNA genes, 2 rRNA genes, 12 protein-encoding genes and a long non-coding region (NCR). There are two highly repetitive regions in the NCR. We have used concatenated nucleotide sequences of all 36 genes to perform the phylogenetic analysis using Bayesian inference and maximum likelihood approaches. As expected, the two dactylogyrids, D. lamellatus (Dactylogyrinae) and Tetrancistrum nebulosi (Ancyrocephalinae), were closely related to each other. These two formed a sister group with Capsalidae, and this cluster finally formed a further sister group with Gyrodactylidae. Phylogenetic affinity between Dactylogyrinae and Ancyrocephalinae was further confirmed by the similarity in their gene arrangement. The sequencing of the first Dactylogyrinae, along with a more suitable selection of outgroups, has enabled us to infer a much better phylogenetic resolution than recent mitogenomic studies. However, as many lineages of the class Monogenea remain underrepresented or not represented at all, a much larger number of mitogenome sequences will have to be available in order to infer the evolutionary relationships among the monogeneans fully, and with certainty.

Introduction

Monogenea, a group of largely ectoparasitic members of the flatworm phylum Platyhelminthes, mainly found on skin or gills of fish, is composed of two major subgroups, Polyopisthocotylea and Monopisthocotylea, distinguished by the morphology of attachment organs (Justine, 1998). Sometimes an alternative nomenclature is also used: Polyonchoinea and Heteronchoinea (including Oligonchoinea and Polystomatoinea) (Boeger & Kritsky, 2001). Although some studies use the family name Ancyrocephalidae (Bychowsky & Nagibina, 1978; Lebedev, 1988), it was shown to be a junior synonym for Dactylogyridae, based on both morphological characters (Kritsky & Boeger, 1989) and molecular data (Šimková et al., 2003). With more than 900 nominal species, Dactylogyrus (Dactylogyrinae subfamily) is the largest helminth genus, with a globally widespread distribution (Gibson et al., 1996). They are highly host-specific, commonly infecting only one host species or several congeneric hosts, but mainly cyprinids (Šimková et al., 2004). Dactylogyrus lamellatus (Achmerov, 1952) is a parasite commonly found on the gills of grass carp (Ctenopharyngodon idella Valenciennes, 1844). As in other monogeneans, its single-host life cycle can be completed easily in a closed system (Schäperclaus, 1991), and it can cause high mortality in cultured grass carp fry and fingerlings (Shamsi et al., 2009). Establishment of D. lamellatus parasites on the gill lamellae gives rise to local and general lesions (Molnár, 1971), and may incur ensuing secondary bacterial, viral and fungal infections (Thoney & Hargis, 1991).

Previously, single molecular markers, such as partial mitochondrial rrnS, cox1 and nad2 genes, 18S and 28S rDNA genes, and internal transcribed spacer (ITS 1 and 2) of the rRNA gene, have been used to study the Dactylogyridae family. These studies focused mainly on the phylogeography (Wang et al., 2014), species identification (Borji et al., 2012), co-speciation in host–parasite assemblages (Šimková et al., 2004) and phylogeny of the Dactylogyridae (Šimková et al., 2006). Regarding the phylogeny, poly- and/or paraphyly of the ‘catch-all’ Ancyrocephalidae family (or Ancyrocephalinae subfamily) have been widely discussed (Šimková et al., 2003, 2006; Mendoza-Palmero et al., 2015). As such, single gene sequences may not provide sufficient resolution. Thus, Delsuc et al. (2008) have argued that 18S rRNA (and other single molecular markers) have a limited resolving power for phylogeny, and suggested phylogenomics as a much more powerful tool. Huyse et al. (2008) have also reported that both the V4 hypervariable region of 18S and the ITS rRNA region did not provide sufficient phylogenetic resolution to distinguish between Gyrodactylus thymalli (Zitnan, 1960) and Gyrodactylus salaris (Malmberg, 1957), so they developed a mitogenomic approach for the identification of Gyrodactylus species and strains. Nevertheless, within the Dactylogyridae family, the complete mitogenomic sequence is currently available only for one representative of the Ancyrocephalinae subfamily, whereas such sequences remain unavailable for the entire large Dactylogyrinae subfamily.

Mitogenomes have several useful characteristics: haploidy, relatively high mutation rates, a lack of recombination and a rapidly increasing set of available orthologous sequences (Gissi et al., 2008), so they have been widely used as genetic markers in population genetics (Shao & Barker, 2007), phylogenetics (Park et al., 2007; Perkins et al., 2010) and diagnostics (Huyse et al., 2008). Notably, based on the phylogenetic analysis of the six available monogenean mitogenomes and all other published platyhelminth mitogenomes, Perkins et al. (2010) rejected the Cercomeromorphae theory (an assertion that Monogenea is more closely related to Cestoda than to Trematoda) and found that Monogenea is paraphyletic.

Comparisons of mitochondrial gene arrangements are also emerging as a powerful tool for investigating phylogeny and systematics (Boore & Brown, 1998). They might be particularly useful for resolving the phylogeny of Platyhelminthes, which are undergoing a very fast nucleotide substitution rate (Lavrov & Lang, 2005), which can lead to mutational saturation, whereas the much slower rate of gene rearrangements preserves the phylogenetic signal for much longer periods of time. Indeed, mitochondrial gene order has proven useful for inferring phylogenetic relationships in certain groups of Platyhelminthes, including the relationships between major lineages of monogeneans, trematodes and cestodes (Park et al., 2007), African and Asian schistosomes (Huyse et al., 2007), as well as monopisthocotyleans and polyopisthocotyleans (Zhang et al., 2012).

So far, 16 complete monogenean mitogenome sequences have been published, including 13 monopisthocotylid species (9 Gyrodactylidae, 3 Capsalidae and 1 Dactylogyridae) and 3 polyopisthocotylid species (2 Microcotylidae and 1 Chauhaneidae). In the present study, we have sequenced and characterized the first mitogenome of any representative of the Dactylogyrinae subfamily, D. lamellatus. Furthermore, we have used all 36 genes to infer its phylogenetic relationships with the remaining 16 monogeneans for which a mitogenome sequence is available.

Materials and methods

Specimen collection and DNA extraction

Monogeneans were collected from the gills of grass carp specimens obtained from an earthen pond in Wuhan, Hubei Province, China. Dactylogyrus lamellatus was identified morphologically by the hard parts of the haptor (anchors, dorsal and ventral connective bars, marginal hooks) and reproductive organs (male copulatory organ and vaginal armament) (Gusev, 1985), under a stereomicroscope and a light microscope. The parasites were preserved in 99% ethanol and stored at 4°C. The total genomic DNA of about 120 individual parasites was extracted from the entire parasites using TIANamp Micro DNA Kit (Tiangen Biotech, Beijing, China) according to the manufacturer's instructions, and stored at −20°C. Taxonomic identity of the specimen was further verified by amplifying a fragment of 18S rDNA and the complete ITS1 region using the upstream primer S1 (5′-ATTCCGATAACGAACGAGACT-3′) and downstream primer H7 (5′-GCTGCGTTCTTCATCGATACTCG-3′) (Sinnappah et al., 2001).

PCR and DNA sequencing

Partial sequences of cox1, cob, rrnS, nad5, nad1, cox2, nad4 and cox3 genes were initially amplified via polymerase chain reaction (PCR) using eight primer pairs (see supplementary table S1). Based on these fragments, we have designed specific primers for subsequent PCR amplification (supplementary table S1). PCR reactions were conducted in a 20-μl reaction mixture, containing 7.4 μl double-distilled water (dd H₂O), 10 μl 2 × PCR buffer (Mg²⁺, dNTP plus; Takara, Dalian, China), 0.6 μl of each primer, 0.4 μl rTaq polymerase (250 U, Takara) and 1 μl DNA template. Amplification was performed 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&h;s, 68°C for 1 min/kb; and a final extension at 68°C for 10 min. PCR products were sequenced bidirectionally at Sangon Company (Shanghai, China) using the primer-walking strategy.

Sequence analyses

The complete mitogenomic sequence of D. lamellatus was assembled manually and aligned against the mitogenomic sequences of other published monogeneans using the program MAFFT 7.149 (Katoh & Standley, 2013) to determine approximately the boundaries of genes. Protein-coding genes were inferred with the help of BLAST and ORF Finder tools (both available from the National Center for Biotechnology Information (NCBI)), employing the echinoderm mitochondrial genetic code (Codon table 9), and translated into amino acid sequences in MEGA 5 (Tamura et al., 2011). A majority of the tRNAs were identified using tRNAscan-SE web tool (Lowe & Eddy, 1997), and the rest were found by visual comparison with other monopisthocotylids. rrnL and rrnS were found by alignment with other published monopisthocotylean mitogenomes, and their ends were assumed to extend to the boundaries of their flanking genes. Their secondary structures were predicted by Mfold software (Zuker, 2003). Tandem Repeats Finder (Benson, 1999) was used to identify tandem repeats in the non-coding region (NCRs), and their secondary structures were predicted by Mfold software. Base composition, amino acid composition of protein-encoding genes (PCGs) and codon usage were computed with MEGA 5. Rearrangement events in the mitogenomes and pairwise comparisons of gene orders of seven monogeneans were calculated with the CREx program (Bernt et al., 2007) utilizing the breakpoint dissimilarity measurement. Due to the limitations of the program, the taxa with multiple NCRs were removed from the CREx analysis; these included: Aglaiogyrodactylus forficulatus (Kritsky et al. 2007), Gyrodactylus kobayashii (Hukuda, 1940), G. gurleyi (Price, 1937), G. salaris, G. thymalli, G. derjavinoides (Malmberg et al., 2007), G. brachymystacis (Ergens, 1978), G. parvae (You, Easy & Cone, 2008), Benedenia hoshinai (Ogawa, 1984) and B. seriolae (Yamaguti, 1934). Linear comparison of the 17 studied monogenean mitogenomes was visualized using EasyFig2.2 (Sullivan et al., 2011), with the E-value threshold set to 0.001 for BLASTn. Non-synonymous (dN)/synonymous (dS) mutation rates among the 12 PCGs of D. lamellatus and the only remaining Dactylogyridae species with available mitogenomic sequence, Tetrancistrum nebulosi (Young, 1967), were calculated with KaKs_Calculator (Zhang et al., 2006) using a modified Yang–Nielsen algorithm.

Phylogenetic analyses

Phylogenetic analyses were undertaken on the newly sequenced mitogenome of D. lamellatus and the 16 monogenean genomes available in GenBank (see supplementary table S4). Two Tricladida (order) species, Crenobia alpina (Dana, 1766) (KP208776) and Obama sp. MAP-2014 (NC_026978), were used as outgroups. A Fasta file with the nucleotide sequences for all 36 genes (12 PCGs, 2 rRNAs and 22 tRNAs) was extracted from the GenBank files using a home-made GUI-based program: MitoTool (https://github.com/dongzhang0725/MitoTool). This program was also used to generate the *.sqn file (for GenBank submission) by parsing the Word document's comment box, yielding the original *.csv format files for table 1 and supplementary tables S2, S3 and S4. All the genes were aligned in batches with MAFFT integrated in another home-made GUI-based program, BioSuite (https://github.com/dongzhang0725/BioSuite), wherein codon-alignment mode was used for the 12 PCGs, and normal-alignment mode for the remaining genes (rRNAs and tRNAs). BioSuite was then used to concatenate these alignments into a single alignment and generate phylip and nexus format files for the phylogenetic analyses, conducted using maximum likelihood (ML) and Bayesian inference (BI) methods. The selection of the most appropriate evolutionary models for the dataset was carried out using ModelGenerator v0.8527 (Keane et al., 2006). Based on the Akaike information criterion, GTR + I+G was chosen as the optimal model of nucleotide evolution. ML analysis was performed by RaxML GUI (Silvestro & Michalak, 2012) using the ML + rapid bootstrap (BP) algorithm with 1000 replicates. BI analysis was performed in MrBayes 3.2.1 (Ronquist et al., 2012) with default settings, and 1 × 10⁷ metropolis-coupled MCMC generations.

Results and discussion

Genome organization and base composition

The circular mitogenome of D. lamellatus was 15,187 bp in length (GenBank accession number: KR871673), which is close to the longest monogenean mitogenome (15,527 bp) described so far – Polylabris halichoeres (Wang & Yang, 1998) (supplementary table S4). Apart from lacking the Atp8 gene, which is common in flatworms (Le et al., 2002), the mitogenome contained the standard 36 genes: 22 tRNA genes, 2 rRNA genes and 12 protein-encoding genes (PCGs), as well as a major NCR (fig. 1). All genes were transcribed from the same strand. Similar to other monogeneans (supplementary table S4), it had a high A + T content (70.6%) (supplementary table S2). Nine overlapping regions were found in the genome (table 1), indicating that it is highly compacted, which is typical for flatworm mitogenomes. The overlap between nad4L and nad4, although variable in length, is common in flatworm mtDNAs, with the exception of two Benedenia species, B. hoshinai and B. seriolae, whose nad4L and nad4 genes are separated by a short NCR (Perkins et al., 2010).

Fig. 1. Map of the complete mitochondrial genome of Dactylogyrus lamellatus. PCGs are red, rRNAs green and tRNAs blue. tRNA genes are labelled with the one-letter amino acid code, where L1 = trnL1 (uag), L2 = trnL2 (uaa), S1 = trnS1 (gcu) and S2 = trnS2 (uga). NCR (grey) is the main non-coding region.

Table 1. The annotated mitochondrial genome of Dactylogyrus lamellatus.

Protein-coding genes and codon usage

The total length of the concatenated 12 protein-coding genes was 9931 bp, with an average A + T content of 68.9%, varying from 63.4% (cox2) to 74.5% (nad6) (supplementary table S2). The most frequent start codon was ATG (for seven PCGs), followed by GTG (five genes). Among the terminal codons, six were TAG, five were TAA, and cox2 had an abbreviated stop codon (T–), which seems to be exclusive to D. lamellatus among the monogenean genomes published to date (supplementary table S3). Codon usage, relative synonymous codon usage (RSCU) and codon family proportion (corresponding to the amino acid usage) of the two dactylogyrids (D. lamellatus and T. nebulosi) are presented in fig. 2. Leucine (15.52%), serine (10.64%) and phenylalanine (10.43%) were the most frequent amino acids in the PCGs of D. lamellatus, while glutamine (1.24%), arginine (1.52%) and lysine (1.61%) were relatively scarce (fig. 2). The most frequent codons were TTT (phenylalanine, 10.00%) and TTA (leucine, 9.00%), whereas the CGC codon for arginine was absent. Such a preference for codon and amino acid usage (as in D. lamellatus) was also found in other Gyrodactylids (Huyse et al., 2007, 2008; Plaisance et al., 2007; Ye et al., 2014; Bachmann et al., 2016), capsalids (Perkins et al., 2010; Zhang et al., 2014b) and polyopisthocotylids (Zhang et al., 2012). From fig. 2, one cannot fail to observe that the codons ending in A or T were predominant, which corresponds to the high A + T content of the third coding position of all PCGs in D. lamellatus (76.2%) and T. nebulosi (69.1%) (supplementary table S2).

Fig. 2. Relative Synonymous Codon Usage (RSCU) of the complete mitochondrial genome of Dactylogyrus lamellatus and Tetrancistrum nebulosi. Codon families are labelled on the x-axis. Values on the top of the bars refer to amino acid usage.

The ratios of non-synonymous (dN) to synonymous (dS) substitutions (ω) for all 12 PCGs of D. lamellatus versus T. nebulosi ranged from 0.03 to 0.26 (supplementary fig. S1). All of the PCGs were under negative (purifying) selection (dN/dS < 1), suggesting the existence of functional constraints affecting the evolution of these genes. Among them, functional constraints on nad4L, nad2, nad6 and atp6 genes were the most relaxed, which was also detected in previous monogenean mitogenomic studies (Huyse et al., 2008; Zhang et al., 2012; Ye et al., 2014). On the contrary, the cox1 gene (ω = 0.03 in this study) was repeatedly found to evolve under strong selective pressure (Huyse et al., 2008; Zhang et al., 2012, 2014b; Ye et al., 2014).

Transfer and ribosomal RNA genes

All 22 commonly found tRNAs were present in the mitogenome of D. lamellatus, ranging from 58 bp (trnS1 (gcu)) to 67 bp (trnC (gca), trnS2 (uga) and trnI (gau)) in size, and adding up to 1412 bp in total concatenated length (supplementary table S2). All of the tRNA sequences could be folded into the conventional cloverleaf structure, including an amino-acyl stem of seven nucleotide pairs (ntp), a dihydrouridine (DHU)-stem of 3–4 ntp with a 4–9 nt loop, an anticodon stem of 5 ntp with a loop of 7 nt, and a TΨC stem of 3–6 ntp with a loop of 3–6 nt (Park et al., 2007). trnS1 (gcu), which lacked the DHU arms, was an exception. rrnL and rrnS were 947 and 725 bp in size, with 69.8 and 66.3% A + T content, respectively (supplementary table S2). They were separated only by trnC (gca), whereas trnT (ugu) was found between rrnL and cox1, which is the standard arrangement for monopisthocotylids (Huyse et al., 2007, 2008; Ye et al., 2014; Zhang et al., 2014a, b) (fig. 1). Benedenia seriolae was an exception again, with a positional change of trnT (ugu) (Perkins et al., 2010) (fig. 3B). Predicted secondary structures of the two ribosomal RNA genes are displayed in supplementary fig. S2.

Fig. 3. (A) Phylogeny of the five families inferred from almost complete genomic datasets (36 genes: 12 PCGs, 2 rRNAs and 22 tRNAs), using two Tricladida species as outgroups. Scale bar represents the estimated number of substitutions per site. Bootstrap/posterior probability support values of ML/BI analysis are shown above the nodes. (B) Linear comparison of monogenean genomes (corresponding to tip labels in (A)). The more similar the sequence, the darker the grey blocks between them. tRNAs are coloured blue; protein genes, red; rRNAs, green; and the major NCR, grey.

Non-coding regions

A total of 13 short intergenic regions (1–135 bp) were interspersed within the mitogenome, adding up to a total of 297 bp (table 1). The major non-coding region (NCR), 1926 bp in size and located between nad5 and trnL2 (uaa), had a much higher A + T content (82.6%) than any other part of the mitogenome (supplementary table S2). Within the major NCR there were two highly repetitive regions (HRRs): HRR1 and HRR2, 587 and 416 bp in size, 78.% and 92.3% A + T content, respectively. HRR1 was composed of six tandem repeats, ranging from 31 to 115 bp in length (fig. 4). Repeat units 2 and 3 were identical in nucleotide composition (both 110 bp). In comparison to these two repeat units, unit 1 was 1 bp longer (111 bp) and differed in two nucleotides; unit 4 was identical in size (110 bp) and differed in one nucleotide, and unit 5 was 5 bp longer (115 bp) and differed in seven nucleotides. Repeat unit 6 was considerably truncated (31 bp), with approximately 80 bp missing at the 3′-end. HRR2 also possessed six DNA repeats, rich in A + T. Repeat units 2–5 were identical (78 bp), whereas the truncated units 1 and 6 were 64 and 40 bp long, respectively. The consensus repeat units of HRR1 (2–3) and HRR2 (2–5) were capable of forming stem-loop structures: −9.48 kcal/mol for the former and −12.24 kcal/mol for the latter (fig. 4). Stem-loop secondary structure of the repeat unit was also reported in polyopisthocotylids, and was assumed to be associated with the origin of replication (Park et al., 2007; Zhang et al., 2012).

Fig. 4. Illustration of highly repetitive regions in the major non-coding region of the Dactylogyrus lamellatus mitochondrial genome. ‘HRR1’ and ‘HRR2’ represent two different repeat regions, both containing six repeat units (repeat 1–6). Sequence alignment of the HRR1 repeat units is shown below, and HRR2 is above. The length of each repeat unit is shown on the right of the alignments. Consensus sites are highlighted with exclamation point symbols and non-consensus sites with asterisks. Logos represent the information content of the aligned sequences at a position and the relative frequency of a base. Secondary structure of the consensus repeat units of HRR1 (2–3) is shown below the alignment, and HRR2 (2–5) above the alignment.

Usually, invertebrate mtDNA genomes possess one or two major NCR(s), rich in A + T content, and hence referred to as ‘AT-rich regions’ (Lunt et al., 1998). Tandem repeats (TRs), presumed to be a consequence of strand slippage during replication (Levinson & Gutman, 1987), are common in the major non-coding region of mitogenomes of Cestoda (von Nickisch-Rosenegk et al., 2001; Duan et al., 2015). Zhang et al. (2012) have proposed a hypothesis that, within the Monogenea, monopisthocotylids can be distinguished from polyopisthocotylids by having fewer and smaller (in size) TRs in the major NCR. However, the NCR of D. lamellatus was pronouncedly longer than NCRs in other monopisthocotylid mitogenomes published to date. The structure of the NCR (large in size and containing two HRRs) and the pattern of the two HRRs (rich in iterations) are similar to those observed in a polyopisthocotylid species, Pseudochauhanea macrorchis (Zhang et al., 2012), which indicates that we can reject that hypothesis. Due to the variability in the iteration number and polymorphisms within and between populations, TRs have been exploited extensively in population genetics studies in a broad range of animal species, including insects, fish and molluscs (Shui et al., 2008; Liu et al., 2014; Atray et al., 2015). Hence, these tandem repeats may provide a novel molecular marker for further studies on systematics and population genetics of D. lamellatus.

Phylogeny and gene order

Both phylograms (BI and ML) had very high statistical support for the branch topology: with one exception (99), all bootstrap support values were 100 and all Bayesian posterior probabilities (BPP) were 1.0. Since both phylograms had identical topology, only the latter is shown (fig. 3A). The tree topology indicated the existence of two major clades: Monopisthocotylea subclass (Gyrodactylidae, Capsalidae and Dactylogyridae) and Polyopisthocotylea subclass (Microcotylidae and Chauhaneidae). Monopisthocotylea subclass was also divided into two major clades: one containing only the Gyrodactylidae family, and the other containing Capsalidae and Dactylogyridae families. This supports early taxonomic classifications: Dactylogyridae and Capsalidae were originally placed together within the order Dactylogyridea, whereas Gyrodactylidae was placed within Gyrodactylidea, and later reassigned to the corresponding superorders Dactylogyria and Gyrodactylia (Lebedev, 1988). These classifications mostly relied on morphological characteristics, such as the number of edge hooks, of which there are 14 in Dactylogyridae and Capsalidae, and 16 in Gyrodactylidae. Several other characteristics specific for Gyrodactylidae further support the hypothesis that this family is evolutionarily distant from the Capsalidae and Dactylogyridae cluster: viviparity, the absence of a larval stage, the absence of eyes and the corona of chitinous hooks placed within the copulatory organ (Bychowsky, 1957).

The results of our phylogenetic analysis are in agreement with the genetic relationships inferred from previous studies that used mitogenomic data to reconstruct the phylogenetic relationships among monogeneans (Bachmann et al., 2016; Zhang et al., 2016; Zou et al., 2016). Remarkably, all (but one) nodes of our phylogenetic tree had the maximum possible statistical support (BP = 100, BPP = 1.0), which was not the case in a very recent study that also used 36 concatenated genes for the analysis (Bachmann et al., 2016). This might be a consequence of the selection of outgroups: whereas Bachmann et al. (2016) used two Cestoda species, we used two Tricladida (Turbellaria) species. To test this hypothesis, we have replaced our two outgoup sequences with these two Cestoda sequences and conducted another phylogenetic test on our dataset: the topology was identical, but statistical support values were somewhat lower (supplementary fig. S3). In order to infer the correct topology, an outgroup should be a taxonomic group that diverged before the last common ancestor of the ingroup (Pearson et al., 2013). As Cestoda appears to be younger, and Turbellaria older, than Monogenea (Park et al., 2007; Perkins et al., 2010), rooting with Tricladida (Turbellaria) species should result in a topology more closely resembling the actual evolutionary history of these taxonomic groups. Most importantly, as none of these studies had a Dactylogyrinae mitogenome sequence at their disposal, our results markedly improve the resolution of monogenean relationships.

Although several researchers suggested that Ancyrocephalinae should be elevated to the family level (Bychowsky & Nagibina, 1978; Lebedev, 1988), our phylogenetic results indicate a sister-group relationship between Ancyrocephalinae and Dactylogyrinae with maximum support value, which adds support to the revision of Kritsky & Boeger (1989). Beyond these, mtDNA gene order and transformational pathway analysis in this study has also provided additional supporting evidence for the close relationship between these two subfamilies. Pairwise comparisons of the mitochondrial gene order among seven taxa revealed that, based on the dissimilarity value (breakpoint algorithm), D. lamellatus was the most similar to T. nebulosi. As expected, D. lamellatus was the most dissimilar from the three Polyopisthocotylea species: Microcotyle sebastis (Goto, 1894), P. halichoeres and P. macrorchis (table 2), mostly due to the pronounced differences in gene order at the cox3–nad5 junction between monopisthocotylids and polyopisthocotylids (Zhang et al., 2014b) (also see fig. 3B). The transformational pathway from D. lamellatus to its sister taxon T. nebulosi requires one transposition event; a single transposition and two coupled transposition events are required to Neobenedenia melleni (MacCallum, 1927); and two coupled long-distance transpositions are required to Paragyrodactylus variegatus (You et al. 2014) (supplementary fig. S4). trnL2 (uaa) abutted upstream of the trnG (ucc) in the mitogenome of D. lamellatus, but in other monopisthocotylids it is usually located between trnS2 (uga) and trnR (ucg). Aglaiogyrodactylus forficulatus genome, where trnL2 (uaa) is located between trnR (ucg) and nad5, is another exception. Apart from the position of trnL2 (uaa), the gene order of D. lamellatus was identical to that of T. nebulosi. From N. melleni, it differed with respect to the position of NCR and trnQ (uug) (fig. 3B).

Table 2. Pairwise breakpoint comparison of mitochondrial DNA gene orders among seven monogenean species, based on the order of all 37 elements (36 genes plus NCR). Taxa with multiple NCRs were not used for the analysis. Scores indicate the dissimilarity between gene orders, where ‘0’ represents an identical gene order. N is the number identifying the taxon in the pairwise comparisons.

In conclusion, we have sequenced and analysed the complete mitogenome of D. lamellatus, which contains the standard 22 tRNA genes, 2 rRNA genes and 12 protein-encoding genes (PCGs). There are two highly repetitive regions in the NCR. Using concatenated nucleotide sequences of all 36 genes (12 PCG, 2 rRNA and 22 tRNA genes), the phylogenetic analysis performed using Bayesian inference and maximum likelihood methods revealed that the two dactylogyrids, D. lamellatus (Dactylogyrinae) and T. nebulosi (Ancyrocephalinae), are very closely related to each other. These two then form a sister group with Capsalidae, and finally, this cluster further forms a sister group with Gyrodactylidae. This result is in agreement with some early traditional classifications and is further supported by morphological characteristics. The phylogenetic affinity between Dactylogyrinae and Ancyrocephalinae is further confirmed by the similarity in their mitochondrial gene arrangement. As many lineages of the class Monogenea are still underrepresented (such as Dactylogyridae and Chauhaneidae), or not represented at all (such as Diplectanidae, Monocotylidae and Diplozoidae), our results do not provide a comprehensive phylogenetic hypothesis for the Monogenea. In order to resolve the evolutionary relationships among the monogeneans with confidence, a much larger number of mitogenomic sequences will have to be available. Hence, the publication of this mitogenome will lend support to future molecular, evolutionary and population studies of D. lamellatus and related monogenean parasites.