Characterization of the complete mitochondrial genome sequence of Tracheophilus cymbius (Digenea), the first representative from the family Cyclocoelidae

Abstract

Tracheophilus cymbius (Trematoda: Cyclocoelidae) is a common tracheal fluke of waterfowl, causing serious loss in the poultry industry. However, taxonomic identification of T. cymbius remains controversial and confused. Mitochondrial (mt) genomes can provide genetic markers for the identification of closely related species. We determined the mt genome of T. cymbius and reconstructed phylogenies with other trematodes. The T. cymbius mt genome is 13,760 bp in size, and contains 12 protein-coding genes (cox 1–3, nad 1–6, nad 4L, cyt b and atp 6), 22 transfer RNA (tRNA) genes, two ribosomal RNA genes and one non-coding region. All are transcribed in the same direction. The A + T content is 62.82%. ATG and TAG are the most common initiation and termination codons, respectively. Phylogenetic analyses of concatenated nucleotide sequences show T. cymbius grouping in suborder Echinostomata, and clustering together, with high statistical support, as a sister taxon with Echinochasmus japonicus (Echinochasmidae), the two forming a distinct branch rooted to the ancestor of all Echinostomatidae and Fasciolidae species. This is the first report of the T. cymbius mt genome, and the first reported mt genome within the family Cyclocoelidae. These data will provide a significant resource of molecular markers for studying the taxonomy, population genetics and systematics of trematodes.

Introduction

Tracheophilus cymbius (Diesing, 1850) Skrjabin, 1913, belongs to the family Cyclocoelidae (Stossich, 1902), and is parasitic in the trachea of ducks and other waterfowl, often causing mild clinical symptoms, including cough, malnutrition, weight loss and dyspnoea, but sometimes even causing death (Bisseru, 1957; Tang & Tang, 1978; Scott et al., 1982). Tracheophilus cymbius have been reported in Asia, Europe, America and Africa (Bisseru, 1957; Gu et al., 1973; Tang & Tang, 1978; Scott et al., 1982; Wang et al., 2004). Tracheophilus cymbius-infected domestic ducks and wild aquatic birds have also been widely recorded in China, including Fujian, Hebei, Ningxia, Jiangxi and Heilongjiang Provinces (Gu et al., 1973; Tang & Tang, 1978; Wang et al., 2004). The prevalence of T. cymbius in domestic ducks ranges from 5.8% to 30% in Fujian Province, and from 2.6% to 37.8% in wild aquatic birds in Bai-Yang-Dian Lake, Hebei Province (Gu et al., 1973; Tang & Tang, 1978). Tracheophilus cymbius is particularly common and prevalent in domestic ducks in Heilongjiang Province, causing severe clinical symptomology. Our research shows 5.4% (12/222) T. cymbius prevalence in Heilongjiang Province in 2003, and 22.3% (21/94) in 2017, with an infection intensity of 1–8 (Wang et al., unpublished data), causing significant economic loss to the duck industry. Therefore, effective control and prevention of this trematode is important to the poultry industry.

Traditional approaches for the identification and differentiation of parasites based on morphological features have long been used worldwide. However, morphological approaches have limitations in identifying and distinguishing closely related species (McManus & Bowles, 1996). Although the family Cyclocoelidae was first named over a century ago, the status of the family and several cyclocoelid species remains controversial and uncertain. This is particularly the case with the important tracheal fluke T. cymbius Skrjabin (1913), considering that Tracheophilus Skrjabin and Typhlocoelum Stossich were classified into two different genera based on testes shape (Skrjabin, 1913; Bisseru, 1957). However, Joyeux & Baer (1927) recognized Tracheophilus Skrjabin and Typhlutimum Witenberg as synonyms of Typhlocoelum Stossich, 1902, and the two species were combined into the same genus, namely Typhlocoelum cymbius (Diesing, 1850) (syn. Tracheophilus sisowi Skrjabin, 1913, Typhlocoelum sisowi Skrjabin, 1913) and Typhlocoelum cucumerinum (Rudolphi, 1809) (Joyeux & Baer, 1927; Bisseru, 1957). Tang & Tang (1978) reported that T. cymbius and T. sisowi should be considered the same species based on life cycle, in spite of the two trematodes having one key difference in morphological features (T. cymbius with the posterior vitelline united, versus T. sisowi with a separated posterior vitelline) (Tang & Tang, 1978). Scott et al. (1982) used Typhlocoelum cucumerinum sisowi (Skrjabin, 1913) and Typhlocoelum cucumerinum (Rudolphi, 1809) metacercariae, identified by morphological observation, to infect wild waterfowl, and the results suggested that the two parasites are not separate species, and rather should be considered two separate sub-species (Scott et al., 1982).

Mitochondrial (mt) genome and nuclear ribosomal DNA (rDNA) sequences can effectively identify species of the parasite (Wang et al., 2011; Gao et al., 2017). However, only a partial Typhlocoelum sp. 28S ribosomal RNA (rRNA) sequence (KT956960) is currently available in GenBank, and no complete mt genome data for any parasite belonging to the family Cyclocoelidae is available. This is a crucial limitation into investigations of the systematics and phylogeny of this family.

Thus, the purpose of our study was to determine the complete T. cymbius mt genome, and to analyse the phylogenetic relationships of T. cymbius with other trematodes based on a dataset consisting of the concatenation of 12 protein-coding gene nucleotide sequences. The results will help resolve issues of taxonomy, population genetics and systematics within and beyond the fluke family Cyclocoelidae.

Materials and methods

Parasite and DNA extraction

Adult T. cymbius flukes were collected from the trachea of naturally infected ducks in Daqing, Heilongjiang Province, China. Specimens were washed in physiological saline, and then morphologically identified to the species level (Tang & Tang, 1978). Next, the flukes were fixed in 70% ethanol and stored at −20°C until further use. Total genomic DNA was extracted from individual worms with a TIANamp Genomic DNA Kit (Tiangen, Beijing, China), according to the manufacturer's protocol.

Amplification and sequence analysis

The entire mt genome of a single T. cymbius specimen was amplified in seven overlapping fragments using primers (supplementary table S1) designed from relatively conserved regions of the mt genome nucleotide sequences of closely related species. Polymerase chain reaction (PCR) cycling conditions used to amplify the T. cymbius mtDNA genome were based on a previous report (Li et al., 2019). PCR products were sent to Sangon Biotech Company (Shanghai, China) for sequencing in both directions using the same primers.

Sequences were assembled and aligned manually against the complete Fasciola gigantica (NC024025) and Fasciola hepatica (NC002546) mt genome sequences to identify gene boundaries using the program DNAStar v. 5.0 (Burland, 2000). Each protein-coding gene was translated into amino acid sequences using the trematode mt genetic code in MEGA X (Kumar et al., 2018). The secondary structures of the 22 predicted tRNA genes were estimated using the online program ARWEN (http://130.235.46.10/ARWEN/), and/or by visual identification combined with manual proofreading. The rRNA genes were identified by comparison with the mt genomes of closely related species.

Using the individual protein-coding genes of T. cymbius, Echinochasmus japonicus, F. hepatica and Echinostoma miyagawai, all members of Echinostomata, comparisons were made based on complete mt genome size, gene arrangement, A + T content, AT/GC-skew and nucleotide and amino acid sequence differences. Nucleotide and amino acid sequences differences were calculated using DNAStar v. 5.0 (Burland, 2000). The AT-skew and GC-skew values in both coding and non-coding regions (NCRs) were calculated using the following equations: AT-skew = (A − T)/(A + T), and GC-skew = (G − C)/(G + C).

Phylogenetic analyses

Nucleotide sequences from the 12 protein-coding genes in the T. cymbius mt genome were concatenated and aligned with those from published trematodes mt genomes using MEGA X (Kumar et al., 2018). We used the Gblocks online server (http://molevol.cmima.csic.es/castresana/Gblocks_server.html) to exclude ambiguously aligned regions from the multiple nucleotide sequence alignment, specifying the ‘less stringent’ selection option (Castresana, 2000). We used 38 trematode mt genomes available in GenBank, in addition to our T. cymbius data, to create our phylogenetic datasets (see supplementary table S2). Gyrodactylus salaris (NC008815) was included as an outgroup. Phylogenetic trees were reconstructed using the Bayesian inference (BI) and maximum likelihood (ML) methods on the concatenated nucleotides sequences of the 12 protein-coding genes.

BI was performed using the mixed model in MrBayes v. 3.1.1 (Ronquist & Huelsenbeck, 2003) with 1,000,000 metropolis-coupled Markov chain Monte Carlo generations (Ronquist & Huelsenbeck, 2003). The first 250 trees were omitted as burn-in, and the remaining trees were used to calculate Bayesian posterior probabilities. We used PAUP v. 4.0 Beta 10 (Swofford, 2002), with 100 random addition searches and tree-bisection–reconnection branch swapping, for the nucleotide ML analysis. ML bootstrap reliability estimates were calculated from 100 bootstrap replicates with ten random additions per replicate, also using PAUP (Swofford, 2002). Phylograms were drawn using Tree View v. 1.65 (Page, 1996) for all inferences.

Results and discussion

Tracheophilus cymbius mtDNA features

The whole T. cymbius mt genome (GenBank accession code MK355447) is a typical circular mtDNA molecule 13,760 bp in size (supplementary fig. S1). The mt genome contains 36 genes: 12 protein-coding genes (cox 1–3, nad 1–6, nad 4L, atp 6 and cyt b), 22 tRNA genes, two rRNA genes (rrn L and rrn S) plus one NCR (table 1 and supplementary fig. S1). All genes are transcribed in the anticlockwise direction. The T. cymbius mt genome gene arrangement is identical to those of F. hepatica, Clinostomum complanatum and some Opisthorchiidae species (Opisthorchis felineus, Metorchis orientalis and Clonorchis sinensis) (Le et al., 2001; Cai et al., 2012; Chen et al., 2016; Na et al., 2016). The entire T. cymbius mt genome nucleotide composition is biased toward A and T, with an overall A + T content of 62.82% (table 2). There is a very low base C content (9.75%) in the T. cymbius mt genome. Furthermore, a 40 bp overlap between the nad 4L and nad 4 genes exists in T. cymbius mtDNA, which is consistent with most flukes, but longer than those of Schistosoma haematobium (28 bp) and Schistosoma japonicm (37 bp), and shorter than that of Schistosoma mekongi (64 bp) (Le et al., 2001, 2016; Littlewood et al., 2006; Cai et al., 2012; Liu et al., 2016).

Table 1. Mitochondrial genome organization of Tracheophilus cymbius.

aa, amino acid; nt, nucleotide.

Table 2. Comparison of complete mtDNA among Tracheophilus cymbius and other Echinostomata species.

T. c, Tracheophilus cymbius; Es. m, Echinostoma miyagawai; Ec. j, Echinochasmus japonicus; F. h, Fasciola hepatica; nt s, nucleotides similarity; aa s, amino acid similarity.

The lengths of the 12 protein-coding T. cymbius mtDNA genes arrange in the order nad 5 > cox 1 > nad 4 > cyt b > nad 1 > nad 2 > cox 3 > cox 2 > atp 6 > nad 6 > nad 3 > nad 4L. A total of 3372 amino acids are encoded in the T. cymbius mt genome (table 2). Four genes (nad 1, nad 3, nad 4 and nad 5) use GTG as a start codon, eight genes (atp 6, cyt b, cox 1, cox 2, cox 3, nad 2, nad 4L and nad 6) use ATG as a start codon. Additionally, all genes have complete TAG termination codons (table 1). The most frequent codon is TTT (Phe) in the 12 protein coding genes, with a frequency of 12.5%. The least frequent codons are ACC (Thr), CGA (Arg) and AGC (Ser), all with frequencies of 0.12% (supplementary table S3).

The T. cymbius rrn L gene is located between trn T and trn C, and rrn S is located between trn C and cox 2. The lengths of the rrn S and rrn L genes are 760 bp and 985 bp, respectively (table 1). The A + T contents of the rrn L and rrn S are 61.93% and 59.61%, respectively. A total of 22 tRNA genes were identified in the T. cymbius mt genome, ranging from 63 bp to 76 bp in length (table 1). Twenty-one of the tRNA gene products fold into a predicted secondary structure with the conventional cloverleaf shape. The remaining trn S1 transcript has unpaired D-arms replaced by 7-bp loops (supplementary fig. S2). This is different than Echinostoma hortense, Es. miyagawai and F. hepatica, in which the two trn S transcripts both have unpaired D-arms (Le et al., 2001; Liu et al., 2016; Li et al., 2019). The T. cymbius mt genome NCR is located between trn G and cox 3, with a length of 142 bp (table 1).

Comparative mt genome analysis: T. cymbius and other closely related species

The complete T. cymbius mt genome is shorter than that in Ec. japonicus (15,865 bp), F. hepatica (14,462 bp) and Es. miyagawai (14,413 bp) (table 2). This length discrepancy is mainly due to the NCR lengths: 152 bp in T. cymbius, 931 bp in Es. miyagawai, 2342 bp in Ec. Japonicus and 806 bp in F. hepatica. Full-length mt genome nucleotide sequence identities range from 67.6% to 72.0% among the four Echinostomata trematodes (T. cymbius, Ec. japonicus, F. hepatica and Es. miyagawai) (table 2). Among them, the highest full-length nucleotide sequence identity (69.1%) is between T. cymbius and F. hepatica (70.7% nucleotide identity and 64.8% amino acid identity in the 12 protein-coding gene sequences) (supplementary table S4). The total A + T content (62.82%) is slightly higher than those of Ec. japonicus (61.48%) and F. hepatica (62.18%), but lower than that of Es. miyagawai (65.3%). ATG and TAG are the most frequently used initiation and termination codons, respectively, which is the case in many other trematodes, including F. hepatica, Es. miyagawai and Ec. japonicus (Le et al., 2001, 2016).

The AT/GC-skews in each mt genome gene or region of the four Echinostomata trematodes are listed in supplementary table S5. AT-skew values of the four trematodes are generally negative, and the majority of GC-skew values are positive. However, the AT-skew value for the T. cymbius NCR (−0.500) is lower than the other three (Ec. japonicus, F. hepatica and Es. miyagawai), which range from −0.386 to −0.214. GC-skew values in the four species have a large range of variation, with the lowest (0.238) in Ec. japonicus and the highest (0.676) in T. cymbius. Although some genes in T. cymbius and Ec. japonicus possess different skew values, the overall pattern of nucleotide skew in these two species is quite close.

Phylogenetic analyses

Phylogenetic analyses of the nucleotide sequences of the 12 concatenated coding regions from 39 trematode mt genomes were performed using two approaches (BI and ML).

The phylogenetic results of the two methods are somewhat different (See figs 1 and 2). The BI phylogenetic tree splits into two large clades. One clade contains seven members of Schistosomatidae, and the other clade contains 32 members of another 15 families (fig. 1). Interestingly, although the families Dipostomidae and Clinotomidae are both traditionally placed in the suborder Diplostomata, these species did not group with Schistosomatidae, which is also traditionally placed in the order Diplostomida. Rather, these species cluster basal within the order Plagiorchiida, paraphyletic to Diplostomida. This result is consistent with a previous study using mt sequences and ultra-conserved genomic elements to analyse the validity of the Diplostomoidea and Diplostomida (Locke et al., 2018). However, our result is inconsistent with another earlier study using nuclear small-subunit ribosomal and large-subunit ribosomal DNA sequences to analyse the phylogeny of Digeneans (Olson et al., 2003). In the non-Diplostomata clade of our BI tree, trematodes of the suborder Pronocephalata form a clade basal to all the other suborders (just more recent than the aberrant Dipostomidae and Clinotomidae clade). In the Echinostomata branch, all of the Echinostomatidae species cluster together within the Echinostomatidae clade, except Es. hortense. Echinostoma hortense clusters with the Fasciola, Fascioloides and Fasciolopsis species used in our study. The result is consistent with a recent study by Li et al. (2019), who used mtDNA sequences to analyse the phylogenetic relationships of Es. miyagawai with other species of Echinostomata (Li et al., 2019).

Fig. 1. Genetic relationships of Tracheophilus cymbius with other representative trematodes based on mitochondrial nucleotide data. Phylogenetic analysis based on the concatenated mitochondrial sequence data representing 12 protein-coding genes was conducted using Bayesian inference (BI), with Gyrodactylus salaris as the outgroup.

Fig. 2. Genetic relationships of Tracheophilus cymbius with other representative trematodes based on mitochondrial nucleotide data. Phylogenetic analysis based on the concatenated mitochondrial sequence data representing 12 protein-coding genes was conducted using maximum likelihood (ML), with Gyrodactylus salaris as the outgroup.

The ML phylogenetic (fig. 2) tree differs from the BI phylogenetic tree in a couple of key respects. Suborders Echinostomata, Opisthorchiata, Troglotremata, Xiphidiata and Pronocephalata form independent, monophyletic groups in both phylogenetic trees. Furthermore, the family Schistosomatidae forms a monophyletic clade in both phylogenetic trees, yet the suborder Diplostomata (Dipostomidae + Clinostomidae + Schistosomatidae) is an independent monophyletic clade in the ML tree and paraphyletic in the BI tree. Both methods' trees confidentially group T. cymbius and Ec. japonicus as sister taxa with high support values (BI 1.0, ML 100%), basal to the other Echinosomata. Both phylogenies show T. cymbius to be most closely related to Echinochasmidae, with a common ancestor basal to a clade containing Echinostomatidae and Fasciolidae. This is inconsistent with a previous study using the nuclear 28S rRNA gene (Tkach et al., 2016) in which the Cyclocoelidae are more closely related to Echinostomatidae and Fasciolidae than to Echinochasmidae. Inconsistency between nuclear and mt phylogenies is not uncommon though, perhaps here due to a lower phylogenetic signal/noise ratio in the mt genomes of digeneans than in nuclear genomes, aggravated by the effect of incomplete taxon sampling (Philippe et al., 2011; Locke et al., 2018).

In conclusion, the T. cymbius complete mtDNA genome sequence has been determined and reported for the first time. This represents the first mt genome available from any member of the family Cyclocoelidae. These data will provide novel mtDNA markers for studying the molecular epidemiology and population genetics of Cyclocoelidae and other trematodes.