Complete Mitochondrial Genome of Two Ectoparasitic Capsalids (Platyhelminthes: Monogenea: Monopisthocotylea): Gene Content, Composition, and Rearrangement

The capsalid monogeneans are important pathogens that generally infect marine fishes and have a substantial impact on fish welfare in aquaculture systems worldwide. However, the current mitogenome information on capsalids has received little attention, limiting the understanding of their evolution and phylogenetic relationships with other monogeneans. This paper reports the complete mitochondrial genomes of Capsala katsuwoni and Capsala martinieri for the first time, which we obtained using a next-generation sequencing method. The mitogenomes of C. katsuwoni and C. martinieri are 13,265 and 13,984 bp in length, respectively. Both species contain the typical 12 protein-coding genes, 2 ribosomal RNA genes, 22 transfer RNA genes, and a control region. The genome compositions show a moderate A+T bias (66.5% and 63.9% for C. katsuwoni and C. martinieri, respectively) and exhibit a negative AT skew but a positive GC skew in both species. One gene block rearrangement was found in C. katsuwoni in comparison with other capsalid species. Instead of being basal to the Gyrodactylidea and Dactylogyridea or being clustered with Dactylogyridea, all species of Capsalidea are grouped into a monophyletic clade. Our results clarify the gene rearrangement process and evolutionary status of Capsalidae and lay a foundation for further phylogenetic studies of monogeneans.


Introduction
Monogeneans are a class of important helminth parasites that are commonly found on the gills, skin, branchiostegal membranes, or buccal cavities of marine and freshwater fishes [1]. Monogeneans have direct life cycles and multiple reproductive strategies, resulting in easy propagation in natural and aquaculture environments [2]. Poor fish health, retarded fish growth, and reduced value of the market product are often caused by monogenean flukes feeding on the cells and mucus of host fish [3]. The Capsala, which was established by Bosc (1811), is a genus of the family Capsalidae Baird, 1853, which comprises 9 subfamilies, 48 genera, and about 200 species [4,5]. Currently, 22 species of Capsala have been recognized [6], among which C. martinierei is considered to be the largest (up to 3 cm) known monogenean species.
The evolutionary relationships in monogeneans remain debatable. For example, on the basis of comprehensive morphological characteristics, Capsalidea was considered basal to the Gyrodactylidea and Dactylogyridea [7]. However, it was inferred to be phylogenetically closely related to the Gyrodactylidea and Dactylogyridea based on the evidence of the 28S rRNA [8], 18S rRNA [9], and mitochondrial genes [10]. For the family Capsalidae, its monophyly is currently supported by the presence of accessory sclerites or modified hooklets on the haptor, providing a synapomorphy for the family [11]. Despite being studied for nearly 230 years [4], many aspects of the classification, systematics, validity, and phylogenetic position (for most species) of the capsalid species taxon are still not resolved. For example, the genera Neobenedenia and Benedenia, which belong to the subfamily Benedeniinae, were revealed to be paraphyletic categories, by analyzing the large subunit rDNA sequence data [5].
Though the characterization of Capsala dates back to 1811, this genus was not recognized in China until 2019 [12]. C. martinierei was first reported on the surface of Mola mola in Chili, with the following morphological characteristics: haptoral accessory sclerites absent; dorsomarginal body sclerites consist of multiple scattered bicuspid and multicuspid sclerites [4]. C. katsuwoni was first recorded on the gills of Katsuwonus pelamis and was characterized by a haptor diameter accounting for approximately 20% of the body's length, haptoral accessory sclerites being bifid at one end, and the presence of a small finger-like fringe on the edge of the haptoral marginal valve [4][5][6]. Despite the clarity of the classification and morphology for most capsalid species, only a limited number of genomic resources exist for this family, which limits our understanding of its molecular evolution and phylogenetic relationships.
The mitochondrion is a fundamental place for respiration and energy production and thus plays an important role in cell metabolism [13]. Owing to the abundance of mitochondria in animal tissues, maternal inheritance, and the absence of introns, mitochondrial genomes have become a powerful tool in population genetics, phylogenetics, and diagnostics [14][15][16]. The mt genomes of monogeneans usually have a set of 36 genes, including 12 protein-coding genes, 22 transfer RNA (tRNA) genes, and 2 ribosomal RNA (rRNA) genes, which are organized and oriented in different ways [17]. This diversity makes them a valuable alternative tool for species identification and phylogenetic studies at the genomic level. Several mt genomes of capsalids have been reported including Capsala pricei [18], Benedenia Diesing [19], Neobenedenia melleni [20], and Capsaloides cristatus [21]. The phylogenetic position of the Capsalidae was deduced by either morphological characters or molecular markers.
In this study, the complete mitochondrial genomes of two capsalid monogeneans on fish found in the South China Sea were newly sequenced. We described the details of genome assembly, annotations, codon usage, and amino acid usage. The available complete mitogenomes of monogenean species, retrieved from the GenBank database, provided insight into the phylogenetic relationship between capsalids and other monogenean species. These results will help us better understand the gene arrangement features of monogenean mitogenomes and lay the foundation for further phylogenetic study of this highly diverse flatworm group.

Sampling and DNA Extraction
During two fishery surveys in the South China Sea, 10 specimens of C. martinierei were collected from the skin of Mola mola Linnaeus, 1758 on 10 February 2021 in the Zhongsha Sea area (13 • 52 N, 110 • 98 E), and 22 specimens of C. katsuwoni were collected from the gills of Auxis thazard (Lacepède, 1800) on 15 May 2018 in the Nansha Sea area (9 • 30 N, 114 • 00 E). After collection, the specimens of the two capsalid parasites were preserved in 95% ethanol at −20 • C for long-term storage. Whole-genome DNA was extracted from one specimen for each parasite species using a Steady Pure Universal Genomic DNA Extraction Kit (Accurate Biotechnology, Changsha, China) following the manufacturer's instructions. DNA quality was evaluated through electrophoresis in a 1% agarose gel and a NanoPhotometer ® spectrophotometer (IMPLEN, Westlake Village, CA, USA).

Mitochondrial Genome Sequencing and Assembling
The library was constructed using a TruseqTM RNA sample Prep Kit (Illumina, San Diego, CA, USA) with 1 µg DNA, which was fragmented into 300-500 bp by a Covaris M220. The library was then sequenced using an Illumina Hiseq platform. High-quality clean data were obtained by filtering out low-quality reads and duplicated reads. Clean data were assembled into optimal contigs by the de novo assembler, NOVOPlasty (https: //github.com/ndierckx/NOVOPlasty, accessed on 10 April 2021). The gene map was generated with the online program OGDraw v1.2 [22].

Gene Rearrangement Analysis and Phylogenomic Reconstruction
The gene rearrangements and phylogenomic analysis were analyzed by comparison of the two newly sequenced mitochondrial genomes of Capsala with 24 species of monogeneans retrieved from GenBank, including 8 species of Gyrodactylidae, 4 species of Capsalidae, 2 species of Tetraonchoididae, 2 species of Diplectanidae, 2 species of Ancylodiscoididae, 4 species of Ancyrocephalidae, and 2 species of Dactylogyridea (Table 1). Schistosoma japonicum from Diplostomida (Digenea) was used as an outgroup species for phylogenomic analysis.
MAFFT was used to perform the sequence alignment [24], and then the data sets were trimmed by trimAl [25]. The concatenated set of nucleotide sequences was used for phylogenetic analysis, which was performed with the BI and ML methods using MrBayes v3.2.6 [26]. The optimal evolution model was GTR + I + G in the jModelTest v2.1.7 [27], and the maximum-likelihood method was used to infer the phylogenetic relationship with 1000 bootstrap replicates in MEGA 6.0 [28]. Bayesian inference (BI) was performed using Mrbayes v3.2 [29]. In the Markov chain Monte Carlo (MCMC) analysis, we ran 1 × 10 8 generations. Samples were taken every 1000 generations, and the first 25% were discarded as burn-in. The stationarity was achieved when the average standard deviation of the splitting frequency remained below 0.01. The resulting phylogenetic trees were visualized in FigTree v1.4.2.

Genome Organization
The mitochondrial genomes of C. katsuwoni and C. martinierei are typical circular molecules that are 13,265 and 13,984 bp in length, respectively ( Figure 1). They were deposited in GenBank under accession numbers OL884727 and OL790148, respectively. The length of the C. martinierei mitogenome is close to the boundaries of those previously reported in Capsalidae species (13,270 bp for Neobenedenia melleni to 13,948 bp for Capsaloides cristatus) [20,21]. The mitogenomes of both studied species contain 12 protein-coding genes (PCGs), 2 rRNAs, and 22 tRNAs, which is in accordance with those of other genera in Capsalidae [18]. In accordance with previously studied monogenean species [20,30], the atp8 gene is absent. There are slight differences between C. katsuwoni and C. martinierei in the sizes of PCGs (10,026 vs. 9966 bp), rRNAs (1662 vs. 1665 bp), and tRNAs (1429 vs. 1436 bp, Table 2). Their overlapping sequences and intergenic sequences are also slightly different. C. katsuwoni has 7 overlapping sequences ranging from 1 to 36 bp and 20 intergenic sequences ranging from 1 to 80 bp. On the other hand, C. martinierei has 7 overlapping sequences ranging from 1 to 35 bp and 18 intergenic sequences ranging from 1 to 413 bp. The longer genomes of C. martinierei should be ascribed to two long intergenic fragments between trnV-trnN and trnN-trnA.

Genome Organization
The mitochondrial genomes of C. katsuwoni and C. martinierei are typical circular molecules that are 13,265 and 13,984 bp in length, respectively ( Figure 1). They were deposited in GenBank under accession numbers OL884727 and OL790148, respectively. The length of the C. martinierei mitogenome is close to the boundaries of those previously reported in Capsalidae species (13,270 bp for Neobenedenia melleni to 13,948 bp for Capsaloides cristatus) [20,21]. The mitogenomes of both studied species contain 12 protein-coding genes (PCGs), 2 rRNAs, and 22 tRNAs, which is in accordance with those of other genera in Capsalidae [18]. In accordance with previously studied monogenean species [20,30], the atp8 gene is absent. There are slight differences between C. katsuwoni and C. martinierei in the sizes of PCGs (10,026 vs. 9966 bp), rRNAs (1662 vs. 1665 bp), and tRNAs (1429 vs. 1436 bp, Table  2). Their overlapping sequences and intergenic sequences are also slightly different. C. katsuwoni has 7 overlapping sequences ranging from 1 to 36 bp and 20 intergenic sequences ranging from 1 to 80 bp. On the other hand, C. martinierei has 7 overlapping sequences ranging from 1 to 35 bp and 18 intergenic sequences ranging from 1 to 413 bp. The longer genomes of C. martinierei should be ascribed to two long intergenic fragments between trnV-trnN and trnN-trnA. Abbreviations of PCGs are: cox1-3 for cytochrome oxidase subunits 1-3; nad1-6 and nad4L for NADH dehydrogenase subunits 1-6 and 4 L, respectively; cob for cytochrome b; atp6 for ATP synthase subunits 6; and rrn12 and rrn16 for small and large rRNA subunits, respectively. Both photographs in the genome maps are ventral view of the whole monogenean bodies. Map A is a microscopy photo of one stained specimen of C. katsuwoni, and map B is a hand-painted image of C. martinierei colored with Photoshop. Protein-coding genes (PCGs) are color-coded (cox: light green, nad: dark green, cob: gray-green, atp: yellow), tRNA genes are in blue, and rRNA genes are in red. Abbreviations of PCGs are: cox1-3 for cytochrome oxidase subunits 1-3; nad1-6 and nad4L for NADH dehydrogenase subunits 1-6 and 4 L, respectively; cob for cytochrome b; atp6 for ATP synthase subunits 6; and rrn12 and rrn16 for small and large rRNA subunits, respectively. Table 2. Nucleotide composition and skewness comparison of different elements of the mitochondrial genomes of Capsala katsuwoni and Capsala martinieri. The nucleotide distributions are also different between C. martinierei and C. katsuwoni. Namely, C. martinierei contains less T than C. katsuwoni, while the former has a higher G+C content than the latter. The GC skewness of mitochondrial genomes is slightly positive (0.13 and 0.06), while the AT skewness is negative (−0.24 and −0.19) for C. katsuwoni and C. martinierei, respectively (Table 2). Compared with previously reported capsalid monogeneans, the A+T contents of the two newly sequenced Capsala species are relatively low, while their AT skewness is similar to that of Benedenia hoshinai (Figure 2, Tables S1 and S2). In most monogenean mitogenomes, the strand skew biases are found to have a negative AT skew and positive GC skew [31]. The strand skew biases of monogenean mitogenomes varies between −0.45 and −0.01, and 0.05 and 0.50 for AT and GC, respectively [32]. For both sequenced mitogenomes of the present study, the PCGs exhibited the highest negative AT skewness in comparison with tRNAs and rRNAs, in accordance with those of other previously recorded monogenean species [33].  The nucleotide distributions are also different between C. martinierei and C. katsuwoni. Namely, C. martinierei contains less T than C. katsuwoni, while the former has a higher G+C content than the latter. The GC skewness of mitochondrial genomes is slightly positive (0.13 and 0.06), while the AT skewness is negative (−0.24 and −0.19) for C. katsuwoni and C. martinierei, respectively (Table 2). Compared with previously reported capsalid monogeneans, the A+T contents of the two newly sequenced Capsala species are relatively low, while their AT skewness is similar to that of Benedenia hoshinai (Figure 2, Tables S1 and S2). In most monogenean mitogenomes, the strand skew biases are found to have a negative AT skew and positive GC skew [31]. The strand skew biases of monogenean mitogenomes varies between −0.45 and −0.01, and 0.05 and 0.50 for AT and GC, respectively [32]. For both sequenced mitogenomes of the present study, the PCGs exhibited the highest negative AT skewness in comparison with tRNAs and rRNAs, in accordance with those of other previously recorded monogenean species [33].

Protein Coding Genes and Codon Usage
The 12 protein-coding genes, accounting for 75.6% and 71.3% of the whole mitochondrial genome, encode 2998 and 3140 amino acids of C. katsuwoni and C. martinierei, respectively (Tables 2, S3 and S4). All protein-coding genes are coded on a majority strand (J-strand). The initiation codons and termination codons are different between species (Table 3). There are three initiation codons (ATG, ATA, and GTG) and three termination codons (TAA, TAG, and TGA) in the mitochondrial genome of C. martinierei, while only two initiation codons (ATG and GTG) and two termination codons (TAA and TAG) were found in the mitochondrial genome of C. katsuwoni. The most common initiation codon is ATG, which occurs in ten and nine genes in C. martinierei and C. katsuwoni mitogenomes, respectively. The starting codon is ATA in nad4 and GTG in nad3 in the mitogenome of C. martinierei. In the C. katsuwoni mitogenome, nad2, nad4, and nad4 L use GTG as an initiation codon. The most common stop codon is TAG, which occurs in seven genes in the C. martinierei mitogenome. Four genes (cox1, nad6, cox3, and cob) use TAA as the stop codon. The least frequent termination codon in C. martinierei mitogenome is TGA, which was only seen in nad1. In C. katsuwoni mitogenome, half of the genes use TAA as a stop codon, while the other half use TAG.
The most common codon used in protein-coding genes is leucine (Leu1 + Leu2), while arginine is the least used codon in both Capsala species ( Figure 3A). This is different from the Benedenia diesing mitogenome, which uses glutamine the least in protein-coding genes (Baeza, Sepúlveda, et al., 2019). The relative synonymous codon usage (RSCU) values for the 12 PCGs are shown in Figure 3B,C and Tables S3 and S4. The usage of both twoand four-fold degenerate codons is biased toward codons abundant in T or A, which is in accordance with other species of Capsalidae [33].

Transfer RNA Genes and Ribosomal RNA
Like most monogenean mitogenomes, the Capsala mitogenome contains a set of 22 tRNA genes [31][32][33], and all of them are encoded on the majority strand ( Table 2). The tRNA genes range in size from 42 bp (Ser1 in C. martinierei) to 68 bp (Cys in C. martinierei and Phe in C. katsuwoni), and the total lengths were 1429 and 1436 bp in C. katsuwoni and C. martinierei, respectively. The anticodons of all the tRNA genes in two newly sequenced species are consistent with those found in closely related monogenean mitogenomes [20],

Transfer RNA Genes and Ribosomal RNA
Like most monogenean mitogenomes, the Capsala mitogenome contains a set of 22 tRNA genes [31][32][33], and all of them are encoded on the majority strand ( Table 2). The tRNA genes range in size from 42 bp (Ser1 in C. martinierei) to 68 bp (Cys in C. martinierei and Phe in C. katsuwoni), and the total lengths were 1429 and 1436 bp in C. katsuwoni and C. martinierei, respectively. The anticodons of all the tRNA genes in two newly sequenced species are consistent with those found in closely related monogenean mitogenomes [20], with the exception of C. martinierei, which exhibits the anticodon AGC instead of GCT in the trnS1 gene. The locations of the rRNAs are the same in both Capsala species; the large rRNA subunit is between trnT and trnC, while the small rRNA subunit is located close to the large rRNA subunit, between trnC and cox2. The locations of both rRNA subunits are conserved in the previously reported monogenean species [18,20,21], with the exception of two Ancylodiscoididae species (Thaparocleidus asoti and Thaparocleidus varicus) [31].

Gene Rearrangement
For a better comparison of the gene order among monogenean species, we extracted and visualized the previously reported and sequenced mitogenomes for 24 species of Monogenea. This resulted in a set of 13 unique gene orders (Figure 4). In general, gene orders within Capsalidae mt genomes are relatively conserved. The newly sequenced C. katsuwoni genome has an identical gene order to another reported Capsala species (C. pricei). The other species of Capsalidae mt genomes exhibit only minor variations in the tRNA order. In both Benedenia species, the trnQ gene occurs between trnF and trnM, while the same gene is located between trnA and trnD in C. katsuwoni. The gene order of C. katsuwoni is remarkably similar to those of Ancyrocephalidae species, and the transformational pathway from the former to the latter requires only one transposition of trnQ. In comparison with C. katsuwoni, C. martinierei exhibited a significant rearrangement in one gene block: from trnV-trnA-trnQ-trnD-nad1-trnN to trnD-nad1-trnV-trnN-trnA-trnQ. with the exception of C. martinierei, which exhibits the anticodon AGC instead of GCT in the trnS1 gene. The locations of the rRNAs are the same in both Capsala species; the large rRNA subunit is between trnT and trnC, while the small rRNA subunit is located close to the large rRNA subunit, between trnC and cox2. The locations of both rRNA subunits are conserved in the previously reported monogenean species [18,20,21], with the exception of two Ancylodiscoididae species (Thaparocleidus asoti and Thaparocleidus varicus) [31].

Gene Rearrangement
For a better comparison of the gene order among monogenean species, we extracted and visualized the previously reported and sequenced mitogenomes for 24 species of Monogenea. This resulted in a set of 13 unique gene orders (Figure 4). In general, gene orders within Capsalidae mt genomes are relatively conserved. The newly sequenced C. katsuwoni genome has an identical gene order to another reported Capsala species (C. pricei). The other species of Capsalidae mt genomes exhibit only minor variations in the tRNA order. In both Benedenia species, the trnQ gene occurs between trnF and trnM, while the same gene is located between trnA and trnD in C. katsuwoni. The gene order of C. katsuwoni is remarkably similar to those of Ancyrocephalidae species, and the transformational pathway from the former to the latter requires only one transposition of trnQ. In comparison with C. katsuwoni, C. martinierei exhibited a significant rearrangement in one gene block: from trnV-trnA-trnQ-trnD-nad1-trnN to trnD-nad1-trnV-trnN-trnA-trnQ. Conserved gene order is a typical feature of mitogenomes [34], which is the case of Garodyctylidae and Ancyrocephalidae in our analysis. Our results suggested that extensive gene order rearrangement occurred in the Capsala genus, which further confirmed the hypothesis that gene order in monogeneans is evolving at a relatively rapid rate [33]. Several commonly used mechanisms, including duplication-random loss, duplication-nonrandom loss, and recombination [35][36][37], have been proposed to explain the gene rearrangements of mitogenomes. Here, we only observed recombination in C. martinierei mitogenome compared with C. katsuwoni, where trnD-nad1 is translocated upstream of trnV, and trnN is translocated to the position of the trnV and trnA junction. The underlying mechanism of this recombination among Capsala species needs further study. As for the gene order in monogeneans, the evolution of mitogenomic gene order arrangements is Conserved gene order is a typical feature of mitogenomes [34], which is the case of Garodyctylidae and Ancyrocephalidae in our analysis. Our results suggested that extensive gene order rearrangement occurred in the Capsala genus, which further confirmed the hypothesis that gene order in monogeneans is evolving at a relatively rapid rate [33]. Several commonly used mechanisms, including duplication-random loss, duplicationnonrandom loss, and recombination [35][36][37], have been proposed to explain the gene rearrangements of mitogenomes. Here, we only observed recombination in C. martinierei mitogenome compared with C. katsuwoni, where trnD-nad1 is translocated upstream of trnV, and trnN is translocated to the position of the trnV and trnA junction. The underlying mechanism of this recombination among Capsala species needs further study. As for the gene order in monogeneans, the evolution of mitogenomic gene order arrangements is generally continuous, with the exception of the Thaparocleidus genus. Our conclusions should be interpreted with more available monogenean mitogenomes to produce reliable evolutionary signals.

Phylogenetic Analyses
The phylogenetic analysis included twenty-six species from seven families of the Monogenea and one outgroup species (Schistosoma japonicum) (Figure 4). Overall, the monogeneans divide into two clades: one containing the orders of Capsalidea and Dactylogyridea, and the other with only Gyrodactylidea species. This classification was chosen because, from the view of morphology, both Dactylogyridea and Capsalidea have 14 marginal hooklets on the haptor, while Gyrodactylidea has 16 hooklets [38]. This analysis is in agreement with those of previous studies that used mitogenomic data to reconstruct the phylogenetic relationships among monogeneans [39].
The genetic distances between families from the Capsalidea and Dactylogyridea orders are controversial. Zhang et al. reported that two dactylogyrid monogeneans formed a sister group with the three capsalid species, even closer than those species from Diplectanidae and other families of Dactylogyridea [10]. With two newly sequenced Capsala species, our results clearly classify Capsalidae into Capsalidea and Diplectanidae and Dactylogyridae into Dactylogyridea. The closer phylogenetic relationships between the Dactylogyridae and Diplectanidae than either of the two with Capsalidae are supported by several morphological-and molecular-data-based studies, including spermatozoon ultrastructure, comprehensive morphological characters, 28S rRNA, and 18S rRNA [35,40,41]. Hence, our results provided comprehensive phylogenetic relationships among the monogeneans with more confidence.
Additionally, we found C. martinierei clusters with C. pricei, and further forms a sister group with C. katsuwoni, which is inconsistent with the gene orders of these three species. C. katsuwoni shows the same gene order as C. pricei, while the recombination of one gene block occurs in C. martinierei. Therefore, we realized that gene order alone cannot be used to analyze the phylogenetic relationships of monogeneans.

Conclusions
In the present study, we sequenced and analyzed the mitogenomes of two capsalid species, which are common parasites on oceanic fishes. In summary, the mitogenome of C. katsuwoni shows a relatively conserved gene architecture, while extensive gene order rearrangement occurred in C. martinierei. The monophyly of Capsalidea was strongly supported by the phylogenetic analysis based on the PCG data. Our results provided useful information for further understanding the gene rearrangement process, phylogenetics, and evolution of monogenean parasites.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/genes13081376/s1, Table S1: A+T content based on 13 mitogenomic elements of the 26 selected monogeneans; Table S2: AT skewness based on 13 mitogenomic elements of the 26 selected monogeneans; Table S3: Amino acid composition and relative synonymous codon usage of Capsala katsuwoni; Table S4: Amino acid composition and relative synonymous codon usage of Capsala martinierei.

Data Availability Statement:
The genome sequence data that support the findings of this study are openly available in GenBank of NCBI at (https://www.ncbi.nlm.nih.gov, accessed on 13 April 2022/20 April 2022) under the accession numbers OL884727 and OL790148.