The complete mitochondrial genome of Paragonimus ohirai (Paragonimidae: Trematoda: Platyhelminthes) and its comparison with P. westermani congeners and other trematodes

We present the complete mitochondrial genome of Paragonimus ohirai Miyazaki, 1939 and compare its features with those of previously reported mitochondrial genomes of the pathogenic lung-fluke, Paragonimus westermani, and other members of the genus. The circular mitochondrial DNA molecule of the single fully sequenced individual of P. ohirai was 14,818 bp in length, containing 12 protein-coding, two ribosomal RNA and 22 transfer RNA genes. As is common among trematodes, an atp8 gene was absent from the mitogenome of P. ohirai and the 5′ end of nad4 overlapped with the 3′ end of nad4L by 40 bp. Paragonimusohirai and four forms/strains of P. westermani from South Korea and India, exhibited remarkably different base compositions and hence codon usage in protein-coding genes. In the fully sequenced P. ohirai individual, the non-coding region started with two long identical repeats (292 bp each), separated by tRNAGlu. These were followed by an array of six short tandem repeats (STR), 117 bp each. Numbers of the short tandem repeats varied among P. ohirai individuals. A phylogenetic tree inferred from concatenated mitochondrial protein sequences of 50 strains encompassing 42 species of trematodes belonging to 14 families identified a monophyletic Paragonimidae in the class Trematoda. Characterization of additional mitogenomes in the genus Paragonimus will be useful for biomedical studies and development of molecular tools and mitochondrial markers for diagnostic, identification, hybridization and phylogenetic/epidemiological/evolutionary studies.


INTRODUCTION
Paragonimiasis, a neglected disease caused by lung flukes of the genus Paragonimus (Paragonimidae: Trematoda: Platyhelminthes), is a health threat in many tropical and subtropical regions. About seven species in the genus are common pathogens of humans, who usually acquire infection by eating undercooked or raw freshwater crabs or crayfish containing encysted metacercariae. Infection may also occur by eating raw meat of mammalian paratenic hosts (Yoshida et al., 2016).
Genomic information from mitogenomes of animals can be utilized for taxonomic identification, phylogenetic and population studies (Boore, 1999;Hardman & Hardman, 2006;Bernt et al., 2013;Wey-Fabrizius et al., 2013;Solà et al., 2015). A typical trematode mitogenome is a circular DNA molecule consisting of 12 protein-coding, 22 transfer RNA (tRNA) and two ribosomal RNA genes. Gene order tends to be highly conserved within the Order Plagiorchiida (that contains most trematodes, including paragonimids), but variable features, mainly repeats, do occur in the long non-coding regions of these genomes (Wey-Fabrizius et al., 2013;Le, Blair & McManus, 2002;Le et al., 2016). Currently, about 50 complete or near-complete mitogenomes are available from trematodes (Table S1). Many important families are not yet represented at all (Wey-Fabrizius et al., 2013), while other families, such as the Paragonimidae, are represented by data for only a few species. GenBank lists four mitogenomes for P. westermani of different ploidies and/or geographical origins. Coding regions are complete for all of these, but not all include the full non-coding region. Two of these mitogenomes are from eastern Asia, both from South Korea, a diploid (GenBank: AF540958, incomplete non-coding region) and a triploid (AF219379). Another two P. westermani mitogenomes available in Genbank, KM280646 (incomplete non-coding region) and KX943544, are from Northeast India.
Here, we provide genomic information from the complete mitochondrial genome for Paragonimus ohirai, determining the gene content, arrangement and composition, to compare features of this genome with those of the four available mitogenomes from P. westermani and of other members of the genus. We also assess the value of mitogenomes for resolution of the relationships between the trematode species.

Genomic DNA extraction
Following morphological identification, worms were preserved in 70% ethanol. Three adult worms derived from Kinosaki Prefecture and three from Nagoya Prefecture, were individually used for genomic extraction. The mitogenome of one specimen from Kinosaki was completely sequenced, remaining specimens were used to investigate variation in structure of the non-coding region.
Total genomic DNA was extracted from individual specimens using the GeneJET TM Genomic DNA Purification Kit (Thermo Fisher Scientific Inc., Waltham, MA, USA) according to the manufacturer's instructions. Genomic DNA was eluted in 50 µL of the elution buffer provided in the kit and stored at −20 • C until use. The concentration of DNA was estimated using a GBC UV/visible 911A spectrophotometer (GBC Scientific Equipment Pty. Ltd., Braeside, Australia). A working concentration of 50 ng/µL was prepared and 1 µL of this was used as template for short PCRs and 2 µL for long PCRs, in a 50 µL reaction volume.

Molecular identification of specimens
A molecular approach was also used to confirm the identity of our specimens as P. ohirai. Partial sequences of the cox1 gene were obtained from PCR products amplified using the primer pairs JB3F/JB4.5R (Table 1). A tree was inferred from an alignment of 224cox1 nucleotide sequences (309 bp) of 19 Paragonimus species available in GenBank and in previous publications. Phylogenetic reconstruction was performed using maximumlikelihood analysis (ML) with the Tamura-Nei model in the MEGA 7 package (Kumar, Stecher & Tamura, 2016). This model had the best Bayesian information criterion value, as determined using MEGA. Bootstrap support for each node was evaluated using 100 bootstrap resamplings. The tree inferred showed all P. ohirai sequences together in one group with those of the closely related Japanese P. iloktsuenensis (Fig. S1). Paragonimus iloktsuenensis and P. sadoensis have been regarded as synonyms of P. ohirai (Blair, Agatsuma & Watanobe, 1997).

Amplification of the mitochondrial genome
The complete mitogenome was sequenced from one specimen from Kinosaki Prefecture. Initially, platyhelminth-universal primers (URNLF/URNSR for the 16S rRNA-12S rRNA, TRECOBF/TRECOBR for cytochrome b (cob) and JB3F/JB4.5R for the cox1 region), were used for amplification and subsequent sequencing of the corresponding genes/regions (Le, Blair & McManus, 2002;Le et al., 2016). The sequence data obtained were used to design further specific P. ohirai primers (Table 1).
Some specific primers in nad5 or cob regions were used with the GLYF or GLYR platyhelminth-universal primers in long-PCRs. Long PCR was used to obtain mediumlength and long fragments (5-8 kb) of overlapping regionsand conventional PCRs were used to amplify the rest of the mitogenome. All long-PCRs were prepared using AccuPower R ProFi Taq PCR PreMixKit (Bioneer, Daejeon, South Korea) and short PCRs using Thermo Scientific DreamTaq PCR Master Kit (Thermo Fisher Scientific Inc., Waltham, MA, USA). All PCRs of 50 µL volume were performed in a MJ PTC-100 thermal cycler with initial denaturation at 95 • C for 5 min, followed by 35 cycles, each consisting of denaturation step for 30 s at 94 • C, annealing at 56 • C for 30 s, extension at 68 • C (in some cases, at 72 • C) for 3 to 10 min, depending on the lengths of the expected amplicons; and a final extension at 68 • C (or 72 • C for some reactions) for 7 or 10 min. A negative (no-DNA) control was included in some cases. The PCR products (5-10 µL of each) were examined on a 1% agarose gel, stained with ethidium bromide, and visualized under UV light (Wealtec, Sparks, NV, USA).
Notes. a Platyhelminth-universal primers used for inital amplification of the corresponding genes/regions to get sequences to design further primers.
Sequencing, assembling and annotation of the mitogenome PCR products were purified using a GeneJET Gel Extraction and DNA Cleanup Micro Kit (Thermo Scientific Inc.). In some cases, PCR products were cloned using a CloneJET PCR Cloning Kit (if <2 kb) or a TOPO R XL PCR Cloning Kit (Invitrogen Inc., Carlsbad, CA, USA) (if >2 kb). Plasmid DNA was extracted using a GeneJET Plasmid Miniprep Kit (Thermo Scientific, Inc.). A primer-walking approach to sequencing was used. PCR fragments and/or recombinant plasmid DNA were sequenced on automated sequencers using specific or M13 universal sequencing primers, respectively (Macrogen Inc., Seoul, South Korea). Both strands were completely sequenced and at least six sequences (three from each strand) were aligned to obtain the final sequence for characterization.
To explore size variation in the non-coding region (NR), we used forward primer GLYF (binding in tRNA Gly ) and reverse primer PO17R (located in the cox3 gene) (Table 1), spanning the full NR, for amplification and sequencing of this region from three additional P. ohirai individuals each from Kinosaki and from Nagoya.
Raw and edited sequences (edited using SeqEd v.1.3) were assembled and aligned using AssemblyLIGNv 1.9c and analyzed using the MacVector 8.2 package (Accelrys Inc.). Preliminary identity of a sequence or a region was assigned by searching the GenBank database (https://blast.ncbi.nlm.nih.gov/) using Basic-BLAST or two-sequence-BLAST, or by comparison with the corresponding platyhelminth sequences earlier obtained by us (Le, Blair & McManus, 2002;Le et al., 2016). Repeat sequences were detected in the long non-coding region (LNR) using the Tandem Repeats Finder v3.01 (Benson, 1999).
To confirm the length and to determine the composition and codon usage of each gene, comparisons of nucleotide and amino acid sequences of individual genes were done using MacVector 8.2 (Accelrys Inc.). The echinoderm/flatworm mitochondrial genetic code (translation table 9) was used in all programs for sequence analysis of mitogenomes of trematodes, including P. ohirai. Protein-coding genes were identified by the sequence similarities of translated amino acid sequences from their open reading frames to those already available in the GenBank database. Initiation codons other than ATG (specifically GTG) were considered. Codon usage and nucleotide composition were analyzed using MEGA 7 (Kumar, Stecher & Tamura, 2016). Codon usage for all the protein-coding genes was determined with the online program GENE INFINITY (Codon Usage: http://www.geneinfinity.org/sms/sms_codonusage.html). TreePuzzle v5.3 (Schmidt et al., 2002) was used to confirm (chi-square tests) that the amino-acid composition of each species did not significantly different from the average over the entire alignment (thus potentially biasing phylogenetic inference).
The transfer RNA genes (tRNA were identified using online software, primarily tRNAscan-SE v1.21 (Lowe & Eddy, 1997) with parameters specified for mitochondrial/chloroplast DNA. Three other separate software packages were also used online; ARWEN v 1.2 (Laslett & Canback, 2008); DOGMA (Wyman, Jansen & Boore, 2004) and MITOS Alpha version (Bernt et al., 2013). Final sequences and secondary structures were based on comparisons using all these programs. All transfer RNAs proposed by these programs were checked to confirm their typical ''clover-leaf'' structure or known variants of this. Any tRNAs not detected by these programs were found by inspection of the sequences, based on the alignment with other trematode tRNA sequences, to form tRNA structures.
Skew values (ranging from −1 to +1) were according to the formula defined by Perna & Kocher (1995) as: AT skew = (A − T)/(A + T) and GC skew = (G − C)/(G + C), where the letters stand for the usage of the corresponding nucleotides in the sequences (see also Solà et al., 2015).

Phylogenetic reconstruction
Concatenated aligned sequences of 12 mitochondrial protein-coding genes from 50 strains/species of 14 families of trematodes listed in Table S1 were used in a phylogenetic analysis. In addition to P. ohirai (Kinosaki strain, Japan), the four available P. westermani sequences (two from South Korea and two from India) were included. Also included were two mitogenomes of P. heterotremus (accession numbers MH059809 and KY952166) and one of P. kellicotti (MH322000). Protein-coding nucleotide sequences were conceptually translated and the amino acid sequences for each gene were aligned separately using MAFFT (available at https://www.ebi.ac.uk/Tools/msa/mafft/) (Katoh & Standley, 2013). Subsequently, the alignments of all 12 genes were concatenated (final alignment length 3,660 amino acids) for phylogenetic analysis. The inferred amino acid sequences of individual species ranged in length from 3,314 amino acids (Orientobilharzia turkestanicum (GenBank: HQ283100; Schistosomatidae)) to 3,405 amino acids (Metagonimus yokogawai (GenBank: KC330755; Heterophyidae)). Most sequences ranged from 3,351 to 3,361 amino acids. Prior to phylogenetic analysis, regions of questionable alignment quality were removed using GBlocks (Castresana, 2000). This reduced the number of positions to 2,509 amino-acid residues.
Phylogenetic reconstruction was performed in MrBayes 3.2 program (Ronquist et al., 2012). The amino-acid model prior was set to ''mixed''. This setting allows the program to infer the best fitting (from ten options) amino-acid substitution model. Two runs, each of four chains, were run for 1,000,000 generations, after which the standard deviation of split frequencies across runs was <0.01. Trees were sampled every 500 generations, and the first 25% discarded as burnin. The Jones, Taylor and Thornton model (Jones, Taylor & Thornton, 1992) was the only amino-acid model supported in the analysis (posterior probability 1.000 and standard deviation 0.000). The tree produced was a 50% majority-rule tree to which all compatible groupings were added. The tree was rooted by outgroup (i.e., members of the Schistosomatidae), a family known to be near the base of the trematode phylogeny in Olson et al. (2003).

Gene organization, content and genomic features
The entire mitochondrial genome sequence obtained from the single fully sequenced P. ohirai adult (designated as Pohi-Kinosaki-JP) was deposited in GenBank with the accession number KX765277. This mtDNA was a circular molecule of 14,818 bp (Fig. 1). Locations of genes and other features in this genome are given in Table 2. This is within the length range commonly reported in trematodes, including P. westermani (Paragonimidae), Echinochasmus japonicus (Echinochasmidae), Haplorchis taichui (Heterophyidae), Homalogaster paloniae (Gastrodiscidae), Dicrocoelium spp. (Dicrocoeliidae) (Le, Blair & McManus, 2002;Le et al., 2016;Lee et al., 2013;Liu et al., 2014a;Brabec et al., 2015;Ma et al., 2016a;Yang et al., 2016a), but slightly shorter than in some schistosomes such as Schistosoma spindale (Schistosomatidae) and longer than in opisthorchiids, fasciolids and paramphistomids such as Clonorchis sinensis, Opisthorchis viverrini, Fasciola hepatica,  Boore (1999) with slight modification by Le, Blair & McManus (2002). Genes are transcribed in a clockwise direction from one strand. The transfer RNA genes (tRNA) are designated by the single-letter code for the corresponding amino acid abbreviations, with the exception of those coding for Serine (tRNA Ser1(AGN ) and tRNA Ser2(UCN ) ) ( Table 2). The long non-coding region consists of two long repeats (LR1 and LR2) separated by tRNA Glu , and followed by six short tandem repeats (STR). A unique sequence is between the last STR6 and cox3.

Characteristics of protein-coding genes in Paragonimus species
In the P. ohirai mitogenome, ten protein-coding genes used ATG as the start codon and TAG as the stop codon and the remaining two genes were initiated with GTG and terminated with TAG and TAA codons (Table 2). Cox1 and cox2 genes were terminated by an incomplete TA-codon, which was presumably completed by the adjacent nucleotide, G (at the start of tRNA Thr andnad6, respectively). Translations in all the remaining protein-coding genes were consistent with translation Table 9. Unusual initiation codons, such as ATA or ATT, reported from some metazoans, and TTG or GTT, reported from nematodes, were not found to initiate genes in P. ohirai (Table 2). Base composition in protein-coding genes is very different between P. ohirai and the four mitogenomes of P. westermani (Table 3). In the former species, frequency of A is 17.00% (13.33-14.89% in P. westermani) and of T is 46.00% (38.21-40.50% in P. westermani). For both these bases, but particularly T, the base composition in P. ohirai is far closer to that typical of the other trematodes we included than that of P. westermani (Le, Blair & McManus, 2004) (Table 3). The G+C value in P. ohirai is lower than in the available P. westermani mitogenomes (Table 3), but closer to that in other trematodes (Le, Blair & McManus, 2004). In all the Paragonimus mitogenomes (as in trematodes generally), there is a strong AT skew: (−0.46 to −0.48). In P. ohirai, the preponderance of G over C is more marked (GC-skew 0.41) than in P. westermani (GC skew: 0.26 to 0.32) (Table 3), and closer to that seen in parasitic flatworms overall (Le, Blair & McManus, 2004). Base composition statistics for P. heterotremus and P. kellicotti, were intermediate between those for P. ohirai and P. westermani. All 64 codons in the mitochondrial genetic code table were used in the mtDNA protein-coding genes in P. ohirai (Table S1). There are 3,356 codons (excluding stop codons) in mitogenomes of P. ohirai and Indian P. westermani Type I (GenBank: KM280646), 3,355 codons in both Korean diploid and triploid P. westermani from East Asia and P. heterotremus and 3,354 in P. kellicotti. The remaining Indian P. westermani (from Arunachal Pradesh, KX943544) uses 3,368 codons (Table S1). An insertion of 45 nucleotides (15 amino acids) relative to other Paragonimus mitogenomes (and indeed all other trematodes) was found in the cytochrome b gene of this mitogenome. It remains to be determined whether this is a true feature unique to this Indian strain or whether it is due to mis-assembly. This region was omitted from phylogenetic analyses.
Codon usage indicated that codons containing T and/or A are far more common in P. ohirai and those containing C and/or G are less frequent than in P. westermani (Table  S1). For example, for phenylalanine (Phe), the TTT codon is much more frequently used in P. ohirai (325 codons/9.65%) than in diploid and triploid P. westermani (208/6.18-6.19%), P. westermani type I (255/7.57%) and P. westermani AP (253/7.49%). Similar differences in frequencies are noted for the ATT codon for isoleucine (Ile), TTA for leucine (Leu), TCT for serine (Ser), and TAT for tyrosine (Tyr) ( Table S1). These differences in base composition greatly influence the skew value calculation (Table 3).

Ribosomal and transfer RNA genes
The mitochondrial large and small ribosomal subunit RNA genes (16S rRNA and 12S rRNA, respectively) were located between the tRNA Thr andthe cox2 genes and were separated from each other by tRNA Cys , as in all platyhelminths reported to date (Le, Blair & McManus, 2002;Le et al., 2016;Liu et al., 2014a;Ma et al., 2016a;Ma et al., 2016b); 2017; (Shekhovtsov et al., 2010;Na et al., 2016;Biswal et al., 2014;Chen et al., 2016). The lengths of 16S rRNA and 12S rRNA, in the mitogenome of P. ohirai, were 974 bp and 736 bp, respectively ( Table 2). As in the protein-coding genes, usage of T was high relative to A and usage of C was low relative to G in ribosomal genes. Similarly, the AT-skew was comparable for this region of all Paragonimus mitogenomes (−0.23 to −0.26) but the GC-skew was slightly higher in P. ohirai (0.35) than in P. westermani (0.29-0.34) ( Table 3).
Twenty-two transfer RNAs were identified, ranging in size from 61 nucleotides (for tRNA Tyr ) to 72 (for tRNA Gly ): most are within the range 64-67 nucleotides ( Table 2). All anticodon usage is consistent with those described for other platyhelminth species. Twenty of the tRNAs can be folded into typical 'cloverleaf' secondary structures. The remaining two tRNAs for serine, tRNA Ser1(AGN ) (62 bp) and tRNA Ser2(UCN ) (64 bp), contain the full T C arm but not the dihydrouridine (DHU) arm (instead, it is replaced by a loop of ten nucleotides for the former and 13 for the latter) (Fig. 2). A DHU arm-lacking form for the former transfer RNA is usual in metazoans and is usual for the latter tRNA in platyhelminths and occurs in some other animal taxa, especially nematodes (Boore, 1999;Le, Blair & McManus, 2001a;Le, Blair & McManus, 2002;Le et al., 2016;Lee et al., 2013;Brabec et al., 2015;Ma et al., 2016a;Yang et al., 2016a;Littlewood et al., 2006;Wolstenholme, 1992). In P. ohirai, the tRNA Thr gene uses the last nucleotide G of the protein-coding cox1 gene termination codon; tRNA Pro overlaps tRNA Asn by 2 nucleotides; and tRNA Leu1(CUN ) overlaps tRNA Ser2(UCN ) by 3 nucleotides in their sequences (Table 2). tRNA Gly is located downstream of nad5 and tRNA Glu is located in the non-coding region, separating this region into two parts ( Table 2) Table S1). High A and T nucleotide use (17.00% and 46.00%) in P. ohirai is bolded. S.D. standard deviation nt nucleotide. a The overlap between nad4L and nad4 is counted twice in this calculation.

Polymorphism in non-coding regions feature long and short repeats
The major non-coding region (NR) in P. ohirai from Kinosaki is 1,465 bp in length and divided into two parts: the short NR (SNR, 328 nucleotides) bounded by tRNA Gly and tRNA Glu , and the long NR (LNR, 1,137 nucleotides) between tRNA Glu and cox3. There are two long identical repeats (LR1 and LR2, 292 bp each) and six short, identical, tandem repeats (STR1-6) each of 117 bp, in the fully sequenced individual (Pohi-Kinosaki-JP) ( Table 2; see GenBank: KX765277). LR1 is located within the SNR while the LR2 overlaps the 3 end of tRNA Glu by 29 bp at the start of the long non-coding region. STR1-6 immediately follow and a 157-nucleotide unique sequence region occurs between STR6 and cox3 (Table 2; Fig. S2).
We amplified and sequenced this non-coding region from additional specimens from Kinosaki and Nagoya (Fig. S3) (GenBank: KX765277; MF510407; MG214475-MG214478). In most cases, their NRs were identical in structure and sequence to that in the fully sequenced specimen from Kinosaki (Pohi-Kinosaki-JP). However, in one individual from Nagoya, this region was 2,602 bp long, 1,404 bp longer than that of the Pohi-Kinosaki-JP  sample. This individual from Nagoya (GenBank: MF510407) possessed 18 identical 117-bp STRs, 12 more than in the Pohi-Kinosaki-JP individual (Fig. S3).

T-A C-G T-A G-C A-T T-A A-T T CTTAC A A T G
Length variability in non-coding regions in most trematodes, including P. ohirai, is mainly due to different numbers of repeats, which are common features of this region. This polymorphism, seen also in P. westermani, Echinochasmus japonicus, Schistosoma spp. and other trematodes, contributes to a high level of mitochondrial inter-individual and inter-species variation (Blair et al., 2016;Bernt et al., 2013;Le, Blair & McManus, 2001a;Le, Blair & McManus, 2002;Le, Blair & McManus, 2004;Zarowiecki, Huyse & Littlewood, 2007;Le et al., 2016;Oey et al., 2019).
Some caution is required when comparing lengths and structures of non-coding regions from published sources. The long NRs are often not fully represented in mitogenomes in GenBank, despite statements to the contrary. This means that reported lengths of mitogenomes may be very different from actual lengths. When common next-generation sequencing (NGS) methods are used, the presence of multiple and long repeats causes assembly difficulties (Oey et al., 2019). PacBio long-reads were used to establish the structure of the NR in the mitogenome of an Indian P. westermani (Oey et al., 2019). This genome was 20.3 kb in length, with 6.3-6.9 kb of repetitive NR region containing three distinct types of repeats of lengths ranging between 229 and 406 bp. Oey et al., 2019 found evidence of length variation, that might occur at the intra-individual level. In BLAST searches, the closest match was with an Indian P. westermani Type 1 mitogenome (KM280646: 14,103 bp) that had been assembled from short SOLiD reads. Despite high sequence similarity elsewhere in these two genomes, there was only partial similarity with the repeat motifs.
A further two mitogenomes of Paragonimus species have been published recently. These are of P. kellicotti (13,927 bp by Wang et al. (2018) and P. heterotremus (13,927 bp by Qian et al. (2018). Both were assembled from short-read NGS data: the NR may therefore not be complete. In the NR of P. kellicotti, there are three copies (complete or nearly so) of a 111 bp region. In the published genome (MH322000) these are at positions 13, 345-13,355, 13,456-13,566 and 13,567-13,362. This repeat does not closely resemble those in other Paragonimus species. A putative tRNA Glu , not identified by Wang et al. (2018), starts at 13,289 and ends at 13,348. In the NR of P. heterotremus (MH059809) there are three complete or near-complete copies of a 37 bp region (13,539,13,614,13,717). Interestingly, this motif is a near-palindrome (with BLAST matches at almost the same locations in the reverse orientation). There are also two copies of a 52 bp region (13,614-13,665 and 13,762-13,813). No similarity was found with other Paragonimus species. tRNA Glu is from 13,343 to 13,408.

Phylogenetic analysis
Phylogenetic trees were constructed from concatenated amino-acid sequences from 50 taxa of trematodes (see Table S1). The Bayesian analysis recovered a monophyletic Paragonimidae with the cluster of P. westermani sequences appearing as sister to the remaining members of the genus. This matches previous phylogenetic schemes for Paragonimus (e.g., Nawa et al., 2014). The two Indian P. westermani mitogenomes, originating from the same region of India, despite clustering together, were more distinct from each other than were the two Korean (2n and 3n) mitogenomes (Fig. 3). Three distinct genotypes belonging to the P. westermani complex are known from northeastern India (Devi et al., 2010;Devi et al., 2013;Oey et al., 2019). The two genomes of P. heterotremus Figure 3 Phylogenetic tree based on concatenated amino-acid sequence data for the 12 mitochondrial proteins from 50 strains of 42 digenean trematode species (Table S1). Phylogenetic reconstruction was performed in MrBayes 3.2 program (Ronquist et al., 2012) with a final alignment of 2,509 amino-acid residues in length after removing poor-quality regions of the alignment. The tree produced (left-hand side of figure) was a 50% majority-rule tree to which all compatible groupings were added and rooted by outgroup (i.e., members of the Diplostomata -schistosomes and related taxa). Bayesian posterior support values for each node were 100% in every case except one (within the Opisthorchiidae, value of 57% indicated). Accession numbers are given at the end of each sequence label. Paragonimus ohirai in this study is marked by a star symbol. Full names of each species and of families where they belong are provided. The scale bar represents the number of substitutions per site. The tree on the right-hand side of the figure shows the relationships of the same families according to Olson et al. (2003).
are near-identical. Paragonimus ohirai is distinct from the other species, with P. kellicotti as a rather distant sister. The relationships depicted in the Bayesian tree matched reasonably closely those found by Olson et al. (2003) in a study using nuclear ribosomal sequences. Each included family was distinguished and grouped into the suborders proposed in Olson et al. (2003) (Fig. 3). For example, fasciolids and echinostomes are sister groups within the Echinostomata, as are paramphistomes and notocotylids (within the Pronocephalata), and also heterophyids and opisthorchiids (within the Opisthorchiata) (Fig. 3). The suborders mostly appear in the tree where expected from Olson et al. (2003), with two striking exceptions. The first is that the dicrocoeliids appear as rather basal in the tree, whereas they should belong close to Paragonimus within the derived suborder Xiphidiata (Olson et al., 2003). We can offer no convincing explanation for this discrepancy. The other exception is that, in our tree, the echinostomes and fasciolids (Echinostomata) are basal to the Pronocephalata, a situation reversed in Olson et al. (2003).

CONCLUSION
The present study provides the fully annotated mitogenome of Paragonimus ohirai (from an individual worm from Kinosaki, Japan) and a description of its genomic features in comparison with those of other paragonimids. There is an emphasis on comparisons with members of the P. westermani species complex. Although genomic organization is the same in P. westermani and P. ohirai, their mitogenomes differ remarkably in base composition, influencing codon usage and skew values. Despite this, phylogenetic analysis of the trematodes recovered a monophyletic family Paragonimidae. Fully characterized mitogenomes of additional paragonimid species will be useful for diagnostic, taxonomic, epidemiological, systematic, phylogenetic and population genetic studies.