Evidence for Adaptive Selection in the Mitogenome of a Mesoparasitic Monogenean Flatworm Enterogyrus malmbergi

Whereas a majority of monogenean flatworms are ectoparasitic, i.e., parasitize on external surfaces (mainly gills) of their fish hosts, Enterogyrus species (subfamily Ancyrocephalinae) are mesoparasitic, i.e., parasitize in the stomach of the host. As there are numerous drastic differences between these two environments (including lower oxygen availability), we hypothesized that this life-history innovation might have produced adaptive pressures on the energy metabolism, which is partially encoded by the mitochondrial genome (OXPHOS). To test this hypothesis, we sequenced mitochondrial genomes of two Ancyrocephalinae species: mesoparasitic E. malmbergi and ectoparasitic Ancyrocephalus mogurndae. The mitogenomic architecture of E. malmbergi is mostly standard for monogeneans, but that of A. mogurndae exhibits some unique features: missing trnL2 gene, very low AT content (60%), a non-canonical start codon of the nad2 gene, and exceptionally long tandem-repeats in the non-coding region (253 bp). Phylogenetic analyses produced paraphyletic Ancyrocephalinae (with embedded Dactylogyrinae), but with low support values. Selective pressure (PAML and HYPHY) and protein structure analyses all found evidence for adaptive evolution in cox2 and cox3 genes of the mesoparasitic E. malmbergi. These findings tentatively support our hypothesis of adaptive evolution driven by life-history innovations in the mitogenome of this species. However, as only one stomach-inhabiting mesoparasitic monogenean was available for this analysis, our findings should be corroborated on a larger number of mesoparasitic monogeneans and by physiological studies.


Background
In contrast to a vast majority (more than 95% in estimation [1]) of other monogenean flatworms, which parasitize on external body surfaces, such as fins, skins, and gills of fish hosts, species belonging to the mesoparasitic genus Enterogyrus (Ancyrocephalinae, Dactylogyridae) parasitize in the digestive system (stomach) of their fish hosts [2,3]. This genus currently contains eight known species [4]. They insert the sharp curved terminal ends of their anchors into the epithelial tissue (stomach wall) [3,5]; thus, causing wounds and chronic morbidity to the host [5]. Their peculiar mesoparasitic lifestyle enables

Specimen Collection and Identification
A. mogurndae was obtained from the gills of Siniperca chuatsi (Basilewsky, 1855) (Centrarchiformes: Sinipercidae) bought at a local market in the Wuhan city, Hubei Province on 17th of January 2017. E. malmbergi was obtained from the stomach of Oreochromis niloticus (Linnaeus, 1758) (Cichliformes: Cichlidae), caught in a freshwater pond in the Sun Yat-sen University (Guangzhou city, Guangdong province). Morphological identification of the latter species was conducted as described previously [3], so further details are not given here, whereas A. mogurndae was identified morphologically by the hard parts of the haptor and reproductive organs, as described in Wu et al. [26]. Additionally, to confirm the taxonomic identity of A. mogurndae, its 28S rRNA gene was amplified using a primer pair designed before [27]: C1 (5 -ACCCGCTGAATTTAAGCAT-3 ) and D2 (5 -TGGTCCGTGTTTCAAGAC-3 ). It shares a very high identity of 99.74% with conspecific homologs available in GenBank: 779/781 identical bp with AY841871. All sampled and identified parasites were first washed in 0.6% saline and then stored in 100% ethanol at 4 • C.

DNA Extraction, Amplification, and Sequencing
Due to a small size of these parasites, to ensure a sufficient amount of DNA for amplification and sequencing, we used two kinds of genomic DNA: mixture DNA (20 parasite specimens) and individual DNA (a single parasite specimen). Both were extracted using TIANamp MicroDNA Kit (Tiangen Biotech, Beijing, China). The mixture DNA was first used to amplify the whole mitogenome. First, we selected 14 monogenean mitogenomes from the GenBank, aligned them using ClustalX [28], and designed degenerate primer pairs (Supplementary Table S1) matching the generally conserved regions of mitochondrial genes. On the basis of these obtained fragments, specific primers were designed using Primer Premier 5 [29], and the remaining mitogenome was amplified and sequenced in several PCR steps (Supplementary Table S1). PCR products were sequenced bi-directionally using both degenerate and specific primers mentioned above on an ABI 3730 automatic sequencer using the Sanger method. All obtained fragments were BLASTed [30] to confirm that the amplicon was the actual target sequence. To address the possibility of intraspecific sequence variation present in the mixture DNA, we inspected all chromatograms for double peaks or any other sign of the existence of two different sequences, and then we used individual DNA and long-range PCR to further verify the obtained sequences. If we found two different sequences, we used the DNA extracted from a single specimen to assemble the mitogenomic sequence, thereby ensuring that the entire sequence belongs to a single specimen.

Sequence Annotation and Analyses
Both mitogenomes were assembled and annotated following the procedure described before [31][32][33][34] using DNAstar v7.1 software [35], and MITOS [36] and ARWEN [37] web tools. The annotated genome was recorded in a Word document, and PhyloSuite software [38] was then used to parse and extract the annotation information, as well as create GenBank submission files and organization tables for mitogenomes. Data for the comparative genomics analyses were also extracted and processed using PhyloSuite. Tandem Repeats Finder [39] was employed to find tandem repeats in the long non-coding regions (LNCR), and their secondary structures were predicted by the Mfold software [40]. Genetic distances (identity) among mitogenomic sequences were calculated with the "DistanceCalculator" function in Biopython using the "identity" model.

Phylogenetic Analyses
Phylogenetic analyses were conducted using the two newly sequenced ancyrocephaline mitogenomes and 28 monogenean mitogenomes available in the GenBank (22/June/2019). Among these species, six polyopisthocotylean monogeneans were used as outgroups (Supplementary Table S2). Amino acid alignment of 12 protein-coding genes was used for the phylogenetic analyses. Best-fit homogeneous model (mtZOA+F+I+G4) was calculated with ModelFinder [41], integrated in PhyloSuite. On the basis of the selected model, phylogenetic analyses were performed using two different algorithms: maximum likelihood (ML) and Bayesian inference (BI). ML analysis was conducted using RAxML [42] using an ML+rapid bootstrap (BS) algorithm with 1000 replicates. Bayesian inference analysis with the empirical MTZOA model was conducted using PhyloBayes (PB) MPI 1.5a [43]. For each analysis, two MCMC chains were run after the removal of invariable sites from the alignment, and the analysis was stopped when the conditions considered to indicate a good run (PhyloBayes manual) were reached: maxdiff <0.1 and minimum effective size >300. Phylograms were visualized and annotated using iTOL [44], with the help of several dataset files generated by PhyloSuite.

Selection Analyses
Based on the tree reconstructed above, mutation pressures acting on the 12 mt OXPHOS genes were tested by comparing the nonsynonymous/synonymous substitution ratios (ω = dN/dS), through which positive (ω > 1), neutral (ω = 1), and purifying (ω < 1) selection can be distinguished. To avoid mutational saturation, which can blur the evolutionary signal at higher taxonomic levels, only monopisthocotyleans monogeneans, excluding Gyrodactylidae, were included in the selection analysis (15 species in total). First, we invoked the CODEML package in PAML 4.7 [45] to conduct a branch-site evolutionary analysis of the 12 PCGs. The modified branch-site model A allowed ω to vary, both among sites in the protein and across branches on the tree. After specifying foreground-lineages (E. malmbergi), the null model fixes ω to 1, whereas the alternative model assumes that sequences may have experienced positive selection (ω > 1). The significance of differences between the two models was assessed using likelihood ratio tests (LRTs). Subsequently, positively selected sites of genes with p < 0.05 in LRTs were further evaluated using the Bayes Empirical Bayes (BEB) [45] analysis with posterior probabilities ≥ 0.95. A mixed-effects model of evolution (MEME) and fixed-effect likelihood (FEL) methods implemented in HYPHY [46] were also used to infer the positively selected sites on the same foreground-lineage. Sites were considered as candidates under positive selection when they met the following conditions: β + > α, significant likelihood ratio test (p < 0.05) in MEME, and p < 0.05 in the likelihood ratio test of FEL.
The optimal 3D structures of proteins were predicted using I-TASSER [47] on the basis of the confidence score (C-score) value. Superimposition, visualization, and manipulation of the 3D structures were conducted using PYMOL [48], and the prediction of protein structure and function using the PredictProtein server [49].

Genome Organization and Base Composition
The complete mitochondrial genomes of E. malmbergi (MN095193) and A. mogurndae (MN095192) are circular molecules of 14,760 bp and 14,107 bp, respectively ( Figure 1). The E. malmbergi mitogenome contains the standard [50] 36 flatworm mitochondrial genes (12 protein-encoding genes (PCGs; atp8 is absent), 22 tRNA genes, and 2 rRNA genes), whereas we managed to identify only 35 genes in the mitogenome of A. mogurndae, from which the trnL2 gene appears to be missing (Table 1 and Figure 1). The PCGs encoded by both mitogenomes mostly used canonical start and stop codons (start: ATG and GTG; stop: TAA and TAG), but the nad2 gene of A. mogurndae uses non-canonical codons: the abbreviated version of TAA (T-) as the stop codon (otherwise relatively common, but unique in our dataset), and it proved difficult to determine the initial codon (Table 1 and Supplementary  Table S3). The architecture and similarity of orthologous sequences of the two newly sequenced mitogenomes are summarized in Table 1. The average sequence similarity of PCGs between the two studied mitogenomes ranged from 43.14% (nad4L) to 66.57% (cytb) ( Table 1). The gene order of E. malmbergi (Figure 1) was identical to that of E. johnii, Tetrancistrum nebulosi, Cichlidogyrus sclerosus and Paragyrodactylus variegatus, whereas the gene order of A. mogurndae was almost identical, with the exception of a transposition of trnK and the missing trnL2.

Non-Coding Regions
The long non-coding region (LNCR) was located between nad5 and trnG genes in both mitogenomes ( Figure 2). The two segments were 983 bp and 1430 bp in size in E. malmbergi and A. mogurndae, respectively. The LNCR of E. malmbergi contained two highly repetitive regions (HRR): HRR1 was composed of four successive tandem repeats (TRs), where repeat units 1-3 were identical (122 bp), whereas unit 4 lost one nucleotide at the end; HRR2 was comprised of 23 successively short TRs, in which units 1-22 were identical (16 bp), whereas the unit 23 was 3 bp shorter at the 3 end (Figure 2). Similarly, we found one HRR in the LNCR of A. mogurndae; it contained five successive long TRs, where repeat units 1-3 were identical (253 bp), whereas unit 4 was longer, with one nucleotide insertion at the 102nd position, and unit 5 exhibited one nucleotide deletion at the 253rd position ( Figure 2). The consensus repeat patterns of the HRR1 in E. malmbergi and the HRR in A. mogurndae are capable of forming stem-loop structures ( Figure 2).

Non-Coding Regions
The long non-coding region (LNCR) was located between nad5 and trnG genes in both mitogenomes ( Figure 2). The two segments were 983 bp and 1430 bp in size in E. malmbergi and A. mogurndae, respectively. The LNCR of E. malmbergi contained two highly repetitive regions (HRR): HRR1 was composed of four successive tandem repeats (TRs), where repeat units 1-3 were identical (122 bp), whereas unit 4 lost one nucleotide at the end; HRR2 was comprised of 23 successively short TRs, in which units 1-22 were identical (16 bp), whereas the unit 23 was 3 bp shorter at the 3′ end ( Figure 2). Similarly, we found one HRR in the LNCR of A. mogurndae; it contained five successive long TRs, where repeat units 1-3 were identical (253 bp), whereas unit 4 was longer, with one nucleotide insertion at the 102nd position, and unit 5 exhibited one nucleotide deletion at the 253rd position ( Figure 2). The consensus repeat patterns of the HRR1 in E. malmbergi and the HRR in A. mogurndae are capable of forming stem-loop structures (Figure 2).

Phylogeny and Selection Analyses
BI and ML produced an identical topology (Figure 3), wherein A. mogurndae was resolved as the basal species within the Dactylogyridae, followed by E. malmbergi. The remaining dactylogyrids were divided into two sister-clades: (E. johni, Dactylogyrus lamellatus) and (T. nebulosi, (C. mbirizei, (C. sclerosus, C. halli))). However, the support values within the Dactylogyridae clade in the ML tree were very low (Figure 3). The subfamily Ancyrocephalinae was rendered paraphyletic by the embedded Dactylogyrinae species, D. lamellatus.

Phylogeny and Selection Analyses
BI and ML produced an identical topology (Figure 3), wherein A. mogurndae was resolved as the basal species within the Dactylogyridae, followed by E. malmbergi. The remaining dactylogyrids were divided into two sister-clades: (E. johni, Dactylogyrus lamellatus) and (T. nebulosi, (C. mbirizei, (C. sclerosus, C. halli))). However, the support values within the Dactylogyridae clade in the ML tree were very low (Figure 3). The subfamily Ancyrocephalinae was rendered paraphyletic by the embedded Dactylogyrinae species, D. lamellatus. With the mesoparasitic E. malmbergi selected as the test branch within the phylogenetic tree described above (only the topology of monopisthocotylean monogeneans, excluding gyrodactylids, was used), branch-site specific analyses in PAML, and MEME and FEL analyses in HYPHY, all detected evidence for positive selection in this branch ( Figure 3, Table 2, and Supplementary Table  S4). Under the branch-site model (BSM) test, cox1, cox2, cox3, nad1, and nad5 were identified as genes under positive selection (p < 0.05; Supplementary Table S5). However, BEB analyses only identified several codons in cox2 and cox3 genes as undergoing positive selection (p-value < 0.05). In comparison, the MEME method identified episodic diversifying positive selection (p-value < 0.05) in eleven PCGs (all PCGs except nad1) (Supplementary Table S4). The FEL method found evidence of pervasive positive diversifying selection in all 12 PCGs except for nad3, which only had 30 sites under negative selection (at p ≤ 0.05; Supplementary Table S4). Intriguingly, nad1 was special again, with only one site under pervasive positive diversifying selection and 137 sites under negative selection (Supplementary Table S4).
As positively selected sites are more reliable when they are supported by two or more methods [51,52], only cox2 and cox3 genes were used in subsequent analyses of positive selection on proteins. Positively selected sites that were consistently identified by three methods (BSM, MEME, and FEL) were mapped onto the predicted 3D structures of Cox2 and Cox3 proteins of E. malmbergi (Figure 4). PredictProtein analysis showed that these positively selected sites were located in, or close to, the functional regions: 13 positively selected sites were located within a helical transmembrane region, and one positively selected site (codon 171 in cox2) was located within a protein-binding region ( Figure 4, Table 2). With the mesoparasitic E. malmbergi selected as the test branch within the phylogenetic tree described above (only the topology of monopisthocotylean monogeneans, excluding gyrodactylids, was used), branch-site specific analyses in PAML, and MEME and FEL analyses in HYPHY, all detected evidence for positive selection in this branch ( Figure 3, Table 2, and Supplementary Table S4). Under the branch-site model (BSM) test, cox1, cox2, cox3, nad1, and nad5 were identified as genes under positive selection (p < 0.05; Supplementary Table S5). However, BEB analyses only identified several codons in cox2 and cox3 genes as undergoing positive selection (p-value < 0.05). In comparison, the MEME method identified episodic diversifying positive selection (p-value < 0.05) in eleven PCGs (all PCGs except nad1) (Supplementary Table S4). The FEL method found evidence of pervasive positive diversifying selection in all 12 PCGs except for nad3, which only had 30 sites under negative selection (at p ≤ 0.05; Supplementary Table S4). Intriguingly, nad1 was special again, with only one site under pervasive positive diversifying selection and 137 sites under negative selection (Supplementary Table S4).
As positively selected sites are more reliable when they are supported by two or more methods [51,52], only cox2 and cox3 genes were used in subsequent analyses of positive selection on proteins. Positively selected sites that were consistently identified by three methods (BSM, MEME, and FEL) were mapped onto the predicted 3D structures of Cox2 and Cox3 proteins of E. malmbergi (Figure 4). PredictProtein analysis showed that these positively selected sites were located in, or close to, the functional regions: 13 positively selected sites were located within a helical transmembrane region, and one positively selected site (codon 171 in cox2) was located within a protein-binding region ( Figure 4, Table 2). Protein binding region 9 Protein binding region 11 Helical transmembrane region 15 Helical transmembrane region 18 Helical transmembrane region 19 Helical transmembrane region 20 Helical transmembrane region 28 Helical transmembrane region 31 Helical transmembrane region 32 Helical transmembrane region 39 44 45 Protein binding region 46 47 49 Helical transmembrane region 51 Helical transmembrane region 53 Helical transmembrane region 54 Helical transmembrane region 57 Helical transmembrane region 60 Helical transmembrane region 61 Helical transmembrane region 62 Helical transmembrane region 64 Helical transmembrane region 66 Helical transmembrane region 68 Helical transmembrane region 69 Helical transmembrane region  71   72    73    74    79    Protein binding region  80    Protein binding region  81   Protein binding region  82    88     91   99   117   118   119   126    128   129   133     136   148    152   154   165     170   171 Protein binding region 177

Discussion
Mitogenomic architecture of E. malmbergi is mostly standard for monogeneans, but the mitogenome of A. mogurndae possesses some intriguing features: the absence of trnL2 gene, low AT content (the lowest among the Monogenea; Supplementary Table S1, Supplementary Figure S1), exceptionally long tandem-repeats within the NCR, and non-canonical start/stop codons of the nad2 gene. On the basis of results reported in other related species [12,53], as a working hypothesis, we propose ATA as the initial codon of nad2 in A. mogurndae. The gene order shared between E. malmbergi and other four monogeneans were previously identified as the plesiomorphic gene order for the subphylum Neodermata [54]. This finding further confirms the hypothesis that identical gene orders found in phylogenetically distant lineages are an indicator of shared ancestry [55]. TRs with high repeat numbers and large size have been reported before in the subclass Monopisthocotylea [32,56], but with 253 bp, the TR in A. mogurndae is the longest reported so far. Along with a few previous studies [32,56], these findings consistently reject the hypothesis that monopisthocotylean monogeneans possess fewer and smaller (in size) TRs in the LNCR than polyopisthocotylean monogeneans [57]. Since the presence of tandem repeats forming a stable secondary structure is often associated with replication origin in mitochondria [12,58,59], it appears likely that the repeat regions found in the LNCR are embedded within the control region.
Results of our phylogenetic analyses are partially consistent with the relationships obtained using 18S rDNA and 28S rDNA, in which subfamilies Dactylogyrinae, Ancylodiscoidinae, and Pseudodactylogyrinae were embedded within the subfamily Ancyrocephalinae [60][61][62][63], thus rendering it paraphyletic. However, this result supports the argument that, regardless of the exact status and interrelationships of these subfamilies, all of them originated within the family Dactylogyridae [64]. As nuclear and mitogenomic data reject the monophyly of Ancyrocephalinae, and even the morphological data do not provide full support for it [64], we can conclude with some confidence that a taxonomic revision shall be needed with regard to this taxon. However, the resolution of our approach is too limited (small number of samples) to officially propose it here, so future studies should address this issue using a larger number of samples and combined datasets (ideally all three data types: nuclear, mitochondrial, and morphological).
Although only a small proportion of monogeneans are mesoparasitic (less than 5% in estimation), they have a wide host range (fishes, mammals, and amphibians), exhibit a wide range of parasitic lifestyles with respect to the exact location in the body where they parasitize (stomach, heart, esophagus, cloaca, urinary bladder, and inside the eyelid), and they form phylogenetically distant clades (different subclasses and families) [1,6]. This indicates multiple independent origins of the mesoparasitic lifestyle within the Monogenea. Due to limited availability of data, here, we only compared the selection pressures between ectoparasitic and stomach-inhabiting mesoparasitic monogeneans. Among the 15 species that we included in the selection analyses, 14 are ectoparasitic monogeneans found on the gills, skin, and fins of fish hosts, whereas E. malmbergi is a mesoparasitic monogenean found in the stomach of fish hosts. We hypothesized that starkly different environments of ectoparasitic and mesoparasitic monogeneans (especially the supposed reduced availability of oxygen in the stomach in comparison to the gills) would require metabolic adaptations, which might be reflected on the mitochondrial genome. Our analyses do show that dozens of codons of mitochondrial OXPHOS genes were subjected to strong positive selection pressures (Table 2 and Supplementary Table S4) in E. malmbergi. Particularly, several positive selection sites in cox2 and cox3 were consistently validated by three different methods (Table 2), which provides strong evidence that these genes may have sustained adaptive evolution in the mesoparasitic E. malmbergi. This was further supported by the evidence that many positively-selected codons were posited in, or close to, the predicted functional regions (helical transmembrane region and protein binding region) in the structure of mitochondrial OXPHOS proteins ( Figure 4 and Table 2). cox2 and cox3, both of which were confirmed to have undergone positive selection in E. malmbergi by three methods, belong to the multi-subunit OXPHOS enzyme complex IV [9]. As this complex is believed to consume a majority of the O 2 [65], this finding is indicative of metabolic adaptations to a comparatively (in relation to gills) hypoxic environment in the stomach in E. malmbergi. However, it remains unclear whether this adaptation is a reflection of an enhanced capacity for oxygen usage or a switch to a predominantly anaerobic metabolism. Some of the genes in other enzymatic complexes were also predicted to exhibit positive selection by one or two methods: Complex I, which is the largest OXPHOS enzymatic complex, complex III, which transfers electrons, and complex V, which mainly synthesizes the ATP [66]. Therefore, there is some evidence for positive selection in all four mitochondrial OXPHOS complexes, which suggests life-history innovation-driven adaptive evolution in the mitogenome of E. malmbergi.
The interpretation of these findings is hampered by the scarcity of physiological studies on these parasitic species. Notably, the exact concentration of oxygen in the stomach of their fish hosts remains unknown, it is unclear whether mesoparasitic species employ aerobic or anaerobic metabolic pathways, and we do not know whether these parasites may be able to gain oxygen from the host's tissues and body fluids. However, as the majority of parasites do not use the oxygen available within the host [25], we can tentatively reject the last hypothesis. Furthermore, in mammals, the oxygen content in the small intestine is about 25% of that in the environment (air) [25], and in humans, the gas in the stomach contains approximately 15% oxygen (as opposed to 21% oxygen in the air) [67]. Therefore we can also safely assume that oxygen availability is significantly different between gills (which are constantly flushed by water) and the stomach environment, which is isolated from the surrounding water in freshwater fish most of the time (apart from when they are ingesting food) [68,69]. It should be noted that O. niloticus is a euryhaline fish that ingests water in high salinity environments, which may increase the availability of oxygen in the stomach. However, studies in a closely related species Tilapia mossambica, which exhibits much higher salinity tolerance than O. niloticus [70], have shown that its ingestion of water is fairly low in a freshwater environment [71]. This implies that any stomach-dwelling parasite of freshwater populations would face limited accessibility of oxygen. Finally, as only one stomach-inhabiting mesoparasitic species was sequenced, we could not establish whether the amino-acid replacements observed in E. malmbergi are unique to this species or fixed in all species within this genus (all are mesoparasitic species).

Conclusions
The existence of several unique architectural characteristics of the mitogenome of A. mogurndae is intriguing and indicates that the sequencing of further mitogenomes belonging to this and closely related genera may be interesting from the aspect of the evolution of mitogenomic architecture in the Monogenea. With regards to the status of the subfamily Ancyrocephalinae, the addition of two new mitogenomes (now seven available in total) did not manage to stabilize the topology, and the subfamily was paraphyletic. Although this is partially in agreement with previous molecular studies, support values for the topology were low, so we cannot draw any conclusions with confidence. More data is needed to resolve the status of this subfamily. Finally, in agreement with our hypothesis, we found evidence for adaptive evolution in the mitogenome of E. malmbergi. Although we tentatively attributed it to the adaptation to the evolution from the ectoparasitic to stomach-inhabiting mesoparasitic lifestyle, where the hypoxic environment of the stomach required adaptive changes in some mitochondrial OXPHOS genes, it remains unclear whether this adaptation is a reflection of an enhanced capacity for oxygen usage or a switch to a predominantly anaerobic metabolism. Also, due to a lack of physiological studies on these species, we cannot exclude other, less obvious explanations for these adaptations. In order to investigate the adaptive evolution of stomach-inhabiting mesoparasitic monogeneans comprehensively, we will have to sequence more data from other mesoparasitic species and simultaneously study their physiological and genetic adaptations to the mesoparasitic lifestyle.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4425/10/11/863/s1, Figure S1: Hierarchical clustering of monogeneans based on their A + T content. The Family of each species is shown on the right. Table S1: Primers used to amplify and sequence mitochondrial genomes of Ancyrocephalus mogurndae and Enterogyrus malmbergi. Table S2: The Monogenea species used for comparative mitogenomic and phylogenetic analyses. Table S3: General statistics for the protein-coding and rRNA genes of the 30 selected monogenean species. Table S4: Codons under positive or negative selection inferred by two HYPHY methods: FEL and MEME. P: positive selection; N: negative selection.