Mitochondrial genome of Rhabdosynochus viridisi (Monogenea: Diplectanidae), a parasite of Pacific white snook Centropomus viridis

Abstract

We report the nearly complete mitochondrial genome of Rhabdosynochus viridisi – the first for this genus – achieved by combining shotgun sequencing of genomic and cDNA libraries prepared using low-input protocols. This integration of genomic information leads us to correct the annotation of the gene features. The mitochondrial genome consists of 13,863 bp. Annotation resulted in the identification of 12 protein-encoding genes, 22 tRNA genes and two rRNA genes. Three non-coding regions, delimited by three tRNAs, were found between the genes nad5 and cox3. A phylogenetic analysis grouped R. viridisi with three other species of diplectanid monogeneans for which mitochondrial genomes are available.

Introduction

Monogeneans are parasitic flatworms (Platyhelminthes) found mostly on freshwater and marine fish, although some species can infect aquatic or semi-aquatic sarcopterygians such as lungfish, and amphibians, freshwater turtles and hippopotamuses. This group of parasites is economically important as some species can cause disease and mortality in fish farmed around the world (Schelkle et al., 2009; Ogawa, 2015). The short and direct life cycle of these parasites favours their proliferation within production tanks and cages. Monogeneans of the Diplectanidae family in particular pose a problem for the production of marine fish (Tokşen et al., 2013; Andree et al., 2015). Such is the case for Rhabdosynochus viridisi Montero-Rodríguez, Mendoza-Franco & López Téllez, 2020, a diplectanid species recently described as a threat to the pisciculture of Pacific white snook (Centropomus viridis Lockington, 1877) in Mexico (Montero-Rodríguez et al., 2020; Morales-Serna et al., 2020).

Although there are more than 250 species of diplectanids (Domingues & Boeger, 2008), mitochondrial genomes (mitogenomes) have been reported for only three species: Pseudorhabdosynochus yangjiangensis Wu & Li, 2005 [also known as Laticola paralatesi (Nagibina, 1976)] (Tingbao et al., 2006), Lepidotrema longipenis (Yamaguti, 1934) Kritsky, Jiménez-Ruiz & Sey, 2000 and Lamellodiscus spari Zhukov, 1970 (Zhang et al., 2018a). Mitogenomes are considered useful for reconstructing phylogenies and inferring population history due to haploid maternal inheritance (Boore et al., 2005; Stewart & Larsson, 2014). However, some studies have questioned the basic assumptions about the evolutionary dynamics of mitochondrial DNA (mtDNA), that is clonality, near-neutrality and evolutionary rate constancy. They suggest that mtDNA is not a valid molecular marker of population genetics and phylogenetics, and that data on mtDNA might instead provide insight into many aspects of mitogenome evolution (Galtier et al., 2009). Generating mitogenomes is important for supporting comparative mitogenomics studies at the boundary between evolutionary biology and medicine (Galtier et al., 2009). In the present study, we assembled and annotated the mitogenome of R. viridisi using short-read sequencing, and inferred its phylogenetic relationship with other monogeneans for which mitogenomes are publicly available.

Materials and methods

Monogenean R. viridisi were obtained from the gills of juvenile C. viridis maintained in the Laboratory of Parasitology at the Centro de Investigación en Alimentación y Desarrollo (CIAD), Mazatlán, northwestern Mexico (Morales-Serna et al., 2020). Prior to the present study, 10 samples of monogeneans (n = 10) were randomly taken from ten individual fish contained in a single 400-l tank. These specimens were fixed with formaldehyde and stained with Gomori's trichrome or partially digested with proteinase K to observe their morphological characteristics under a microscope (BX53; Olympus). Monogeneans were identified as R. viridisi according to the method of Montero-Rodríguez et al. (2020). Voucher specimens were deposited in the Colección de Parásitos de Peces del Noroeste del Pacífico (CPPNP-1379) at CIAD-Mazatlán. Rhabdosynochus viridisi morphologically resembles Rhabdosynochus alterinstitus Mendoza-Franco, Violante-Gonzalez & Vidal-Martinez, 2008; however, the former can be distinguished by the presence of tegumental scales and an accessory piece with a proximal bifurcation of the vagina. For the present study, 260 live specimens were isolated from a single fish host, washed in distilled water and then fixed in RNAlater (Thermo Fisher, Waltham, MA). Of these specimens, 60 were first observed for tegumental characteristics and then partially digested with proteinase K to observe the sclerotized structures. All the specimens examined were identified as R. viridisi.

The DNA and RNA were isolated from the remaining 200 individuals of R. viridisi. Reamplification, sequencing, assembling, annotation and phylogenetic analysis were performed as described by Caña-Bozada et al. (2021). However, in the present study, mitochondrial genes of the diplectanids P. yangjiangensis, L. longipenis, and L. spari (GenBank accession numbers JQ038231.1, NC_039617.1 and MH328204.1, respectively) were used as references to obtain DNA contigs and for annotation. The DNA and cDNA assembled using MEGA X software (Kumar et al., 2018) and TRINITY v. 2.8.6 (Grabherr et al., 2011), respectively, failed to find a region within the mitogenome. However, an assembly of the cDNA paired reads obtained with MEGAHIT (Li et al., 2016) (parameter: –k-list 21, 29, 39, 59, 79, 99, 119, 141) and posteriorly reassembled with CAP3 (Huang & Madan, 1999), with an identity cut-off of 95% (setting: –p N > 95), enabled such a missing region to be found, and therefore a nearly complete mitogenome was obtained. The identity and boundaries of protein-coding genes (PCGs) and rRNA were determined using ORF Finder tools [available from the National Center for Biotechnology Information (NCBI)] with the genetic code 9, and a comparison with alignments of the mitochondrial genes of species closely related to R. viridisi (species of the family Diplectanidae; see the codes in supplementary table S1) in MEGA X. OGDRAW v. 1.1 (Greiner et al., 2019) was used to visualize the resulting mitogenome. In addition, coverage of the mitochondrial genome was visualized using Integrative Genomics Viewer (Thorvaldsdóttir et al., 2013) and BLAST Ring Image Generator (Alikhan et al., 2011).

The PCGs and rRNA genes of the new mitogenome of R. viridisi and 36 monogenean mitogenomes available in GenBank (supplementary table S1) were used for a phylogenetic analysis. Two planarian species (Crenobia alpina and Obama sp.) were used as outgroups (Zhang et al., 2018a). The best partition scheme and the optimal models of molecular evolution as determined by PartitionFinder 2 (Lanfear et al., 2017) are presented in supplementary table S2. Based on the corrected Akaike information criterion, the optimal model of nucleotide evolution was GTR + G + I. A maximum likelihood (ML) tree was constructed in RAxML v. 8.1.12 (Stamatakis, 2014), with 1000 bootstrap replicates and joint branch length optimization. In addition, Bayesian inference (BI) was performed using MrBayes v. 3.2.2 (Ronquist et al., 2012) over ten million generations, sampling the Markov chain at a frequency of 100 generations and using the default settings. The best partition scheme and the optimal models of molecular evolution for the analysis of BI are presented in supplementary table S3. The CREx program (Bernt et al., 2007) was used to analyse rearrangement events in the mitogenomes and to make pairwise comparisons of gene orders of Diplectanidae utilizing the breakpoint dissimilarity measurement.

Results and discussion

Mitogenome characterization

The sequencing runs yielded totals of 15.14 million and 12.10 million reads for the DNA and cDNA, respectively. Of these, 8.63 million and 3.81 million reads, respectively, were retained after quality trimming and filtering. The raw reads were deposited into the NCBI's SRA database under accession numbers SRR13359482 and SRR11563382 (BioProject PRJNA689569).

The DNA reads were assembled into 21 contigs with lengths varying from 234 to 5892 bp and an average coverage of 1128 X. The cDNA reads were assembled into 1 contig with an average coverage of 3153 X (supplementary fig. S1). The assemblies of DNA and cDNA were aligned to correct errors in the frame of translation of the DNA assembly. Specifically, the nad4, nad1 and nad6 genes each presented a region with 1 bp, which we considered false. The correction of errors in these three genes seemed to be related to the assembly obtained by MEGAHIT rather than to coverage (coverage >1000 times). Furthermore, alignment helped to predict the trnC gen and the first 123 bases in the cox3 gene. The use of MEGAHIT and CAP3 to assemble the cDNA enabled identification of the trnS1 gene and two non-coding regions (NCRs) of 585 and 239 bp located between the trnG and trnC genes. For trnS1, cDNA coverage was greater than 300 bp, whereas for DNA, coverage was less than 7 bp. The greater coverage of certain regions in the mitochondrial genome from the cDNA reads, together with the use of MEGAHIT and TRINITY for their assembly, allowed us to obtain all the genes in the mitogenome of R. viridisi. It is important to note that the sequences used from the pool of individuals enabled identification of 120 nucleotides (0.86%) in the mitogenome in which at least 20% of the nucleotide coverage differed from the consensus. Whereas the ribosomal genes presented the fewest nucleotide variations (rrnL 3 bp, 0.32%; rrnS 1 bp, 0.14%), the NCRs had the most variations (33 bp, 3.9%). These bases might represent the most common variations within R. viridisi (supplementary fig. S2).

The nearly complete mitogenome of R. viridisi is 13,863 bp in size and includes 22 tRNA genes (including two for the amino acids serine and leucine), 12 PCGs (atp8 gene is absent) and two rRNA genes (table 1 and supplementary fig. S3). The total base composition is 30% A, 46.9% T, 7.6% C and 15.5% G (supplementary table S4). The annotated sequence was deposited in the NCBI GenBank (accession number MW565922). The majority of PCGs use start codons ATG or GTG, which are canonical for genetic code 9, and either TAA or TAG as the stop codon, which has been also reported for other monogeneans (e.g. Vanhove et al., 2018; Zhang et al., 2019). In L. longipenis and L. spari, TTG and ATT have also been reported as start codons for the nad2 gene (Zhang et al., 2018a). However, in R. viridisi, we only found ATG for nad2, and in nad4L, and TTG also appeared as start codon. The cox3 gene terminates with the abbreviated T codon, and thus was the only gene that does not use a canonical stop codon.

Table 1. Features of the mitochondrial genome of Rhabdosynochus viridisi.

A peak was found in the DNA coverage on the first NCR, which had an AT content of 78.5%, but without repeating regions. This might indicate that the length of this NCR was underestimated; hence we refer to the R. viridisi mitogenome reported here as a nearly complete.

Three NCRs were found. The first consists of 585 bp located between the trnG and trnS1 genes and has an AT content of 78.5%. The second NCR is 239 bp in size and is located between the trnS1 and trnC genes. It has an AT content of 91.2% and has a 182-bp region that consists of four tandem repeats (TRs) of between 39 and 48 bp each. The third NCR consists of 22 bp located between the trnC and cox3 genes, and has an AT content of 90.9%. The AT content of the second and third NCRs is much higher than in other parts of the mitogenome (supplementary table S4). The high number of TRs and the high AT content in these three NCRs, which are delimited by three tRNAs, posed a challenge to completing the mitogenome of R. viridisi. To overcome this problem, we combined the results obtained using three assemblers. It seems that a long NCR between nad5 and cox3 is common in diplectanids (Zhang et al., 2018a). The mitogenome of R. viridisi contains fewer TRs than those of L. longipenis and L. spari. Such differences in the number of TRs have also been observed among representatives of the Dactylogyridae family (Zhang et al., 2018a, b).

Phylogenetic and gene order analyses

The phylogenetic analysis included seven species from five families of the Polyopisthocotylea and 30 species from six families of the Monopisthocotylea. The topologies of the ML and BI trees were identical; nonetheless, the BI tree was better supported statistically (fig. 1). Each subclass formed a monophyletic group. Within the Monopisthocotylea, R. viridisi appeared clustered with diplectanids, as a sister group of P. yangjiangensis, with strong bootstrap support. The diplectanids formed a sister group with the clade formed by the Dactylogyridae and Capsalidae, which is consistent with mitogenome trees previously reported for monogeneans (Zhang et al., 2019).

Fig. 1. Phylogenetic analysis of the mitochondrial genomes of seven species from five families (Polystomatidae, Diplozoidae, Chauhaneidae, Mazocraeidae and Microcotylidae) of the Polyopisthocotylea and 30 species from six families (Gyrodactylidae, Tetraonchidae, Tetraonchoididae, Diplectanidae, Capsalidae and Dactylogyridae) of the Monopisthocotylea. Statistical support values of maximum likelihood and Bayesian analyses are shown below the nodes (ML/BI), except for nodes with maximum support. The scale bar corresponds to the estimated number of substitutions per site. Monogenean families are shown in different colours. Gene order is displayed to the right of the tree.

The gene order in the mitogenome of R. viridisi was different from the other members of the Diplectanidae family. Three transposition events of tRNA genes were found between R. viridisi and L. longipeni, and two events were found between R. viridisi and both P. yangjiangensis and L. spari (supplementary fig. S4). Gene order rearrangement has been previously observed in monogeneans. For instance, Ye et al. (2014) reported different arrangements in tRNA genes in monopisthocotylids, which was explained by the duplication–random loss model and recombination model together with a parsimonious scenario. In another example, Zhang et al. (2019) reported that the order of the two rRNAs in the mitogenomes of two species of Thaparocleidus was different with respect to all Neodermata. They suggested that evolution of mitogenomic gene order arrangements is discontinuous in monogeneans, as gene orders in a proportion of monogenean taxa are highly variable, whereas the remaining are conserved. The present study revealed that, although in a different order, there are three tRNA genes (trnG, trnS1 and trnC) located between nad5 and cox3 in the diplectanids L. longipenis and R. viridisi. This is in contrast with the gene order in the other monogenean families, in which only trnG is typically located between nad5 and cox3 (fig. 1). It is possible that this gene order is conserved among diplectanids.