Introduction
More than 50 species have been described in the genus Paragonimus (Trematoda: Paragonimidae) (Blair, Reference Blair, Bray, Gibson and Jones2008). As adults, Paragonimus species mature in the lungs of mammals. Freshwater snails and crustaceans act as first and second intermediate hosts, respectively. Eight species or subspecies of Paragonimus cause zoonotic diseases in humans, acquired by consumption of crustacean hosts containing metacercariae. These include Paragonimus westermani, Paragonimus heterotremus, Paragonimus skrjabini miyazakii, Paragonimus skrjabini skrjabini, Paragonimus kellicotti, Paragonimus mexicanus, Paragonimus uterobilateralis and Paragonimus africanus (see Chai and Jung, Reference Chai, Jung, Rose and Jiménez-Cisneros2018; Blair, Reference Blair, Toledo and Fried2019).
Some of these names (in particular P. westermani, P. heterotremus and P. skrjabini) refer to species complexes with multiple genotypes and strains, often with a geographical pattern (Blair, Reference Blair, Bray, Gibson and Jones2008, Reference Blair, Toledo and Fried2019; Blair et al., Reference Blair, Nawa, Mitreva and Doanh2016). The P. westermani complex has received the most attention. Its range extends from India and Sri Lanka eastwards to China, Japan and Korean Peninsula. Its northern limits are in the Russian Far East and its southern limits are in the Philippines, with a possible case in Papua New Guinea (Blair, Reference Blair, Toledo and Fried2019, Voronova et al., Reference Voronova, Vainutis, Tabakaeva, Sapotsky, Karaeka, Volkov, Galkina and Shchelkanov2022). Human cases, which are mostly confined to parts of eastern Asia and the Philippines, are typically of pulmonary paragonimiasis: pairs of worms form capsules in the lungs from which their eggs pass into bronchioles and are voided from the body via sputum or feces (Blair, Reference Blair, Toledo and Fried2019). Molecular studies have revealed great diversity across this geographic range. There is also diversity within single regions: 3 distinct genotypes (Types 1–3) occur in northeastern India alone, 2 in Sri Lanka and 2 in Thailand (Blair, Reference Blair, Toledo and Fried2019). One case of pulmonary paragonimiasis has been reported from northeast India, likely caused by the Indian Type 2 form (Singh et al., Reference Singh, Hiromu, Devi and Singh2015).
The P. skrjabini complex includes 2 nominal subspecies, P. skrjabini skrjabini Chen, 1959 in China, and P. skrjabini miyazakii Kamo, Nishida, Hatsushika and Tomimura, Reference Kamo, Nishida, Hatsushika and Tomimura1961, in Japan and parts of eastern China (Blair et al., Reference Blair, Chang, Chen, Cui, Wu, Agatsuma, Iwagami, Corlis, Fu and Zhan2005). Additional named taxa may also belong within this complex (Blair et al., Reference Blair, Chang, Chen, Cui, Wu, Agatsuma, Iwagami, Corlis, Fu and Zhan2005). The complex extends from northeast India into eastern and southern China, and Japan. Its members rarely mature in human lungs. Instead, worms wander into extra-pulmonary sites causing ectopic paragonimiasis with a wide range of signs and symptoms (Blair, Reference Blair, Toledo and Fried2019).
For accuracy in species identification, taxonomic classification and inferring phylogenetic relationships of Paragonimidae, markers from the mitogenomes are among the best tools (Boore, Reference Boore1999; Hardman and Hardman, Reference Hardman and Hardman2006; Bernt et al., Reference Bernt, Braband, Schierwater and Stadler2013a; Wey-Fabrizius et al., Reference Wey-Fabrizius, Podsiadlowski, Herlyn and Hankeln2013; Solà et al., Reference Solà, Álvarez-Presas, Frías-López, Littlewood, Rozas and Riutort2015). The mitogenome is a circular molecule with a length ranging from 13.5 kb to over 20 kb in trematodes (Le et al., Reference Le, Blair and McManus2002). As is typical for trematodes, the Paragonimus mitogenome consists of 12 protein-coding genes (PCGs) named cox1–3, nad1–6 and nad4L. atp6 and cytochrome b or cob or cytB, 2 mitoribosomal RNA genes (MRGs) as rrnL (or 16S) and rrnS (or 12S), 22 transfer RNA genes (tRNAs), and a variable non-coding region (NCR) that possesses multiple repeats (Biswal et al., Reference Biswal, Chatterjee, Bhattacharya and Tandon2014; Le et al., Reference Le, Nguyen, Nguyen, Doan, Agatsuma and Blair2019, Reference Le, Pham, Doan, Le, Saijuntha, Rajapakse and Lawton2020; Oey et al., Reference Oey, Zakrzewski, Narain, Devi, Agatsuma, Nawaratna, Gobert, Jones, Ragan, McManus and Krause2019b; Rosa et al., Reference Rosa, Choi, McNulty, Jung, Martin, Agatsuma, Sugiyama, Le, Doanh, Maleewong, Blair, Brindley, Fischer and Mitreva2020).
Although many complete mtDNA genome sequences have been reported from trematodes (Table S1 lists some examples, and see https://www.ncbi.nlm.nih.gov/genome/browse#!/organelles/Trematoda for a current list), relatively few are available for Paragonimus species (all are listed in Table S1). No mitogenome has been published for any member of the P. skrjabini complex. Next-generation sequencing (NGS) has become increasingly applied to generate mitogenomes for a range of species (Biswal et al., Reference Biswal, Chatterjee, Bhattacharya and Tandon2014; Oey et al., Reference Oey, Zakrzewski, Narain, Devi, Agatsuma, Nawaratna, Gobert, Jones, Ragan, McManus and Krause2019b). Single-molecule real-time (SMRT) sequencing from Pacific Biosystems (‘PacBio’) offers very long read lengths with high accuracy, overcoming errors caused by problematic genomic regions (Rhoads and Au, Reference Rhoads and Au2015).
The aim of this paper is to present the complete mtDNA sequence, which was obtained using NGS, from the Asian lung fluke P. s. miyazakii, representing the first for the P. skrjabini complex, and to use the mitogenomic markers for taxonomic and phylogenetic studies of the Paragonimidae. We also provide a detailed comparative analysis of the mitogenomic sequences and a phylogenetic analysis, emphasizing their utility as molecular markers for estimation of pairwise genetic distance and molecular evolutionary studies. The data are useful for further studies of the family Paragonimidae.
Materials and methods
Parasite samples and DNA extraction
The adult P. s. miyazakii samples used in this study were collected from experimentally infected albino rats (Rattus norvegicus) fed with metacercariae collected from the crab Geothelphusa dehaani hosts near Okuyanai, Kochi Prefecture, Japan (34° 39′ 42.01″ N, 133°56′ 6″ E). The metacercariae from the crabs and the adults obtained from the infected rats were identified by morphology and molecular analysis. The cox1 and ITS-2 sequences from samples in this study were 99–100% identical with these markers from all P. s. miyazakii samples collected previously in different geographical localities in Japan, including Kochi, Miyazaki, Tokushima, Shizuoka and Yamaguchi (Kamo et al., Reference Kamo, Nishida, Hatsushika and Tomimura1961; Agatsuma et al., Reference Agatsuma, Ketudat, Thaithong, Shibahara, Sugiyama, Habe, Terasaki and Kawashima1992; Blair et al., Reference Blair, Agatsuma, Watanobe, Okamoto and Ito1997, Reference Blair, Chang, Chen, Cui, Wu, Agatsuma, Iwagami, Corlis, Fu and Zhan2005; Ohari et al., Reference Ohari, Suzuki, Shibahara and Itagaki2019; Urabe and Marcaida, Reference Urabe and Marcaida2022).
Total genomic DNA was extracted from individual adult worms using the GeneJET™ Genomic DNA Purification Kit (Thermo Fisher Scientific Inc., MA, USA) according to the manufacturer's instructions. The genomic DNA was eluted in 50 μL of the elution buffer and kept at −20°C until use. The DNA content was quantified using a NanoDrop® ND-1000 UV-Vis Spectrophotometer. A working concentration (50 ng μL−1) was prepared and an amount of 2 μL was used in each LPCR in a 50 μL reaction volume.
Targeted enrichment of the mitogenome by long polymerase chain reaction
Trematode-PCR primers based on the alignment of all mtDNA of trematodes available and specific primers binding to nad5, cox3 and cytochrome b (cob) were designed for amplifying the whole mtDNA of P. s. miyazakii in 5 overlapping fragments (Tables S2 and S3). Among these, 2 amplicons (obtained by primer pairs PMI12F/PMI15R and GLYF/PMI20R, respectively) were used for comparative validation of the length of the NCR as well as these amplicons were used for the NGS (Table S3). Long-range PCRs (LPCR) were performed in a 50 μL volume reaction in a MJ PTC-100 Thermal Cycler. Each reaction contained 25 μL 2× LongAmp Master Mix (New England Biolabs, Ipswich, MA, USA), 2 μL each primer (10 pmol μL−1), 2 μL of template DNA, and 19 μL DEPC-water). Long-PCRs were conducted with an initial denaturation at 94°C for 1 min, followed by 30 cycles, each consisting of a denaturation step for 30 s at 94°C, an annealing/extension step at 50°C for 30 s, extension at 65°C for 8 min, and a final extension at 65°C for 10 min. The PCR products (10 μL of each) were examined on a 1% agarose gel, stained with ethidium bromide and visualized under UV light (Wealtec, Sparks, NV, USA)
All dsDNA products were purified using the GeneJET PCR Purification Kit (Thermo Fisher Scientific, USA) and the amplicon length was verified via 1.5% agarose gel electrophoresis. Three amplicons from the coding mtDNA and 2 from the NCR (Table S3) were pooled for NGS. The complete mitogenome of P. s. miyazakii was sequenced using the PacBio SEQUEL system (https://www.pacb.com/) with a targeted long-read sequencing approach at the PacBio facility at the Institute of Biotechnology (Hanoi, Vietnam).
Library preparation
The dsDNA products of each sample from 5 overlapping amplicons were pooled into one tube and purified with AMPure XP beads (Pacific Biosciences, Menlo Park, CA, USA). Input dsDNA was quantified using the Qubit fluorometer 3.0 and Qubit dsDNA HS Assay reagents (Thermo Fisher Scientific, Waltham, MA, USA). SMRTbell Libraries were prepared using the Express Template Prep Kit 2.0 with multiplexing amplicons protocol with low DNA input (100 ng) (Pacific Biosciences, Menlo Park, CA, USA) for sequencing on the PacBio SEQUEL system according to the manufacturer's instructions. The SMRTbell templates were purified once with 1.2 volumes of AMPure PB beads and the size and amount of the library were checked again using the Bioanalyzer 2100 system (Agilent, CA, USA) and the Qubit fluorometer 3.0 with Qubit™ dsDNA HS Assay reagents, respectively. The libraries of all amplicons were then pooled before long-read sequencing.
Sequencing and de novo assembly
The pooled library was bound to polymerase using the Sequel Binding and the Internal Control Kit 3.0 (Pacific Biosciences, Menlo Park, CA, USA) and purified using Ampure PB beads. The DNA Control Complex 3.0 and the Internal Control Kit 3.0 from Sequel Binding and Internal Ctrl Kit 3.0 were used to control the sequencing procedure. The final library was loaded onto Sample Plate (Pacific Biosciences, Menlo Park, CA, USA). The run design was created by the Sample Setup software included in the SMRTLink portal v5.1 version 9.0 with an insert size of 1200 base pairs (bp). The sequencing signals were processed, evaluated and converted into raw data by the Primary Analysis Computer server. All data were automatically transferred to the Secondary Analysis Server system via the intranet. High-quality sequence data were proofread and generated by PacBio's circular consensus sequencing (CCS), then de novo assembled using Canu software v2.0 (Koren et al., Reference Koren, Walenz, Berlin, Miller, Bergman and Phillippy2017), and the quality of the assembly was checked by using Quast software v5.0.2 (Gurevich et al., Reference Gurevich, Saveliev, Vyahhi and Tesler2013).
Mitogenome annotation and characterization
Twelve PCGs were identified by comparative alignment with the available mitogenomes of the Paragonimus strains and species in GenBank or recently reported (Biswal et al., Reference Biswal, Chatterjee, Bhattacharya and Tandon2014; Qian et al., Reference Qian, Zhou, Li, Wang, Miao and Hu2018; Wang et al., Reference Wang, Wang, Xu, Li, Qu, Song, Tang and Lin2018; Oey et al., Reference Oey, Zakrzewski, Narain, Devi, Agatsuma, Nawaratna, Gobert, Jones, Ragan, McManus and Krause2019b; Le et al., Reference Le, Nguyen, Nguyen, Doan, Agatsuma and Blair2019) (Table S1). ATG/GTG as start and TAA/TAG as stop codons were used to define individual PCG gene boundaries. The echinoderm and flatworm mitochondrial genetic code (Translation Table 9 in GenBank) was used for the translation of the PCGs. The nucleotide composition for PCGs, MRGs and the mtDNA coding region (5′-cox3-nad5-3′, designated as mtDNA*) and codon usage of the PCGs were analysed using MEGA X (https://www.megasoftware.net/) (Kumar et al., Reference Kumar, Stecher, Li, Knyaz and Tamura2018). Codon usage for all concatenated 12 PCGs from each strain/species was determined with the online program GENE INFINITY (Codon Usage: http://www.geneinfinity.org/sms/sms_codonusage.html). Genetic distance was estimated based on the pairwise amino-acid comparisons of concatenated 12 PCGs among 11 strains of 5 species (or species complexes) of the family Paragonimidae and was determined using GENEDOC 2.7 (https://softdeluxe.com/GeneDoc-180568/download/) for alignment and MEGA X for percentage estimation.
Transfer RNA genes were identified using the tRNAscan-SE 1.2.1 program (www.genetics.wustl.edu/eddy/tRNAscan-SE/) (Lowe and Chan, Reference Lowe and Chan2016), ARWEN (http://mbio-serv2.mbioekol.lu.se/ARWEN/) (Laslett and Canback, Reference Laslett and Canback2008), and the MITOS Alpha version (http://mitos.bioinf.uni-leipzig.de/index.py) (Bernt et al., Reference Bernt, Donath, Jühling, Externbrink, Florentz, Fritzsch, Pütz, Middendorf and Stadler2013b). 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, and any tRNAs not detected by these programs were found by inspection of the sequences based on the alignment with other trematode tRNA sequences known to form tRNA structures. Repeat sequences were detected in the NCR using the Tandem Repeats Finder v3.01 (Benson, Reference Benson1999). The circular map and gene abbreviations on the map were generated using GenomeVx (Conant and Wolfe, Reference Conant and Wolfe2008) with the manual drawing of the tandem repeat units.
Skew values (ranging from − 1 to + 1) were determined by calculating the percentage of AT and GC nucleotide usage using Perna and Kocher's (Reference Perna and Kocher1995) formula: AT skew = (A − T)/(A + T), and GC skew = (G − C)/(G + C), where the letters represent the absolute usage of the corresponding nucleotides in the sequences [see also Solà et al. (Reference Solà, Álvarez-Presas, Frías-López, Littlewood, Rozas and Riutort2015)].
Phylogenetic reconstruction
Concatenated aligned sequences of 12 PCGs from 36 strains/species of 5 families of trematodes listed in Table S1 were used in a phylogenetic analysis. Schistosoma haematobium (Schistosomatidae) was used as an outgroup. In addition to P. s. miyazakii (OkuST1 strain, Japan), the 10 available Paragonimus mtDNA sequences were included. Among these were 2 strains of P. heterotremus (LC, Vietnam and GX, China), 1 strain of P. ohirai (Kino strain, Japan) and 6 representatives of the P. westermani complex (2 from South Korea, 1 from China and 3 from India). All of the Indian sequences appear to belong to Type 1 of Devi et al. (Reference Devi, Narain, Mahanta, Nirmolia, Blair, Saikia and Agatsuma2013), based on comparisons of available cox1 sequences. One strain of P. kellicotti (Ozark, Missouri, the United States) was also included. The paper reporting this sequence has been retracted by the journal (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7748660/) but the sequence data appear to match expectations: the reader can choose whether to interpret these data with caution. The concatenated protein-coding nucleotide sequences were imported into GENEDOC 2.7 (available at: https://softdeluxe.com/GeneDoc-180568/download/) and translated using Translation Table 9 in GenBank. The inferred amino-acid sequences ranged in length from 3,354 residues (Paragonimus kellicotti, Pkel-Ozark-US; GenBank: MH322000; and P. s. miyazakii; GenBank: ON782295; Paragonimidae) to 3403 residues (Metagonimus yokogawai, GenBank: KC330755; Heterophyidae). The alignment was performed by GENEDOC 2.7, confirmed by the MAFFT 7.122 program (available at https://www.ebi.ac.uk/Tools/msa/mafft/) (Katoh and Standley, Reference Katoh and Standley2013), and lightly trimmed by Gblocks 0.91b for the best-quality final sequence block to avoid excessive ambiguity (http://molevol.cmima.csic.es/castresana/Gblocks_server.html) (Talavera and Castresana, Reference Talavera and Castresana2007). The phylogenetic tree was reconstructed using the maximum-likelihood method in MEGA X with 1000 bootstrap resamplings with a final alignment of 3,146 amino-acid residues after removing gap sites (Kumar et al., Reference Kumar, Stecher, Li, Knyaz and Tamura2018). The substitution model with the best score, according to the Bayesian information criterion, was the Jones, Taylor and Thornton + F + G + I model, with residue frequencies estimated from the data ( + F), rate variation along the length of the alignment ( + G) and allowing for a proportion of invariant sites ( + I).
Results
Paragonimus skrjabini miyazakii mitogenomic organization and features
Paragonimus skrjabini miyazakii, strain OkuST1 from Kochi, Japan, possessed the complete mitogenome of 17 591 bp in length (GenBank accession no. ON782295) (Table 1; Fig. 1). The NCR in the mitogenome of P. s. miyazakii, which is about 4.3 kb, was successfully obtained using the long-reading PacBio system. The circular mtDNA molecule is comprised of 12 PCGs (cox1–3, cob, nad1–6, nad4L, atp6), 2 MRGs (16S or rrnL and 12S or rrnS) and 22 tRNAs or trn, and an NCR, rich in long and short tandem repeats (LRUs and SRUs). The length of the mitogenome was similar to that of members in the Echinostomatidae, e.g. Artyfechinostomum malayanum (17 175 bp, strain EMI3, Thailand) (Pham et al., Reference Pham, Saijuntha, Lawton and Le2022), Echinostoma revolutum (17 030 bp, strain MSD15, Thailand) (Le et al., Reference Le, Pham, Doan, Le, Saijuntha, Rajapakse and Lawton2020), shorter than Paragonimus westermani (20 273 bp, strain IND2009, India, GenBank: CM017921) (Oey et al., Reference Oey, Zakrzewski, Narain, Devi, Agatsuma, Nawaratna, Gobert, Jones, Ragan, McManus and Krause2019b) and Echinostoma paraensei (20,298 bp, GenBank: KT008005) but longer than some other Paragonimus species (Table S1). The linearized map of the P. s. miyazakii mtDNA is 5′-cox3-H-cob-nad4L-nad4-QFM-atp6-nad2-VAD-nad1-NPIK-nad3-S1W-cox1-T-rrnL-C-rrnS-cox2-nad6-YL1S2L2R-nad5-G-NCR[LRU1(E)–8.6#]-[SRU1–6.8#]-3′. The NCR was located between tRNAGlu and cox3. The NCR of P. s. miyazakii mtDNA was bounded by tRNAGlu (trnE) and the cox3 gene, was relatively long (4351 bp) and was divided into 2 subregions. The first subregion contained 8.6 LRUs (8 perfect LRUs of 384 bp each−1 and a partial one of 241 bp) and the second contained 6.8 SRUs (6 perfect SRUs of 118 bp each−1 and a partial one of 93 bp) (Table 1; Fig. 1 and Fig. S1). The 8 LRUs were identical and each LRU (384 bp) contained a 66 bp complete sequence of tRNAGlu (trnE1–E8) in between each pair of LRUs (Fig. S1; Table 1).
A schematic drawing of a circular map of the mitogenome of Paragonimus skrjabini miyazakii (GenBank: ON782295). The protein-coding (PCGs) and mitoribosomal large and small subunit genes (MRGs) are abbreviated as in our previous publications (Le et al., Reference Le, Nguyen, Nguyen, Doan, Agatsuma and Blair2019, Reference Le, Pham, Doan, Le, Saijuntha, Rajapakse and Lawton2020). The transfer RNA genes (tRNAs) are marked with 3-letter amino acid abbreviations (see Table 1). The NCR contains tandem repeat units, and is located between tRNAGly and cox3 with 8.6 LRUs (long, LRU1–8.6) including tRNAGlu in each (8 tRNAGlu) and 6.8 SRUs (short repeat units, SRU1–6.8).
Locations of genes and other features in the mitochondrial genome of Paragonimus skrjabini miyazakii (17 591 bp; GenBank: ON782295)
bp, basepair; aa, amino acid; Int. seq., intergenic sequence ( + , number of nucleotides before start of the following gene; −, number of nucleotides overlapping with the following gene); LRU, long tandem repeat unit, which contains a tRNAGlu; SRU, short tandem repeat unit; Sterisk (*), indicating tRNAs lacking DHU-arm.
Ribosomal and transfer RNA genes
The 16S or rrnL gene is 976 bp and the 12S or rrnS gene is 738 bp. These were located in the region between tRNAThr (trnT) and cox2 and separated from one another by tRNACys (trnC) (Table 1). Twenty-two tRNAs ranged from 61 (trnS1) to 71 nucleotides (trnN) in length. Except for 2 tRNAs that specify serine, tRNASer1(AGN) (trnS1) and tRNASer2(UCN) (trnS2), which possessed a special form, missing the DHU-arm (dihydrouridine-arm), the remaining 20 tRNAs had the common ‘cloverleaf’ secondary structures with the complete 4 arms (Table 1; Fig. S2).
PCGs and codon usage
ATG was used as the start codon in 9 PCGs, and the other 3 genes used GTG. Three PCGs (cox1, cox2 and nad6) used TAA as the stop codon; the other 9 genes used TAG. An overlap of 40 bp exists between nad4L and nad4. There is no atp8 gene (Table 2).
Base composition and skew/skewness value for AT and GC of the PCGs, ribosomal genes (MRGs), and the coding region overall (5′-cox3-nad5-3′, designated as mtDNA*) of the mitogenomes of Paragonimus skrjabini miyazakii and other members of Paragonimidae
PCGs, protein-coding gene; MRGs, mitoribosomal gene; mtDNA*, mitochondrial coding nucleotide sequence (from 5′ cox3 to 3′ nad5).
The codon usage values for 12 PCGs in the mitogenomes of P. s. miyazakii and the other 10 strains/species of the Paragonimidae are summarized in Table S4. The most frequently used codons are TTT (for phenylalanine), followed by TTG (for leucine) and GTT (for valine) in all paragonimids. In P. s. miyazakii, TTT/Phe (10.13%, n = 341), TTG/Leu (8.02%, n = 270), and GTT/Val (5.67%, n = 191) were the most commonly used. The least frequently used codons were CGC/Arg (n = 5) in P. s. miyazakii and P. kellicotti, and CCA/Pro (n = 5) in P. westermani (2 strains from Korea and 1 from China). Only 1 CGC/Arg codon was used in P. ohirai. The least-used codons are CGA/Arg in P. heterotremus (1 strain from Vietnam and 1 from China) and CAA/Gln in P. westermani (2 strains from Korea and 1 from China) (Table S4).
Base composition and skew values
Table 2 presents the base composition of A, T, G and C and skewness values of AT and GC content for PCGs, MRGs and mtDNA* (the coding region from 5′ cox3 to 3′ nad5) of 11 strains/species of Paragonimidae. Paragonimus skrjabini miyazakii used 15.83% A and 44.76% T (A + T = 60.59%), 27.70% G and 11.71% C (G + C = 39.41%) for 12 PCGs; and 17.33% A and 42.90% T (A + T = 60.23%), 27.72% G and 12.05% C (G + C = 39.77%) for its coding mtDNA* region (13 240 bp). The AT-skew was highly negative ( − 0.477/PCGs and −0.425/mtDNA*), indicating that T was more frequently used than A. The GC-skew was highly positive (0.406/PCGs and 0.394/mtDNA*), indicating that greater numbers of C were used than G. The base composition and skewness values for PCGs and mtDNA* seemed to be different between P. s. miyazakii and other paragonimids, for example, P. ohirai (Le et al., Reference Le, Nguyen, Nguyen, Doan, Agatsuma and Blair2019), and also slightly different between Korean/Chinese and Indian P. westermani strains. The AT skew values were the most negative (−0.491 to −0.492/PCGs and −0.441/mtDNA*) in P. heterotremus (Vietnamese and Chinese strains) and (−0.492/PCGs and −0.442/mtDNA*) in P. kellicotti (Table 2). The pattern of the base usage for PCGs and mtDNAs* is T > G > A > C in P. s. miyazakii, and is similar in P. heterotremus, P. ohirai and P. kellicotti, whereas it is T > G > C > A in P. westermani (Korean and Chinese strains, the East-Asian P. westermani subgroup). However, in the Indian P. westermani, the T > G>A > C is for the mtDNAs*, but it is T > G > C > A for the PCGs.
For MRGs, P. s. miyazakii used 22.11% A and 37.11% T (A + T = 59.22%), and 27.77% G and 13.01% C (G + C = 40.78%) with an AT skew of −0.253 and a GC skew of 0.362. This pattern was similar to that of P. heterotremus and P. kellicotti but different from that of P. westermani. The East-Asian (Korean and Chinese) strains used A + T less frequently than the Indian P. westermani (51.58–51.73% vs 54.85–55.42%, respectively) (Table 2). For MRGs, the T > G > C > A pattern of nucleotide usage is recorded in all strains/species of Paragonimidae.
Taxonomic and phylogenetic relationships
Genetic distances between species of Paragonimidae
Pairwise genetic distances (P-distances) were estimated based on the analysis of concatenated amino acid sequences inferred from 12 PCGs each among 11 strains/5 species in the family Paragonimidae (Table 3). Paragonimus skrjabini miyazakii had the lowest genetic distance (14.71–13.99%) to P. heterotremus and, P. kellicotti (16.21%), and the highest (24.34–24.72%) to Indian P. westermani [strains Pwes-AP-IN (KX943544) and Pwes-(typeI)-IN (KM280646)]. There was a low intra-specific genetic distance within strains in the P. heterotremus complex (1.25%) and between the Korean and Chinese P. westermani (0.23–1.21%), but interestingly, a relatively high distance between members of the P. westermani complex within India (3.71–7.35%). The P. westermani complex seemed to have 2 evolutionary sublineages with different genetic distances among the geographical strains (Table 3).
Estimation of pairwise genetic distances (%) among Paragonimus skrjabini miyazakii and strains/species of the family Paragonimidae inferred from the concatenated amino acid sequences of 12 mitochondrial PCGs
Information for strain/species is given in Table S1. Pairwise genetic distance (%) among P. s. miyazakii and species of Paragonimidae showing the lowest (numbers 2–3, and 11) and the highest (numbers 8–9) are bolded. The intra-specific distances between strains within the P. heterotremus complex and within the P. westermani complex (the East-Asian groups and the Indian groups) are boxed.
Phylogenetic relationships
To resolve the phylogenetic relationships between P. s. miyazakii and the other congeners within Paragonimidae and between this family and other families in suborders Opisthorchiata, Troglotremata and Echinostomata, a maximum likelihood phylogenetic analysis was performed on the alignment of the concatenated amino acid sequences of 12 PCGs for 36 strains/26 species, including P. s. miyazakii (listed in Table S1). Schistosoma haematobium (family Schistosomatidae) was used as an outgroup species. In the phylogenetic tree (Fig. 2), the family Paragonimidae representing the suborder Troglotremata was rendered as a monophyletic clade and was placed as a sister clade to Opisthorchiata, comprising 2 families, Heterophyidae and Opisthorchiidae. The topology clearly showed 2 subclades with absolute nodal support (100%). One contained all the sequences of the P. westermani complex (Korean, Chinese and Indian strains), and the other was formed by sequences from the remaining Paragonimus species (P. s. miyazakii, P. heterotremus, P. ohirai and P. kellicotti) (Fig. 2). The P. westermani clade itself is comprised of 2 subclades, one of the East-Asian and the other of the Indian Type 1 strains, indicating the complicated phylogenetic relationships within this complex.
A maximum likelihood phylogenetic tree showing the position of Paragonimus skrjabini miyazakii (diamond symbol) based on the analysis of the complete concatenated amino acids of 12 PCGs of 36 trematode species/strains. Schistosoma haematobium (Platyhelminthes: Schistosomatidae) is included as an outgroup. Nodal support values evaluated using 1000 bootstrap resamplings are shown on each branch. The basal node of the Paragonimidae is shown by an arrow. Accession numbers are given for each species/strain and country of origin (in brackets) (where available) at the end of each sequence. The scale bar represents the number of substitutions per site.
Discussion
The complete mitogenome of Paragonimus skrjabini miyazakii was obtained by long-read NGS using the PacBio system and showed gene order and content similar to those of other trematodes (Biswal et al., Reference Biswal, Chatterjee, Bhattacharya and Tandon2014; Ma et al., Reference Ma, He, Liu, Leontovyč, Kašný and Zhu2016; Oey et al., Reference Oey, Zakrzewski, Narain, Devi, Agatsuma, Nawaratna, Gobert, Jones, Ragan, McManus and Krause2019b; Le et al., Reference Le, Nguyen, Nguyen, Doan, Agatsuma and Blair2019, Reference Le, Pham, Doan, Le, Saijuntha, Rajapakse and Lawton2020; Li et al., Reference Li, Qiu, Zeng, Diao, Chang, Gao, Zhang and Wang2019; Rajapakse et al., Reference Rajapakse, Pham, Karunathilake, Lawton and Le2020; Rosa et al., Reference Rosa, Choi, McNulty, Jung, Martin, Agatsuma, Sugiyama, Le, Doanh, Maleewong, Blair, Brindley, Fischer and Mitreva2020; Pham et al., Reference Pham, Saijuntha, Lawton and Le2022). The NCRs of trematode mtDNAs typically contain long and short repeat units but obtaining complete sequences of these is difficult using conventional Sanger and/or second-generation sequencing methods (Kinkar et al., Reference Kinkar, Korhonen, Cai, Gauci, Lightowlers, Saarma, Jenkins, Li, Li, Young and Gasser2019, Reference Kinkar, Young, Sohn, Stroehlein, Korhonen and Gasser2020; Rosa et al., Reference Rosa, Choi, McNulty, Jung, Martin, Agatsuma, Sugiyama, Le, Doanh, Maleewong, Blair, Brindley, Fischer and Mitreva2020). However, the PacBio SMRT sequencing can yield the long, complicated, and repetitive region of the NCR, which was 4.3 kb for P. s. miyazakii in this study; 6.9 kb for P. westermani (Oey et al., Reference Oey, Zakrzewski, Narain, Devi, Agatsuma, Nawaratna, Gobert, Jones, Ragan, McManus and Krause2019b); 4.4 kb for Echinococcus granulosus (Kinkar et al., Reference Kinkar, Korhonen, Cai, Gauci, Lightowlers, Saarma, Jenkins, Li, Li, Young and Gasser2019); and 4 kb for Schistosoma bovis (Oey et al., Reference Oey, Zakrzewski, Gravermann, Young, Korhonen, Gobert, Nawaratna, Hasan, Martínez, You, Lavin, Jones, Ragan, Stoye, Oleaga, Emery, Webster, Rollinson, Gasser, McManus and Krause2019a). The P. s. miyazakii mitogenome, at 17 591 bp in length, is one of the longest Paragonimus mtDNAs obtained to-date (Oey et al., Reference Oey, Zakrzewski, Narain, Devi, Agatsuma, Nawaratna, Gobert, Jones, Ragan, McManus and Krause2019b; Rosa et al., Reference Rosa, Choi, McNulty, Jung, Martin, Agatsuma, Sugiyama, Le, Doanh, Maleewong, Blair, Brindley, Fischer and Mitreva2020).
A notable finding was the presence of multiple LRUs and SRUs in the NCR and the repeated 66 bp sequence of tRNAGlu (trnE1–E8) included between each pair of LRUs. Multiple long-repetitive sequences tandemly arranged in the NCRs are also found in other trematodes, as summarized in Le et al. (Reference Le, Pham, Doan, Le, Saijuntha, Rajapakse and Lawton2020). The presence of long repeats in the mitogenome has been reported not only in P. s. miyazakii but also in the mtDNAs of Paragonimus westermani, Clonorchis sinensis, Schistosoma bovis and the cestode Echinococcus granulosus G1 (Oey et al., Reference Oey, Zakrzewski, Gravermann, Young, Korhonen, Gobert, Nawaratna, Hasan, Martínez, You, Lavin, Jones, Ragan, Stoye, Oleaga, Emery, Webster, Rollinson, Gasser, McManus and Krause2019a, Reference Oey, Zakrzewski, Narain, Devi, Agatsuma, Nawaratna, Gobert, Jones, Ragan, McManus and Krause2019b; Kinkar et al., Reference Kinkar, Korhonen, Cai, Gauci, Lightowlers, Saarma, Jenkins, Li, Li, Young and Gasser2019, Reference Kinkar, Young, Sohn, Stroehlein, Korhonen and Gasser2020). The multiple copies of trnE in the LRUs raise the question as to which copy might be functional, or whether all the tRNAs are active for transferring Glu. In the cestode Echinococcus granulosus G1, tRNAGly was found in each of the 10 repeat units of the NCR (Kinkar et al., Reference Kinkar, Korhonen, Cai, Gauci, Lightowlers, Saarma, Jenkins, Li, Li, Young and Gasser2019).
The mito-ribosomal gene order and location were conserved in paragonimids, as seen in P. heterotremus, P. ohirai, P. westermani and P. kellicotti and many other trematodes reported to date (Biswal et al., Reference Biswal, Chatterjee, Bhattacharya and Tandon2014; Qian et al., Reference Qian, Zhou, Li, Wang, Miao and Hu2018; Wang et al., Reference Wang, Wang, Xu, Li, Qu, Song, Tang and Lin2018; Oey et al., Reference Oey, Zakrzewski, Narain, Devi, Agatsuma, Nawaratna, Gobert, Jones, Ragan, McManus and Krause2019b; Le et al., Reference Le, Nguyen, Nguyen, Doan, Agatsuma and Blair2019, Reference Le, Pham, Doan, Le, Saijuntha, Rajapakse and Lawton2020; Rosa et al., Reference Rosa, Choi, McNulty, Jung, Martin, Agatsuma, Sugiyama, Le, Doanh, Maleewong, Blair, Brindley, Fischer and Mitreva2020). The absence of a DHU-arm, which was replaced by a DHU-replacement loop in both tRNAs specifying serine (trnS1 and trnS2) in the mitogenome of P. s. miyazakii, is commonly found in other Paragonimus species, i.e. P. ohirai (Le et al., Reference Le, Nguyen, Nguyen, Doan, Agatsuma and Blair2019) and in many other trematodes in the Echinostomatidae, Fasciolidae and Opisthorchiidae, i.e. Clonorchis sinensis, Opisthorchis viverrini and O. felineus (Shekhovtsov et al., Reference Shekhovtsov, Katokhin, Kolchanov and Mordvinov2010; Cai et al., Reference Cai, Liu, Song, Wu, Zou, Yan, Yuan, Lin and Zhu2012), Fasciola hepatica, F. gigantica, Fasciola jacksoni and Fascioloides magna (Le et al., Reference Le, Blair and McManus2001, Reference Le, Blair and McManus2002; Liu et al., Reference Liu, Gasser, Young, Song, Ai and Zhu2014; Ma et al., Reference Ma, He, Liu, Leontovyč, Kašný and Zhu2016; Rajapakse et al., Reference Rajapakse, Pham, Karunathilake, Lawton and Le2020), Echinostoma miyagawai, E. revolutum and Echinostoma/Artyfechinostomum malayanum (Li et al., Reference Li, Qiu, Zeng, Diao, Chang, Gao, Zhang and Wang2019; Le et al., Reference Le, Pham, Doan, Le, Saijuntha, Rajapakse and Lawton2020; Pham et al., Reference Pham, Saijuntha, Lawton and Le2022).
The base composition of A, T, G and C, as well as skewness values of AT and GC content for PCGs, MRGs and mtDNA*, indicated the use of T was more frequent than A, and G more than C. These yield the highly negative values for AT-skew and highly positive GC-skew in the mtDNA of all Paragonimus species. However, the nucleotide usage was not similar among the Asian/American and the Indian paragonimids and was even different between the East-Asian and the Indian P. westermani mitogenomes. Whereas the T > G > C > A pattern for nucleotide use was found in the mtDNA* of the Asian/American Paragonimus spp. (P. s. miyazakii, P. heterotremus, P. ohirai and P. kellicotti) and the East Asian P. westermani (Korean and Chinese strains), the Indian Type 1 P. westermani exhibited a different usage pattern (T > G > A > C).
The most frequently used codons (TTT for Phe, TTG for Leu and GTT for Val) appeared similar in all the 11 mitogenomes of the 5 Paragonimus species, while the least frequently used codons were different in each species and between the East-Asian and Indian P. westermani subgroups. The least frequently used codons in P. heterotremus were (CGC/Arg, n = 3), in East-Asian P. westermani (CAA/Gln, n = 3). In P. s. miyazakii and P. kellicotti, 5 codons were equally least common and in P. ohirai, CGC/Arg was only used once.
We inferred a phylogenetic tree from an alignment of concatenated, translated PCGs. This revealed the Paragonimidae to be monophyletic, given the limited mitogenome sampling possible. The separation of Paragonimus species into 2 main groups (the P. westermani complex in one and all remaining species in the other) reflects the pattern seen using other sequence data (e.g. Blair, Reference Blair2022). Similarly, the P. westermani group itself is subdivided, with the 3 Indian Type 1 mitogenomes being distinct from those from eastern Asia (and rather different from each other). There are too few complete mitogenomes from the complex for us to be able to add to previous findings based on shorter sequence markers. But it is clear that relationships within the P. westermani complex defy easy explanation (Doanh et al., Reference Doanh, Tu, Bui, Loan, Nonaka, Horii, Blair and Nawa2016; Blair, Reference Blair, Toledo and Fried2019). For example, within northeastern India, 3 distinct clusters have been recognized: Type 1, which the mitogenomes used here represent, Type 2, which is phylogenetically close to eastern Asian strains infecting humans and Type 3 which might represent Paragonimus siamensis (Devi et al., Reference Devi, Narain, Mahanta, Nirmolia, Blair, Saikia and Agatsuma2013). No mitogenomes are available as yet for either Type 2 or Type 3.
In our phylogenetic tree, P. s miyazakii appeared closest to P. heterotremus. This reflects previous findings based on molecular data, which place these 2 species in a clade with P. macrorchis and P. vietnamensis (e.g. see Blair, Reference Blair2022). Unfortunately, there is no complete mtDNA sequence for P. s. skrjabini and for other related species, such as P. hokuoensis, P. macrorchis or P. vietnamensis, which could assist in the resolution of the phylogenetic position of P. s. miyazakii and P. s. skrjabini in the Paragonimidae.
Conclusion
The present study provides the complete mitogenome of Paragonimus skrjabini miyazakii (from 1 individual worm, from Kochi, Japan), obtained by PacBio long-read sequencing of long-PCR amplicons. Annotations included a description of its genomic features, including base composition, genetic distance and codon usage in comparison with those of other Paragonimus congeners. Mitogenomic organization and features were the same as in other Paragonimus species. The long NCR was completely sequenced for the P. s. miyazakii mtDNA and featured long repeat units. Two major groups were apparent in a phylogenetic tree: members of the P. westermani complex and the remaining members of the genus for which mitogenomes were available (P. s. miyazakii, P. heterotremus, P. ohirai and P. kellicotti). Despite the monophyly of the family Paragonimidae being revealed by the phylogenetic analysis, the taxonomic and phylogenetic relationships of P. s. miyazakii and P. s. skrjabini need to be confirmed. Fully characterized mitogenomes of additional paragonimid species, particularly those in the major Paragonimus complexes, are necessary and will be useful for diagnostic, taxonomic, epidemiological, systematic, phylogenetic and population genetic studies.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0031182022001184
Acknowledgements
We express our thanks to colleagues and technicians for their contribution to our laboratory work. We are grateful to Mr. Le Tung Lam and Ms. Le Thi Dzung for their technical help with NGS by PacBio at the Institute of Biotechnology, Hanoi, Vietnam.
Author's contributions
Thanh Hoa Le conceived the idea of the study, analysed the final data, prepared figures and tables and wrote the paper; Khue Thi Nguyen, Linh Thi Khanh Pham, Huong Thi Thanh Doan conducted laboratory works and sequence analyses. Takeshi Agatsuma collected, experimentally infected, molecularly identified and provided specimens. David Blair reviewed the drafts. Takeshi Agatsuma, David Blair and Thanh Hoa Le completed and approved the paper. All authors read and approved the final paper.
Financial support
This research is funded by the Vietnam National Foundation for Science and Technology Development (NAFOSTED) under grant number 108.02–2020.07) (PI: Thanh Hoa Le).
Conflict of interest
None.
Ethical standards
None.
