Characterization of the complete mitochondrial genomes of Diplodiscus japonicus and Diplodiscus mehari (Trematoda: Diplodiscidae): Comparison with the members of the superfamily Paramphistomoidea and phylogenetic implication.

Diplodiscus japonicus and Diplodiscus mehari (Trematoda: Diplodiscidae) are two important parasites in wood frogs, which have large infection rates and essential importance of ecology, economy and society. In this study, the complete mitochondrial (mt) genomes of D. japonicus and D. mehari were sequenced, then compared with other related trematodes in the superfamily Paramphistomoidea. The complete circular mt sequence of D. japonicus and D. mehari were 14,210 bp and 14,179 bp in length, respectively. Both mt genomes comprised 36 functional subunits, consisting of 12 protein-coding genes (PCGs), two ribosomal RNA (rRNA) genes, 22 transfer RNA (tRNA) genes, and one non-coding region. The mt genes of D. japonicus and D. mehari were transcribed in the same direction, and the gene arrangements were identical to those of Paramphistomoidea trematodes. In the 12 PCGs, GTG was the most common initiation codon, whereas TAG was the most common termination codon. All tRNAs had a typical cloverleaf structure except tRNA Ser1. A comparison with related Paramphistomoidea trematode mt genomes suggested that the cox1 gene of D. mehari was the longest in these trematodes. Phylogenetic analyses revealed that Paramphistomoidea trematodes formed a monophyletic branch, Paramphistomidae and Gastrothylacidae were more closely related than Diplodiscidae. And the further analysis with Pronocephalata branch found that the flukes parasitic in amphibians (frogs) formed one group, and the flukes from ruminants (cattle, sheep, ect) formed another group. Our study demonstrated the importance of sequencing mt genomes of D. japonicus and D. mehari, which will provide significant molecular resources for further studies of Paramphistomoidea taxonomy, population genetics and systematics.

Introduction

Wood frogs, as amphibians, play a vital role in the connect of water and land environment in ecosystems, and they are one of the most essential components of the biological community (Morrison and Hero, 2003). Nowdays, there are more than 30 species of wood frogs of genus Rana in the world and 17 species in China (Lu and Li, 2005; Hu et al., 2021). However, many factors (such as habitat alteration, climate change, chemical contaminants and diseases) would affect the survival of the frogs, and various infectious and parasitic diseases were known to be the main factors that cause the decline in global amphibian population (Jiménez and Sommer, 2016). Parasitic diseases are one of the most serious factors, they not only affected the health of frogs, but also induced human disease (Eamsobhana, 2014). And as for the variations of frog behavior, the parasite of Green Frogs (genus Pelophylax) and Brown Frogs (genus Rana) were different (Okulewicz et al., 2014). Two kinds of Borwn Frogs, Rana amurensis and Rana dybowskii are mainly distributed in northeastern China. Therefore, to know the parasite species of these frogs and the characteristics of parasites is an urgent issue. In our previous study, from October 2019 to December 2021, the parasitic infection of R. amurensis and R. dybowskii in Heilongjiang Province, northeastern China was investigated. The results showed that 11 species of trematodes, five species nematodes and two species tapeworms were detected. The infection rate of trematodes was 37.75%, D. mehari (22.38%) was the highest and D. japonicus (14.19%) was follow in two kinds of wood frogs (Wang et al. unpublished data).

In the past, morphological features played a key role in identifying and differentiating parasites. However, morphological approach existed extensive limitations in the comparison analyses of related species and in-depth analyses, and the emergence of molecular parasitology has resolved many of these problems. For example, Tkach et al. (2003) used 28S rDNA as genetic marker to examine the relationships for superfamily Microphalloidea, the result showed that Prosthogonimus ovatus and Schistogonimus rarus had higher homology than P. ovatus and Prosthogonimus cuneatus, and phylogenetic analysis also showed that P. ovatus was closer to S. rarus than to its congener P. cuneatus, which proved that Schistogonimus is not a valid genus and S. rarus should belong to Prosthogonimus (Tkach et al., 2003). In addition, because of mt genome have high copy number, maternal inheritance, rapid evolution and a virtual lack of recombination, and it usually used as the important molecular marker to study the taxonomy, population genetics and systematics of parasites (Biswal et al., 2014; Le et al., 2020; Gao et al., 2021). For example, Wang et al. (2011) reported that Orientobilharzia turkestanicum belongs to Schistosoma based on mtDNA genome data (Wang et al., 2011). For nematode, the hypothesis that Cylicostephanus minutus represents a complex species has also been supported with mtDNA data (Gao et al., 2020).

The family Diplodiscidae is a small group of paramphistomoids, found predominantly in amphibians (Jones, 2005). Diplodiscidae species are important parasite in the intestine of frogs and common prevalent in China, which can affect the growth of frogs in severe cases. Gao and Zhang (2014) reported that the infection rate of Diplodiscus sp. was 25% in Rana grahami Boulenger in Yunnan Province, southwestern China (Gao and Zhang, 2014). Diplodiscus nigromaculati was the most prevalent species in Pelophylax nigromaculatus in Shanghai, eastern China, and the infection rate was as high as 64.44%, and Diplodiscus sp. also reached 37.78% (Men et al., 2016). Our research also found that D. japonicus and D. mehari were common trematodes of R. amurensis and R. dybowskii in Heilongjiang Province, northeastern China. Although there are currently 11 Diplodiscus species have been reported in the world, however, the studies mainly focused on morphology and epidemiology (Cichy and Żbikowska, 2016; Men et al., 2016; Besprozvannykh et al., 2018). The research on molecular biology is quite limited, especially in mtDNA, although the complete mt genome of D. nigromaculati was sequenced (GenBank: MW698822.1), no relevant papers have not been published so far. Thus, the objective of our study was to amplify the complete D. japonicus and D. mehari mt genome sequences, analyze the features of these mt genomes, and reconstruct the phylogenetic tree of D. japonicus and D. mehari to analyze the relationships with other related Paramphistomoidea trematodes.

Materials and methods

Sample collection, specimen identification, and DNA extraction

Diplodiscus japonicus and D. mehari used in this study were collected from the intestine of R. dybowskii naturally infected in Heilongjiang Province. This study was approved by the Animal Ethics Committee of Heilongjiang Bayi Agricultural University (Approval No. HBAU2018-007). The trematodes were thoroughly washed in physiological saline and identified using morphological characteristics according to previous standards (Besprozvannykh et al., 2018), then fixed in 75% (v/v) ethanol at −20 °C until use. According to the manufacturer's instructions, total genomic DNA was extracted from individual worms using TIANamp Genomic DNA Kit (TIANGEN Biotech, Beijing, China). Molecular identification was performed by amplifying ITS sequences, and the primers were NC5 (5′- GTA GGT GAA CCT GCG GAA CGA TCA TT -3′) and NC2 (5′- TTA GTT TCT TTT CCT CCG CT -3′) (Gasser et al., 2008).

Mitochondrial genome amplification, sequencing, assembly, and annotation

To amplify the mt genomes of D. japonicus and D. mehari, primers were designed by Primer 5.0 software according to the mt sequences of other related trematodes in NCBI. The primers of D. japonicus and D. mehari are shown in Table S1. PCRs were performed in 25 μL volumes that contained 1 μL DNA sample, 12.5 μL of 2 × Tks Gflex PCR Buffer (Mg2+, dNTP plus) (Takara, Dalian, China), 0.5 μL of each primer (10 pmol/μL) (synthesized by Qingke Biotechnology Co. Ltd., Harbin, China), 10 μL of ddH2O, and 0.5 μL of Tks Gflex DNA Polymerase (1.25 U/μL) (Takara, Dalian, China). The reactions were performed in a thermocycler under the following conditions: 95 °C for 1 min (initial denaturation), followed by 35 cycles of 98 °C for 10 s (denaturation), 50–64 °C for 30 s (annealing), and 68 °C for 30 s–2 min, and a final extension of 72 °C for 7 min. Each amplicon was examined in a 1.0% agarose gel, stained with ethidium bromide, and photographed upon transillumination. The DL2000 and DL5000 markers (Takara, Dalian, China) were used to estimate the sizes of the mtDNA amplicons. Representative PCR products were sent to Tianyihuiyuan Biotechnology Co. Ltd (Beijing, China) for sequencing. The mtDNA sequences of D. japonicus and D. mehari were assembled using DNAStar (v. 5.0) (Burland, 2000). The sequencing results were split with DNAStar, then compared with the mitochondrial gene sequences of related trematodes in GenBank to determine the gene boundaries.

Comparative analyses of Diplodiscus japonicus and Diplodiscus mehari with other related Paramphistomoidea species

Twelve protein-coding gene sequences were translated into corresponding amino acid sequences using MEGA X (Kumar et al., 2018). The assembled complete mt genomes were annotated with the MITOS Webserver (http://mitos2.bioinf.uni-leipzig.de/index.py) (Bernt et al., 2013). An online open reading frame finder (https://www.ncbi.nlm.nih.gov/orffinder/) was used to analyze and translate the PCGs. Contents of A + T and G + C were calculated using DNAStar (v. 5.0) (Burland, 2000). AT-skew and GC-skew values were calculated using the equations AT-skew = (A−T)/(A + T) and GC-skew = (G−C)/(G + C) in both coding genes and non-coding regions (NCRs) (Perna and Kocher, 1995). Transfer RNA (tRNA) genes were inferred using tRNAscan-SE (http://lowelab.ucsc.edu/tRNAscan-SE/) (Lowe and Chan, 2016) and ARWEN (http://130.235.46.10/ARWEN/) (Laslett and Canbäck, 2008). Tandem repeats were detected in NCRs using the Tandem Repeat Finder (http://tandem.bu.edu/trf/trf.html) (Benson, 1999) and DNAMAN (v. 9.0) (Lynnon Biosoft Company, Foster City, CA, USA). The relative synonymous codon usage of the 12 PCGs was calculated using the Sequence Manipulation Suite (http://www.detaibio.com/sms2/codon_usage.html) (Stothard, 2000). The relative synonymous codon usage bias (RSCU) for PCGs was determined with the trematode mt genetic code in CodonW. (do Nascimento et al., 2021).

The nucleotide and amino acid sequence differences were calculated using MEGA X and MegAlign (v. 5.01) (Burland, 2000; Kumar et al., 2018). Comparative analyses were conducted among D. nigromaculati, Fischoederius elongatus, Fischoederius cobboldi, Gastrothylax crumenifer, Calicophoron microbothrioides, Orthocoelium streptocoelium, Paramphistomum cervi, Paramphistomum leydeni and Explanatum explanatum in superfamily Paramphistomoidea.

Mutation rate (nonsynonymous/synonymous, dN/dS) among the 12 PCGs of the three Diplodiscus species mitogenomes were calculated with DnaSP (v.5) (Librado and Rozas, 2009). DnaSP was also used to conduct sliding window analysis implementing window size of 300 bp, and the nucleotide divergences between 12 PCGs and two ribosomal RNAs of 11 Paramphistomoidea trematodes were characterized with a step size of 10 bp.

Phylogenetic analysis

The concatenated amino acid sequences of the complete D. japonicus and D. mehari mt genomes were aligned with the corresponding amino acid sequences of 32 trematodes from seven suborders, the sequences available in GenBank are shown in Table S2. Gyrodactylus salaris (NC008815) was selected as an outgroup. MAFFT (v. 7.471) (Katoh and Standley, 2013) was used for the alignment and Gblocks server (http://molevol.cmima.csic.es/castresana/Gblocks_server.html) was used to truncate and exclude ambiguously aligned regions from the multiple amino acid sequence alignment using the “less stringent” selection option (Castresana, 2000).

Phylogenetic trees were reconstructed from the alignment of concatenated amino acid sequences of mtDNA of tremetodes using maximum likelihood (ML) and Bayesian inference (BI) methods. PhyML (v. 3.0) was used to reconstruct the ML tree with the JTT + I + G + F model, and 100 bootstrap replicates were performed (Guindon and Gascuel, 2003). MrBayes 3.1 was used to conduct bayesian analysis. Four independent Markov chain runs were performed for 1,000,000 metropolis-coupled MCMC generations, sampling a tree every 100 generations. The first 25% (2500) trees were omitted as burn-in, and the remaining trees were used to calculate Bayesian posterior probabilities (Ronquist and Huelsenbeck, 2003). The phylograms were drawn using FigTree v1.4.4 (http://tree. bio.ed.ac.uk/software/figtree/).

Results and discussion

Identification of parasites

Under the microscope, the morphological characteristics of the collected parasites were observed. The body of the two trematodes was conical or trapezoidal, non spined, with narrow anterior end and wide posterior end, and the ventral sucker and oral sucker were located at both ends of the trematode body (Fig. S1). These characteristics were agreed with all qualitative and morphometric characteristics of the genus Diplodiscus (Besprozvannykh et al., 2018). To further confirm, we performed ITS sequencing. ITS sequences of D. japonicus (ON248512) and D. mehari (ON248512) obtained in present study were 943 bp and 947 bp in size, which had 100% nucleotide identities with those flukes form Russia deposited in GenBank (D. japonicus, KX506855; D. mehari, KX506857). So, the two species obtained in this study were identified as D. japonicus and D. mehari.

The general features of the Diplodiscus japonicus and Diplodiscus mehari mt genome

The complete circular mt genomes of D. japonicus and D. mehari were 14,210 bp and 14,179 bp in length, respectively, and have been deposited in GenBank (accession: OL961442 and OL961441). The mt genomes of D. japonicus and D. mehari included 36 genes: 12 PCGs (cox1–3, nad1–6, nad4L, atp6, and cytb), 22 tRNA genes, and two rRNA genes (16S rRNA, or rrnL; and 12S rRNA, or rrnS); however, they lacked an atp8 gene (Fig. 1 and Table 1), which were similar to other trematodes (Yan et al., 2013; Ma et al., 2015). The arrangements of mt genes were closely spaced in D. japonicus and D. mehari. Some genes were separated by 1–61 bp intervals, some genes had overlap between 11 and 40 bp, but another genes were adjacent without any spacer and overlap (Table 1). There were 40 bp overlap between nad4L and nad4 in both D. japonicus and D. mehari, which was consistent with most other trematodes, such as Tracheophilus cymbius, Paragonimus westermani, Clonorchis sinensis, and F. hepatica (Shekhovtsov et al., 2010; Biswal et al., 2014; Li et al., 2019; Liu et al., 2014), but was longer than those of Schistosoma spindale (28 bp) and Schistosoma mekongi (37 bp) (Le et al., 2000; Littlewood et al., 2006). Gene transcription and replication were all carried out in the same direction, different from many other metazoan mt genomes (Saccone et al., 2002). The nucleotide composition of D. japonicus and D. mehari mt genomes was biased toward A + T, with an overall A + T content of 62.14% and 59.26% (Fig. 2A). The overall A + T contents of D. japonicus and D. mehari were similar to most digenean trematodes, including C. sinensis (59.8%) and Haplorchis taichui (59.2%) (Shekhovtsov et al., 2010; Lee et al., 2013). However, the D. japonicus and D. mehari mt genome A + T contents were higher than in Brachycladium goliath (55.6%) and P. westermani (51.68%) (Biswal et al., 2014; Briscoe et al., 2016).

Fig. 1

Diplodiscus japonicus and Diplodiscus mehari mitochondrial genomes arrangement. All 22 tRNA genes are designated by the one-letter code with numbers differentiating each of the two tRNAs leucine and serine. All genes are coded by the same DNA strand and are transcribed clockwise. NCR refers to the non-coding region.

Table 1

Mitochondrial genome organization of Diplodiscus japonicus and Diplodiscus mehari.

Genes	Position	Length (bp)	Initiation codon	Stop codon	Anticodon	Intergenic nucleotide
	D. j/D.m	D.j/D.m	D.j/D.m	D.j/D.m	D.j/D.m	D.j/D.m
cox3	1–645/1–645	645/645	ATG/ATG	TAA/TAG		11/13
trnH	657–718/659–714	62/56			GTG/GTG	6/9
Cytb	725–1837/724–1839	1113/1116	ATG/ATG	TAA/TAG		11/9
nad4L	1849–2112/1849–2112	264/264	ATG/GTG	TAG/TAG		−40/–40
nad4	2073–3359/2073–3359	1287/1287	ATG/GTG	TAA/TAA		8/8
trnQ	3367–3429/3368–3431	63/64			TTG/TTG	1/2
trnF	3431–3495/3434–3501	65/68			GAA/GAA	24/17
trnM	3520–3585/3519–3582	66/64			CAT/CAT	3/3
atp6	3589–4104/3586–4101	516/516	ATG/ATG	TAG/TAA		9/18
nad2	4114–4983/4120–4989	870/870	GTGATG/	TAG/TAG		14/34
trnV	4998–5066/5024–5093	69/70			TAC/TAC	33/17
trnA	5100–5162/5111–5174	63/64			TGC/TGC	10/6
trnD	5173–5238/5181–5245	66/65			GTC/GTC	1/1
nad1	5240–6142/5247–6149	903/903	GTG/GTG	TAG/TAG		11/10
trnN	6154–6222/6160–6229	69/70			GTT/GTT	25/29
trnP	6248–6312/6259–6323	65/65			TGG/TGG	0/0
trnI	6313–6373/6324–6384	61/61			GAT/GAT	20/25
trnK	6394–6451/6410–6479	58/70			CTT/CTT	11/0
nad3	6463–6816/6480–6833	354/354	GTG/GTG	TAA/TAA		13/1
trnS1	6830–6889/6835–6894	60/60			AGC/AGC	18/20
trnW	6908–6978/6915–6983	71/69			TCA/TCA	3/–21
cox1	6982–8529/6963–8534	1548/1571	GTG/GTG	TAG/TAA		16/1
trnT	8546–8619/8546–8618	74/73			TGT/TGT	2/2
rrnL	8622–9629/8621–9625	1008/1005				−12/–11
trnC	9618–9683/9615–9680	66/66			GCA/GCA	3/3
rrnS	9685–10451/9682–10448	767/767				−13/–11
cox2	10439–11023/10438–11022	585/585	ATG/ATG	TAA/TAG		6/12
nad6	11030–11485/11035–11490	456/456	GTG/GTG	TAG/TAG		14/15
trnY	11500–11564/11506–11569	65/64			GTA/GTA	25/14
trnL1	11590–11660/11584–11653	71/70			CTA/TAG	3/1
trnS2	11664–11730/11655–11717	67/63			TGA/TGA	61/28
trnL2	11792–11858/11746–11811	67/66			TAA/TAA	13/14
trnR	11872–11938/11826–11889	67/64			TCG/TCG	2/3
nad5	11941–13503/11893–13455	1563/1563	ATG/ATG	TAG/TAA		15/21
trnG	13519–13586/1347713545	68/69			TCC/TCC	19/16
trnE	13606–13673/13562–13629	68/68			TTC/TTC	0/0
NCR	13674–14210/13630–14179	537/550				0/0

Note: D. j, Diplodiscus japonicus, D. h, Diplodiscus mehari.

Fig. 2

A + T content and nucleotide skew of genes, individual elements, and the complete mitogenome of 11 Paramphistomoidea trematodes.

Transfer RNAs are an indispensable part of mitochondria, which play an essential role in amino acid transport. The total lengths of 22 tRNAs were 1,451 bp and 1,449 bp in the D. japonicus and D. mehari mt genomes, respectively. Individual gene lengths vary from 56 bp to 74 bp (Table 1). The lengths of D. japonicus and D. mehari tRNA were similar with to those of P. cervi (59–72 bp), P. leydeni (59–73 bp), G. crumenifer (59–72 bp), and O. streptocoelium (59–72 bp) (Yan et al., 2013; Ma et al., 2015; Yang et al., 2016; Zhao et al., 2017), which suggested that tRNA gene size has changed quite slowly during trematode evolutionary divergence. All tRNA sequences could be folded into the typical cloverleaf structure, except trnS1, which lacked the dihydrouridine (DHU) arm in both Diplodiscus mitogenomes (Fig. S2). The structure of trnS1 was consistent with those in the previous report on trematode, such as Tamerlania zarudnyi (Suleman et al., 2021). The first site of anticodon sequences of all tRNAs was T or G, except for trnM and trnK in two studied mitogenomes (Fig. S2). Nucleotide substitutions observed among the tRNAs of D. japonicus and D. mehari mitogenomes are mainly restricted to the TΨC and DHU loops, whereas the anticodon loop is highly conserved (Fig. S2).

The rrnL genes in the mtDNAs of D. japonicus and D. mehari were located between trnT and trnC, and the rrnS genes were located between trnC and cox2, which were the same to other published Paramphistomoidea trematodes. But it was different from Schistosoma turkestanicum, the rrnL gene was located between trnT and rrnS, the rrnS gene was located between rrnL and cox2 (Wang et al., 2011). The rrnL and rrnS of D. japonicus were 1008 bp and 767 bp in length, and rrnL and rrnS genes of D. mehari were 1005 bp and 767 bp in length, respectively. The rrnL gene lengths of Diplodiscus species were longer than that of other published Paramphistomoidea trematode mt genomes, e.g. F. elongatus (995 bp), and P. cervi (986 bp) (Han et al., 2020; Yan et al., 2013) and were shorter than T. zarudnyi (1244bp) and Tanaisia sp. (1231 bp) (Suleman et al., 2021). The rrnS gene lengths of D. japonicus and D. mehari were also longer than that of other published Paramphistomoidea trematode mt genomes, e.g. D. nigromaculati (730 bp), F. elongatus (751 bp), and P. cervi (749 bp) (Han et al., 2020; Yan et al., 2013) and were shorter than Tanaisia sp. (777bp) (Suleman et al., 2021).

In the mt genomes of D. japonicus and D. mehari, the NCRs were located between trnE and cox3, which were 541 bp and 550 bp in length, respectively (Table 1). There were only one NCR in both D. japonicus and D. mehari, which were similar to D. nigromaculati (MW698822). However, this was different from other trematodes in superfamily Paramphistomoidea, such as F. cobboldi, O. streptocoelium, G. crumenifer, P. cervi, P. leydeni and F. elongatus, which had two NCRs (Yan et al., 2013; Ma et al., 2015; Yang et al., 2016; Zhao et al., 2017; Han et al., 2020). In Paramphistomoidea trematodes, the longest overall NCR length was G. crumenifer (813 bp) and the shortest overall NCR length was O. streptocoelium (310 bp). Meanwhile, the longest complete genomes is G. crumenifer (14,801 bp) and the shortest complete genomes is O. streptocoelium (13,800 bp). It also indirectly proved that the length of the complete mt genome was related to the length of NCRs, which was consistent with previous study in T. zarudnyi (Suleman et al., 2021). The A + T contents of the D. japonicus and D. mehari NCRs were 64.99% and 58.73%, respectively. In addition, two tandem repeats with 43 bp (TATACTATATATATATATATATATGTATATATAT-TATATATA) in size were found in the NCR of D. japonicus, but similar sequences were not found in D. mehari. A previous study showed that repetitive sequences contain much genetic information, and they could regulate gene replication, transcription, and translation (Jurka et al., 2007). And the frequencies of recombination and base mutations of repetitive sequences were much higher than non-repetitive regions, which ensure the stability of coding genes during replication, recombination, and mutation (Araujo and Arias, 2019). Since the significant genetic variability of NCR regions in trematodes (Littlewood et al., 2006), the research on the function of NCRs in detail is very important.

Comparative analysis

The complete mt genome lengths in superfamily Paramphistomoidea trematodes were different, D. japonicus and D. mehari were similar to those of F. elongatus (14,120 bp) and F. cobboldi (14,256 bp), but shorter than G. crumenifer (14,801 bp) and longer than E. explanatum (13,968 bp) (Table 2). The gene transcription direction of the complete mt genomes of D. japonicus and D. mehari were the same as the all Paramphistomoidea trematodes mt genomes. The rank order of the 12 PCGs of D. japonicus by length was as follows: nad5 > cox1> nad4 > cytb > nad1 > nad2 > cox3 > cox2 > atp6 > nad6 > nad3 > nad4L. The rank order of the 12 PCGs of D. mehari by length was as follows: cox1 > nad5 > nad4 > cytb > nad1 > nad2 > cox3 > cox2 > atp6 > nad6 > nad3 > nad4L. Among the 12 PCGs, the length of cox1 in D. mehari was 1,572 bp, which was the longest than those of other Paramphistomoidea trematodes, including D. japonicus (1,548 bp), D. nigromaculati (1,548 bp), F. elongatus (1,542 bp), F. cobboldi (1,542 bp), G. crumenifer (1,542 bp), C. microbothrioides (1,542 bp), O. streptocoelium (1,542 bp), P. cervi (1,542 bp), P. leydeni (1,542 bp) and E. explanatum (1,547 bp) (Table 2). Moreover, the nad4 genes of three Diplodiscus species (D. japonicus, D. mehari and D. nigromaculati) with length of 1287 bp, were 6 bp longer than those of others in the same superfamily (Yan et al., 2013; Ma et al., 2015; Yang et al., 2016; Zhao et al., 2017; Han et al., 2020).

Table 2

Comparison of 12PCGs among Diplodiscus japonicus, Diplodiscus mehari and other superfamily Paramphistomoidea trematodes.

Genes	D. j	D. m	D. n	F. e	F. c	G. c	C. m	O. s	P. c	P. l	E. e	Nucleotides similarity (%)	Amino acid similarity (%)
	No. nucleotides/No. amino acid
cox3	645/133	645/138	645/132	645/124	645/126	645/127	645/128	645/123	645/131	645/129	645/129	62.3–89.3	83.3–98.8
cytb	1113/236	1116/328	1116/236	1113/237	1113/237	1113/234	1116/239	1113/244	1113/240	1113/240	1113/241	72.0–91.3	71.5–82.8
nad4l	264/50	264/50	264/49	264/45	264/45	264/43	264/45	264/42	264/44	264/44	264/43	68.2–92.3	80.5–100.0
nad4	1287/270	1287/267	1287/267	1281/259	1281/263	1281/258	1281/262	1281/258	1281/264	1281/262	1281/260	66.7–90.4	81.5–97.6
atp6	516/91	516/94	516/94	516/97	516/95	516/94	516/96	516/98	516/94	516/98	516/96	69.6–92.2	86.0–97.7
nad2	870/180	870/177	870/175	876/172	873/171	858/169	873/174	858/167	873/174	873/183	876/169	64.8–89.7	70.1–84.9
nad1	903/190	903/189	903/193	897/178	897/178	903/184	897/180	897/181	897/184	897/178	897/182	72.6–90.5	70.0–91.3
nad3	354/72	354/75	354/74	357/71	357/72	378/80	357/69	356/71	357/72	357/70	357/70	60.8–90.7	81.5–97.6
cox1	1548/334	1572/339	1548/337	1542/330	1542/334	1542/335	1542/334	1542/331	1545/336	1545/349	1542/334	72.9–92.5	70.5–92.2
cox2	585/130	585/128	585/128	582/123	582/124	582/122	585/125	582/127	579/123	582/131	582/126	68.5–94.5	84.7–97.4
nad6	456/99	456/98	456/98	501/109	453/94	453/96	453/94	453/94	453/98	453/96	453/95	67.3–90.2	57.3–98.2
nad5	1563/325	1563/325	1563/322	1581/332	1581/334	1581/325	1581/328	1581/332	1581/331	1584/331	1581/331	67.5–89.3	49.3–67.4
Total nt/aa	10104/2110	10131/2117	10107/2105	10155/2077	10089/2073	10116/2067	10110/2074	10088/2068	10104/2091	10110/2111	10109/2076	69.6.–90.2	63.3–96.3
Total size (bp)	14179	14210	14697	14120	14256	14801	14028	13800	14023	14050	13968	62.1–89.9

Note: D. j, Diplodiscus japonicus, D. m, Diplodiscus mehari, D. n, Diplodiscus nigromaculati, F. e, Fischoederius elongatus, F. c, Fischoederius cobboldi, G. c, Gastrothylax crumenifer, C.m, Calicophoron microbothrioides, O.s, Orthocoelium streptocoelium, P. c, Paramphistomum cervi, P. l, Paramphistomum leydeni, E. e, Explanatum explanatum.

A + T contents of cox3, cytb, nad4, atp6, nad2, cox1, rrnL, and nad5 genes of D. mehari were lower than those of the other 10 species in Paramphistomoidea. However, the A + T contents of the nad3 in E. explanatum were higher than those of the other 10 species (Fig. 2A). The D. japonicus and D. mehari AT-skew values were negative, and the GC-skew values were positive. The AT-skew values of the D. japonicus mt genome ranged from −0.11 (rrnL) to −0.53 (nad6), and the GC-skew values ranged from 0.30 (rrnS) to 0.66 (nad2). The AT-skew values the D. mehari mt genome ranged from −0.14 (rrnS) to −0.57 (nad6), the GC-skew values ranged from 0.34 (rrnS) to 0.55 (nad3) (Fig. 2B). The AT-skew values of atp6 and nad6 of D. mehari were lower while its cox3 and nad4 were higher than other 10 species in Paramphistomoidea (Fig. 2B).

Among the 12 PCGs of the present D. japonicus, GTG was the most common start codon (7/12), and TAG was the predominant stop codon (7/12), which were similar with those of other Paramphistomoidea trematodes, such as G. crumenifer, P. cervi and O. streptocoelium (Table 1). Among the 12 PCGs, ATG (6/12) and GTG (6/12) were both the most common start codons, and TAG (7/12) was the predominant stop codon of D. mehari (Table 1). Disregarding the stop codons, there were 3,356 and 3,365 amino acids in 12 PCGs of D. japonicus and D. mehari, respectively. The most frequently used codons in 12 PCGs of D. japonicus were UUU/Phe, UUG/Leu2, and GUG/Val, while the least used codons were GAC/His, CUC/Leu1, and AUC/Ile. The most frequently used codons in 12 PCGs of D. mehari and were UUG/Leu2, UUU/Phe, and GUG/Val, while the least used codons were CUC/Leu, CGC/Arg, and ACC/Thr. The most frequent amino acids in the 12PCGs of D. japonicus and D. mehari were leucine (L1+L2) and serine (S1+S2) Codon usage and relative synonymous codon usage (RSCU) of both Diplodiscus mitogenomes are presented in Fig. 3. Preferable codons were commonly uncovered with important functional gene regions, as those bias codons with silent sites were found to be related to maximize the translation efficiency (Romero et al., 2000). Codons ending with A or T are used more frequently than those ending with G or C. The interesting results may reveal that the PCGs of the D. japonicus and D. mehari mtDNA genomes are biased toward utilizing T-rich amino acid codons, which suggests the nucleotide bias. However, the function of bias of codon usage for the mt system of parasite is still unclear.

Fig. 3

Relative synonymous codon usage (RSCU) of 12 protein coding genes of Diplodiscus japonicus and Diplodiscus mehari. The termination codon is not given.

Sequence identities of the 12 PCGs of the 11 species in Paramphistomoidea were 69.6.%–90.2% at the nucleotide level, and 63.3%–96.3% at the amino acid level. The complete mt genome nucleotide identities among the 11 trematodes ranged from 62.1% to 89.9%. Among the 12 PCGs, cox3 had the fewest similarities among the 11 sampled species, whereas cox1 had the highest similarity (Table 2). This result was similar to the findings in the study of Na et al. (2016), who found that cox1 has the least variation in Opisthorchiidae (Na et al., 2016).

Diversity and mutation rate among Diplodiscus spp. Mitogenomes

To determine the highly conserved and variable mitochondrial genes among 11 Paramphistomoidea trematodes, a sliding window analysis was conducted by the concatenated nucleotide sequence of 12 PCGs. The results showed obvious differences in 12 PCGs of 11 Paramphistomoidea trematodes. By computing the number of variable positions per unit length of gene, the curve indicated that cox1 was the lowest variable gene, while nad2, nad4, nad5 and cox3 showed high sequence divergence (Fig. 4). The gene cox1 was considered to be a useful barcode for metazoans, and widely employed for trematode studies (Hebert et al., 2003; Brabec et al., 2015). The mitochondrial genes cox1 has also been used to study the population genetic structure of Diplostomum spp. on a local and global scale (Brabec et al., 2015). As shown in Fig. 5, non-synonymous/synonymous (dN/dS) mutation rates among the 12 PCGs of the three Diplodiscus mitogenomes were analyzed. It was clearly showed that all PCGs were under negative (purifying) selection (dN/dS < 0.5), and cox1 had the lowest dN/dS value in 12 PCGs, which was similar with Gyrodactylus derjavinoides reported by Huyse et al. (2008). Because the stronger purification selection pressure has the smaller value of dN/dS, which further corroborates that cox1 has higher synonymous variation. Additionally, the relatively looser selection pressure (0.3 < dN/dS < 0.5), which may due to the accu-mulation of non-synonymous substitutions. These results suggest that the cox1 gene should be considered as optimal candidates for genetic marker to be used for population genetics and species identification studies in the superfamily Paramphistomoidea species.

Fig. 4

Sliding window analysis of the complete mt genome sequences of 11 Paramphistomoidea trematodes. A sliding window of 300 bp (in 10 bp overlapping steps) was used to estimate nucleotide diversity Pi (π) across the alignments. Nucleotide diversity was plotted against the mid-point positions of each window. Each gene boundary is identified.

Fig. 5

Proportions between rates of non-synonymous (dN) and synonymous (dS) nucleotide substitutions (dN/dS). Bar chart for pairwise proportions of dN/dS for each of the mitochondrial subunits of the Diplodiscus spp.

Phylogenetic analyses

The topologies of the two trees were identical based on ML and BI (Fig. 6). The phylogeny was divided into two large clades, one clade was Schistosomatidae, and the other clade contained 30 members of 13 families. Except for the paraphyletic suborder Diplostomata, other five suborders form monophyletic branch. In the Diplostomata branch, the families Clinostomidae and Dipostomidae clustered together instead of grouping with Schistosomatidae, which was consistent with previous studies (Locke et al., 2018; Li et al., 2019; Guo et al., 2022). These results indicate that Clinostomidae and Dipostomidae taxonomy require further study.

Fig. 6

Phylogenetic relationships of Diplodiscus japonicus and Diplodiscus mehari with other 30 representative Digenea trematodes based on concatenated amino acid sequences of 12 protein coding genes analyzed by maximum likelihood (ML) and Bayesian inference (BI) using Gyrodactylus salaris as the outgroup. Statistical support values (Bootstrap/posterior probability) of ML/BI analysis are shown above the nodes. Circles indicate ML/BI = 100/1.0, other values are given above the nodes. Suborders and families are highlighted by individual colors. Accession numbers are given for each species at the end of each sequence. The scale bar corresponds to the estimated number of substitutions per site. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

In the suborder Pronocephalata, Diplodiscidae formed a group, while Paramphistomidae and Gastrothylacida formed the other group. It indicated that Paramphistomidae and Gastrothylacidae were more closely related than Diplodiscidae. In the group of Paramphistomidae and Gastrothylacidae, the five species of Paramphistomidae did not cluster together, O. streptocoelium formed a separate branch, which is consistent with previous study (Zhao et al., 2017). Interestingly, the further analysis found that the flukes parasitic in amphibians (frogs) formed one group, and the flukes from ruminants (cattle, sheep, ect) clustered together and formed another group. D. japonicus and D. mehari were more closely related than the D. nigromaculati. However, only three Diplodiscus species complete mt genomes have been published in GenBank. Therefore, further studies are required to research the phylogenetic relationships of Paramphistomoidea by sequencing the mt genomes of other Diplodiscidae flukes from different regions and hosts.

In conclusion, the complete mt genomes of D. japonicus and D. mehari were determined for the first time in our study, and both comparative complete mt genome sequences analyses and phylogenetic analyses suggested that Paramphistomidae and Gastrothylacidae were more closely related than Diplodiscidae. And the flukes parasitic in amphibians formed one group, and the flukes from ruminants gathered together form another group. The mt genome data of D. japonicus and D. mehari provide new and useful genetic markers for further studies of Paramphistomoidea trematodes.

Declaration of competing interest

The authors report no conflicts of interests.