Introduction
Toxocaridae contains only two genera, Porrocaecum and Toxocara (Gu et al. Reference Gu, Guo, Chen, Sitko, Li, Guo and Li2023). Toxocara species can cause toxocariasis and commonly occur in wildlife and domestic animals. The larvae of some Toxocara species can also accidentally infect humans; they are therefore of veterinary, medical, and economic significance. Larvae also can migrate into host tissue, leading to tissue damage, inflammatory reactions, visual impairment, and blindness (Zhou et al. Reference Zhou, Guo, Deng, He, Ouyang and Wu2020). T. canis can cause detrimental damage to the brain of intermediate or paratenic hosts (Chen et al. Reference Chen, Qiu and Mo2022). Therefore, identifying different Toxocara is conducive to preventing and treating nematode infection in humans and animals.
In Toxocara, T. apodemi and T. mackerrasae are non-zoonotic species that are host-specific parasites of Muridae (Olsen Reference Olsen1957; Warren Reference Warren1972; Asakawa et al. Reference Asakawa, Li, Guo, Yang, Huhebateer Liu, Liu, Cao and Chen1994). T. apodemi mainly infects (Apodemus. peninsulae) A. peninsulae and A. agrarius and has only been reported in China and Korea (Ziegler and Macpherson Reference Ziegler and Macpherson2019; Kim et al. Reference Kim, Hong, Ryu, Park, Cho, Yu, Chae, Choi and Park2020). The nematode was first reported in A. peninsulae from Korea as a new species of the genus Neoascaris and was named Neoascaris apodemi by Olsen in 1957 (Olsen Reference Olsen1957). The parasite was then revised as T. apodemi after it was reported in striped field mice (A. agrarius) in China. Their morphologic characteristics are that external prolongations of the labial pulp are asymmetrical and rounded at the anterior border (Ziegler and Macpherson Reference Ziegler and Macpherson2019).
Despite these nematodes’ ubiquity and essential roles in diverse ecological systems, the origin and early evolutionary history of these species have long been matters of debate. In recent years, mitochondrial (mt) genomes, as genetic markers, have been widely used to analyse the taxonomy and diversity of certain taxa or specific groups because of their strict maternal inheritance, apparent lack of recombination, rapid evolutionary rate, and comparatively conserved genomic structure (Zhou et al. Reference Zhou, Guo, Deng, He, Ouyang and Wu2020; Gao et al. 2020). The study of the complete mitochondrial gene in Toxocaridae will help us to understand its evolutionary relationships. In 1992, the first paper on the complete mt genome of Ascaris suum was published, and since then, more and more mt genomes of other Ascaridoidea species have been sequenced and annotated (Okimoto et al. Reference Okimoto, Macfarlane, Clary and Wolstenholme1992; Liu et al. Reference Liu, Wu, Song, Wei, Xu, Lin, Zhao, Huang and Zhu2012). The phylogenetic relationships within the superfamily were gradually refined. For example, based on mt genome sequences, Baylisascaris species (B. schroederi, B. ailuri, and B. transfuga) are more closely related to A. suum than to the 3 Toxoara species (T. canis, T. cati, and T. malaysiensis) and A. simplex (Xie et al. Reference Xie, Zhang, Wang, Lan, Li, Chen, Fu, Nie, Yan, Gu, Wang, Peng and Yang2011).
Currently, the available molecular information for T. apodemi is limited, with only the partial sequence of a small subunit of a ribosomal RNA gene published in GenBank (Kim et al. Reference Kim, Hong, Ryu, Park, Cho, Yu, Chae, Choi and Park2020). Thus, the objectives of this study were to describe and determine its complete mt genome sequences, compare its mt genome with those of other Toxocara species, and reconstruct its phylogenetic relationships to assess its systematic and phylogenetic position.
Materials and methods
Parasites and species identification
Adult nematodes of T. apodemi were obtained from the intestine of a naturally infected wild mouse in Taizhou, Zhejiang Province, China. The parasite was identified at the species level according to previously described morphological features (Asakawa et al. Reference Asakawa, Li, Guo, Yang, Huhebateer Liu, Liu, Cao and Chen1994; Kim et al. Reference Kim, Hong, Ryu, Park, Cho, Yu, Chae, Choi and Park2020; Olsen Reference Olsen1957). The molecular characteristics of 18S and internal transcribed spacers (ITS) were amplified using universal primers NC5 (5′-GTAGGTGAACCTGCGGAAGGATCATT-3′) and NC2 (5′-TTAGTTTCTTTTCCTCCGCT-3′) for ITS sequence (Gasser et al. Reference Gasser, Bott, Chilton, Hunt and Beveridge2008); the genus Toxocara’s universal primers 18S-F (5´-GCTAATACATGCACCAAAGC-3´) and 18S-R (5´-GATCACGGAGGATTTTCAAC-3´) were reported previously for 18S rDNA sequence (Kim et al. Reference Kim, Hong, Ryu, Park, Cho, Yu, Chae, Choi and Park2020).
Polymerase chain reaction (PCR) amplification of mitochondrial genome and sequencing
In total, 14 specific primers were used to amplify the complete mt genome of T. apodemi and were designed based on those of T. canis (Accession: NC_010690.1) and T. cati (NC_010773.1) published in GenBank (Table S1). PCR reactions were carried out under the following conditions: 94°C for 2 min, then 94°C for 1 min, 40–56°C for 30 s, and 72°C (~1 kb region) for 30 s for 35 cycles, with a final extension at 72°C for 7 min. The positive PCR products were cloned to PMD-18T vectors and then transferred into DH5α cells for positive plasmid sequencing at Sangon Biotech Company (Shanghai, China).
Genome sequence assembly and gene annotation
Sequences were assembled manually and aligned against the mt genome sequences from Toxocara species using the program Clustal X (v. 1.83) and the MegAlign procedure within the DNAStar (Burland Reference Burland2000; Larkin et al. Reference Larkin, Blackshields, Brown, Chenna, McGettigan, McWilliam, Valentin, Wallace, Wilm, Lopez, Thompson, Gibson and Higgins2007). The boundaries of 12 protein-coding gene (PCG) sequences, the 16S ribosomal RNA (rrnL) gene, and the 12S ribosomal RNA (rrnS) gene were determined using the mt genome sequences of other Toxocara nematodes available in GenBank using the software MEGA X (Li et al. Reference Li, Lin, Song, Wu and Zhu2008; Kumar et al. Reference Kumar, Stecher, Li, Knyaz and Tamura2018; Xie et al. Reference Xie, Wang, Chen, Wang, Zhu, Hu, Han, Wang, Zhou and Zuo2022). A total of 22 transfer RNA (tRNA) genes were identified using the online tool tRNA scan-SE (http://lowelab.ucsc.edu/tRNAscan-SE) or via visual inspection. Finally, the circular genomic maps of T. apodemi were generated using the CGView online server V1.0 (http://stothard.afns.ualberta.ca/cgview_server/).
Comparative mt genome sequences analysis
Comparisons were made among the complete mt genomes available in GenBank for five Toxocara nematodes, including gene lengths, A+T contents, nucleotides, and amino acid sequence identities. The A+T content of complete mt genomes, 12 PCGs, and the 1st, 2nd, and 3rd coding positions were computed using DNAStar (v. 12.1) (Burland Reference Burland2000). The relative synonymous codon usage (RSCU) values of the 12 PCGs of five Toxocara species were calculated with MEGA X (Kumar et al. Reference Kumar, Stecher, Li, Knyaz and Tamura2018). The p-distance model of MEGA X was used for the genetic distance analysis of 12 PCGs among five Toxocara species. The rate of non-synonymous substitutions (Ka) and the rate of synonymous substitutions (Ks) were used to predict evolutionary processes; Ka/Ks ratios were calculated for the nucleotide sequences of all 12 mt PCGs of T. canis, T. cati, T. malaysiensis, T. vitulorum, and T. apodemi using DnaSP v5 (Librado and Rozas Reference Librado and Rozas2009). A sliding window of 300 bp (in 10 bp overlapping steps) was used to estimate nucleotide diversity Pi (π) across the alignment. Nucleotide diversity was plotted against the mid-point positions of each window.
Phylogenetic analysis
Phylogenetic relationships were reconstructed with Bayesian inference (BI) and maximum likelihood (ML) using the concatenated amino acid sequences of 12 PCGs from mt genomes available in GenBank for 35 Ascaridomorpha species, using Caenorhabditis elegans as the outgroup (Table S2). The amino acid sequences of 12 PCGs of 36 nematodes were aligned using MAFFT 7.471 and then concatenated into a single alignment (Katoh and Standley Reference Katoh and Standley2013). Sites of ambiguous alignment were eliminated using the Gblocks online server (http://www.phylogeny.fr/one_task.cgi?task_type=gblocks). MrBayes 3.1 was used to reconstruct the BI tree and four independent Markov chain runs were performed for 1,000,000 metropolis-coupled MCMC generations, sampling a tree every 100 generations. The mixed model was selected as the best model and was performed using the BI method. The first 25% (2500) of trees were omitted as burn-in, and the remaining trees were used to calculate Bayesian posterior probabilities (Ronquist and Huelsenbeck Reference Ronquist and Huelsenbeck2003). The phylograms were drawn using FigTree v1.4.4 (http://tree.bio.ed.ac.uk/software/figtree/). An ML tree was inferred by using the JTT matrix-based model and performed using MEGA X (Jones et al. Reference Jones, Taylor and Thornton1992; Kumar et al. Reference Kumar, Stecher, Li, Knyaz and Tamura2018). The tree with the highest log likelihood (‒66471.74) is shown. The percentage of trees in which the associated taxa clustered together is shown next to the branches. Initial tree(s) for the heuristic search were obtained automatically by applying the Neighbor-Join and BioNJ algorithms to a matrix of pairwise distances estimated using a JTT model, and then the topology with the superior log likelihood value was selected.
Results and Discussion
Molecular and morphological identification of T. apodemi
Adult females measured 68–75 mm in length and 1.8–2.2 mm in width. Eggs were oval and measured 74–78 μm. The tail terminates in a conical, retractible spine-like structure. The surface of the egg has rough dents (Figure S1).
The 18S rRNA sequences of T. apodemi were 1352 bp in size with 99.7% identity to the corresponding available sequences in GenBank (Kim et al. Reference Kim, Hong, Ryu, Park, Cho, Yu, Chae, Choi and Park2020). The ITS sequences of T. apodemi (Accession no. OR231233) were first determined in the present study; it was 957 bp in length and had 80.94% identity to that of T. canis in GenBank (Chen et al. Reference Chen, Qiu and Mo2022).
General features of T. apodemi mitogenome
The complete mt genome sequences of T. apodemi were the first to be sequenced, with a length of 14303 bp (Accession no. OR241493), which was similar to those of T. canis (14322 bp) and T. malaysiensis (14266 bp), shorter than those of T. vitulorum (15045 bp), and longer than those of T. cati (14029 bp) (Li et al. Reference Li, Lin, Song, Wu and Zhu2008; Xie et al. Reference Xie, Wang, Chen, Wang, Zhu, Hu, Han, Wang, Zhou and Zuo2022). The complete mt genome of T. apodemi contained 12 PCGs (cox1-cox3, cytb, atp6, nad1-nad6, and nad4L), 22 tRNAs, 2 ribosomal RNAs (rRNAs), and 2 non-coding regions (NCRs) (Table 1, Figure 1). All mt genes of T. apodemi were transcribed in the same direction, which was the same as for the other species in Toxocara available in GenBank, but different from those of Ascaridia columbae and Ascaridia galli in superfamily Heterakoidea (Liu et al. Reference Liu, Shao, Li, Zhou, Li and Zhu2013a; Han et al. Reference Han, Yang, Li, Zhou, Zhou, Liu, Lu, Wang, Yang, Shi, Li, Du, Guan, Zhang, Guo, Wang, Chai, Lan, Liu, Liu, Sun and Hou2022). A total of 3,423 amino acids were encoded by 12 PCGs in the complete T. apodemi mt genome. TTG, ATT, and GTT were used as the start codon of 12 PCGs; TAA and TAG were used as the stop codon of 9 PCGs. “T” as the termination codon also appeared in atp6, nad5, and nad4L, which was the same as congeneric species. However, in the nad2 gene, TAG was the termination codon in T. apodemi, and “TA” or “T” were the termination codons in other species in Toxocara (Li et al. Reference Li, Lin, Song, Wu and Zhu2008; Xie et al. Reference Xie, Wang, Chen, Wang, Zhu, Hu, Han, Wang, Zhou and Zuo2022). The mt genomes of T. apodemi had relatively compact structures, with fewer spacer regions and short overlaps between some adjacent genes, in which 11 intergenic spacers fluctuated with 1 to 9 bp, and 5 overlaps were 1–2 bp in length (Table 1). The A+T content of the complete mt genome and the 12 PCGs of T. apodemi were 68.38% and 66.54%, respectively. This is in accordance with those of other Ascaris species, such as Toxocara species, human-type Ascaris, pig-type Ascaris, and hybrid Ascaris (Li et al. Reference Li, Lin, Song, Wu and Zhu2008; Zhou et al. Reference Zhou, Guo, Deng, He, Ouyang and Wu2020; Xie et al. Reference Xie, Wang, Chen, Wang, Zhu, Hu, Han, Wang, Zhou and Zuo2022), and slightly lower than those of Parascaris equorum in Ascaridoidea (70.25%) (Gao et al. Reference Gao, Zhang, Wang, Li, Li, Xu, Gao and Wang2019). The comparison analysis showed that the A+T content of the second coding position of the 12 PCGs was more similar than those of the first and the third coding positions among five Toxocara species (Figure 1).
Table 1. Mitochondrial genome characteristic of Toxocara apodemi
Figure 1. Gene map of T. apodemi complete mt genome.
There were 22 tRNAs in the mt genome of T. apodemi, ranging from 54 bp (trnP) to 63 bp (trnK) in size (Table 1). The estimated secondary structures of 22 tRNAs were identical to those of all other chromadorean nematodes investigated so far, except for T. spiralis (Gao et al. Reference Gao, Zhang, Wang and Zhao2022; Xie et al. Reference Xie, Wang, Chen, Wang, Zhu, Hu, Han, Wang, Zhou and Zuo2022). The lengths of the rrnL and rrnS of the T. apodemi mt genome were 955 bp and 697 bp, respectively, and they were located between trnH and nad3, and between trnE and trnS2, respectively. Moreover, there was also a NCR with 117 bp placed in nad4 and cox1, and an A+T rich with 966 bp placed in trnS2 and trnN.
Comparative mitogenomics
A comparative analysis of the nucleotide sequences of the complete mt genomes of T. apodemi with T. canis, T. cati, T. malaysiensis, and T. vitulorum showed that their identities were 84.0%, 82.7%, 83.2%, and 83.2%, respectively (Table 2). T. apodemi had higher identities with that of T. canis than those of other congeners. This level of mt genome divergence is lower than that between Ophidascaris wangi and Toxascaris leonina (19.77%), and higher than that of T. leonina from cheetah and dog (7.2%) and the Pseudoterranova decipiens species complex (3.8–9.4%) (Liu et al. Reference Liu, Nadler, Liu, Podolska, D’Amelio, Shao, Gasser and Zhu2016; Jin et al. Reference Jin, Li, Liu, Zhu and Liu2019; Zhou et al. Reference Zhou, Ma, Tang, Zhu and Xu2021). Nevertheless, the identities of the amino acid sequence of five species in Toxocara were similar at 89.8–90.8%. The discrepancy in amino acid sequences was similar to that found between O. wangi and T. leonina (12.17%) (Zhou et al. Reference Zhou, Ma, Tang, Zhu and Xu2021).
Table 2. Comparative analysis of mtDNA sequences in genus Toxocara
Note: Identity: T. apodemi vs. T. canis, T. apodemi vs. T. cati, T. apodemi vs. T. malaysiensis, T. apodemi vs. T. vitulorum.
The RSCU reflected the genetic codon usage bias to reveal the relative frequency of synonymous codons (Yang et al. Reference Yang, Huang, Wang, Yang, Zhang and Zheng2023). It is a metric commonly used to measure codon usage bias. In Ascaridomorpha, the most commonly used codons are UUG, ACU, CCU, GCU, AGA, and AGU (Han et al. Reference Han, Yang, Li, Zhou, Zhou, Liu, Lu, Wang, Yang, Shi, Li, Du, Guan, Zhang, Guo, Wang, Chai, Lan, Liu, Liu, Sun and Hou2022). In the complete mt genome of T. apodemi, UUG (RSCU = 3.82), ACU (RSCU = 3.18), CCU (RSCU = 3.18), CGU (RSCU = 3.15), and UCU (RSCU = 3.07) were the most frequently used, and CUC (RSCU = 0), CGA (RSCU = 0), CCA (RSCU = 0.05), and UUC (RSCU = 0.07) were the least used (Figure 2). This indicated that T. apodemi had similar codon bias in Ascaridomorpha. The RSCU of 5 Toxocara species were identical to the most frequently used codon but had a slight difference with the least used codon. For example, the RSCU of AGC was 0.27, 0.29, 0.16, and 0.18 in T. cati, T. malaysiensis, T. vitulorum, and T. apodemi, respectively. However, it was 0.06 in T. canis (Li et al. Reference Li, Lin, Song, Wu and Zhu2008; Meng et al. Reference Meng, Xie, Gu, Zheng, Liu, Li, Wang, Zhou, Zuo and Yang2019; Xie et al. Reference Xie, Wang, Chen, Wang, Zhu, Hu, Han, Wang, Zhou and Zuo2022).
Figure 2. Relative synonymous codon usage (RSCU) of the mitochondrial genomes and A+T contents of complete genomes, 12PCGs, and 1st, 2nd, and 3rd coding position of five Toxocara species.
Ta: T. apodemi, Tc: T. canis, Tca: T. cati, Tm: T. malaysiensis, Tv: T. vitulorum
The sliding window analysis of the nucleotide diversity (Pi values) of the 12 aligned PCGs among 5 Toxocara species revealed a high degree of nucleotide variation within different genes (Figure 3A). Nucleotide diversity values ranged from 0.09975 (cox1) to 0.16268 (nad4). Cox1 and cox2 had relatively low nucleotide diversity values, indicating that they were relatively conserved genes in the 12 PCGs of Toxocara species. Cox1, as the least variable and most slowly evolving mitogenome gene, was identical to that of Ascaris species and Cyathostominae nematodes (Zhou et al. Reference Zhou, Guo, Deng, He, Ouyang and Wu2020; Gao et al. Reference Gao, Zhang, Wang and Zhao2022). Nad4, cytb, nad2, and nad6 presented higher variability in the five Toxocara species. In a previous study, nad4 was the least conserved gene between Parascaris equorum and Parascaris univalens (Gao et al. Reference Gao, Zhang, Wang, Li, Li, Xu, Gao and Wang2019), as was found in the present study. However, nad5 had the fewest conservative sites in the mitochondrial genomes for 17 Ascaris samples (Zhou et al. Reference Zhou, Guo, Deng, He, Ouyang and Wu2020). These indicate that the least variable and most slowly evolving genes are the same in Ascaridomorpha; however, the least conserved gene is different.
Figure 3. (A) Sliding window analysis of the concatenated alignments of 12 PCGs of five Toxocara species. 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. (B) Evolutionary rates of 12 PCGs among five Toxocara species. The ratio of Ka/Ks is calculated for each PCG.
The Ka/Ks ratio is controlled by functionally related sequence contexts, such as encoding amino acids and participating in exon splicing. It is defined as the degree of evolutionary change; when the ratio is greater than 1, positive selection exists, indicating that non-synonymous mutations are more highly favored by Darwinian selection, and they will be retained at a rate greater than that of synonymous mutations (Liu et al. Reference Liu, Gasser, Otranto, Xu, Shen, Mohandas, Zhou and Zhu2013b; Xing et al. Reference Xing, Liang, Wang, Hu and Huang2022). Our Ka/Ks ratio data showed that all mt genes among Toxocara species were subject to purification selection and not a positive selection (Figure 3B). The results were the same as those of the evolutionary rates of PCGs between T. vitulorum and 28 other Rhabditida mitogenomes (Xie et al. Reference Xie, Wang, Chen, Wang, Zhu, Hu, Han, Wang, Zhou and Zuo2022).
Phylogenetic and genetic distance analysis
Complete mitochondrial sequences can provide a great source of species information and can provide new specific molecular markers for species taxonomy, population genetics, and systematics. For some taxonomic groups, genetic distance is the most effective model to quantify sequence divergences among individuals (Chagas et al. Reference Chagas, Ludwig, Pimentel, de Abreu, Nunez-Rodriguez, Leal and Kalapothakis2020). In the present study, the average interspecific genetic distances of 12 PCGs among five Toxocara species were found to range from 0.1544 (nad2) to 0.0259 (cox1). In a comparison with T. apodemi, the largest genetic distance (0.2064) was observed in the nad2 gene with T. vitulorum. The smallest genetic distance (0.0169) was found in the cox2 gene between T. apodemi and T. malaysiensis (Figure 4). Genetic distance analysis provided essential insights into the evolutionary mechanisms and patterns of the 12 PCGs in Toxocara mitochondrial genomes. The present results were larger than those of the 17 Ascaris samples (human-type Ascaris, pig-type Ascaris, and hybrid Ascaris) based on the 12 PCGs and the cox1 sequence reported in a previous study (Zhou et al. Reference Zhou, Guo, Deng, He, Ouyang and Wu2020). These findings indicate that the differences in the mtDNA sequence among the sequenced Ascaris individuals were smaller than those in Toxocara.
Figure 4. Genetic distance analysis of 12 PCGs among five Toxocara species. The genetic distances between T. apodemi and the other four Toxocara species are shown in the first four columns, and the average interspecific genetic distances of 12 PCGs among five Toxocara species are shown in the last column.
Phylogenetic analyses based on the concatenated amino acid sequences of 12 PCGs were used to assess the phylogenetic relationship of order for Ascaridomorpha with BI and ML methods (Figure 5, Figure S2). The two phylograms generated a similar topology, supporting the idea that each genus in the order forms a sister group with high statistical support, which was identical to that reported in a previous study (Xie et al. Reference Xie, Wang, Chen, Wang, Zhu, Hu, Han, Wang, Zhou and Zuo2022). In the branch of Toxocara, T. apodemi, T. canis, and T. cati separately formed a distinct branch and were always sister taxa to other congeneric species. T. malaysiensis and T. vitulorum clustered together, indicating that the two species were more closely related to each other than to other species in Toxocara. These results were inconsistent with findings of previous studies, in which T. malaysiensis and T. cati formed a sister group, and T. vitulorum formed a distinct branch (Li et al. Reference Li, Lin, Song, Wu and Zhu2008; Xie et al. Reference Xie, Wang, Chen, Wang, Zhu, Hu, Han, Wang, Zhou and Zuo2022; Xing et al. Reference Xing, Liang, Wang, Hu and Huang2022). Moreover, the phylogenetic analyses based on 18S showed that T. apodemi was more closely related to T. canis and Toxascaris leonina (Kim et al. Reference Kim, Hong, Ryu, Park, Cho, Yu, Chae, Choi and Park2020). In the present study, the mitochondrial genome sequences of T. apodemi strongly support it being a member of the Toxocara clade.
Figure 5. Phylogenetic analyses reconstructed using concatenated nucleotide sequences of 12 PCGs of complete mt genomes in 35 Ascaridomorpha species. The tree was performed by the BI method. Caenorhabditis elegans is an outgroup. T. apodemi in the current study is underlined.
Conclusions
In the present study, the complete mitogenome sequence of T. apodemi was determined and characterised. Comparative genomics suggested that T. apodemi was more closely related to T. canis in nucleotide and amino acid sequences. Phylogenetic analysis showed that T. apodemi was a member of Toxocara. These results contribute novel genetic markers for the phylogenetic and evolutionary study of Toxocaridae species.
Supplementary material
The supplementary material for this article can be found at http://doi.org/10.1017/S0022149X24000221.
Financial support
This work was supported by the National Natural Science Foundation of China (32302897), “Pioneer” and “Leading Goose” R&D Program of Zhejiang (2022C02031), Zhejiang Provincial Key R&D Program of China (2021C02049), Development Science and Technology Program in Taizhou (23sfa02), and Taizhou City Foundation for Talents (TZ2022-2).
Competing interest
None.
Ethical standard
Not applicable.
