Complete Mitochondrial Genome of Contracaecum sp. (Nematoda: Ascarididae) from night herons in China

Helminth specimens were obtained from the digestive tracts of gray and night herons in Beijing Zoo, China. The species were washed with ultrapure water and physiological saline solution, fixed in 75% ethanol, and stored at -40°C. The specimens were preliminarily identified as Contracaecum based on hosts and primary characteristic morphology (Zhang et al., 2021). For additional examination of molecules, the total genomic DNA was extracted using a QIAampÒ DNA Micro Kit as per the manufacturer’s instructions. Based on polymerase chain reaction (PCR) amplification of partial cox1 (with primers JB3 – JB4.5) (Bowles et al., 1992; Bowles and McManus, 1994) and ITS (including ITS-1, 5.8S, and ITS-2) (with primers NC5 – NC2) (Newton et al., 1998; Chilton et al., 2001), the worms were recognized at the species level. The obtained ITS sequence was totally matched with published Contracaecum sp. (GenBank no. MW538933~36), and the partial cox1 sequence showed 99.7% identity with Porrocaecum reticulatum (GenBank no. MF113244).

The genomic DNA sample was fragmented to a size of 350 bp. The DNA libraries were sequenced using high-throughput sequencing (HTS) on an Illumina Hiseq 6000 platform (Novogene Co. Ltd., Tianjin, China), and 250-bp paired-end reads were generated. The raw data were obtained and recorded in FASTQ format. Then, the reads with low-quality bases (Phred quality <5) or uncertain reads with repetitive “N” bases were filtered to acquire clean data. The partial cox1 sequence was used as the initial reference to assemble complete mt sequence of Contracaecum sp. using Geneious Prime 2022.0.1 (Kearse et al., 2012). The assembly was operated with the following parameters: (i) minimum overlap within the range of 150 bp to 200 bp; (ii) minimum overlap identity among 98% to 100%; and (iii) maximum gap of 5 bp. The assembled mitogenome was verified by long PCR with designed primers (Table S1 and Fig. S1 in Supplementary Materials).

ORF Finder (https://www.ncbi.nlm.nih.gov/orffinder/) was used to identify the start/stop codons and boundaries of protein-coding genes (PCGs). Later, two ribosomal RNAs (rRNAs, rrnL and rrnS) were framed using Tandem repeats finder (Benson, 1999). The 12 protein genes were then further confirmed with previously published Ascarididae sequence (Contracaecum osculatum, GenBank no. JN786330). The tRNAscan-SE 2.0 (Chan et al., 2021) with a cutoff score of 1.0 and MITOs (Bernt et al., 2013) were applied to search 22 potential transfer RNAs (tRNAs).

Based on available mitogenome sequences of the genus Contracaecum in the NCBI, mt sequences were aligned using Clustal X1.83 to a single alignment dataset, including C. osculatum, C. rudolphii, C. ogmorhini, and Contracaecum sp. The nucleotide diversity of Contracaecum species was computed by DnaSP v5 using sliding windows (Librado and Rozas, 2009). The parameters of the sliding window were followed with 300-bp window length and a default 25-bp step site to calculate the nucleotide diversity (Pi or p). Each boundary of protein genes was identified due to mid-point position, and we then graphed nucleotide diversity for 12 protein genes from Contracaecum.

A total of 41 mitogenomes of species from families Ascaridoidea and Heterakoidea were applied to analyze phylogeny with outgroups Enterobius vermicularis (GenBank accession no. EU281143) and Wellcomia siamensis (GenBank accession no. NC_016129) (Table S2 in Supplementary Material). Each amino acid sequence was aligned using a MAFFT computational algorithm (Katoh et al., 2019). The aligned sequences were then concatenated to a single alignment dataset. The ambiguous gaps in the alignment were excluded by Gblocks 0.91b with default parameters “less stringent” (Dereeper et al., 2008). Computational algorithm maximum likelihood (ML) (Guindon et al., 2010) was conducted to perform a phylogenetic tree with the best model “JTT+I+G+F” screened by ProtTest 3.4.2 (Darriba et al., 2011) and applied 1,000 replicates. Bayesian analysis was operated with MrBayes 3.2 (Ronquist et al., 2012), and “GTR + F + G” was selected as the most suitable model by ModelFinder in IQTree v.2.1.3 (Kalyaanamoorthy et al., 2017). Four Markov chains were progressed with 1,000,000 MCMC generations, with sampling analysis tree each 100 generations. The residual trees were calculated with Bayesian posterior probabilities (BPP), burning first 250 trees.

The clean data of Contracaecum sp. are nearly 2 GB with a total of 8,677,194 × 2 clean reads for further assembling. The circular mt genome of Contracaecum sp. (GenBank accession: ON149889) assembled was 14,082 bp in size, shorter than that Zhang et al. (2021) published, with 12 PCGs, 22 tRNAs, 2 rRNAs, and 2 noncoding regions (NCRs) (Table 1 and Fig. 1). A total of 36 genes were transcribed in the forward direction and gene arrangement was recognized as the typical GA3 pattern, which is mostly observed in the worms (Liu et al., 2013). Consistent with previous reports, there was an obvious bias of A + T bases (71.2%). A total of 10 intergenic regions were found among the complete mt genome of Contracaecum sp. ranging from 1 bp to 16 bp (Table 1). One short NCR (122 bp) was located between nad4 and cox1, and one long NCR (691 bp) was placed in tRNA-Ser2 and tRNA-Asn. The values of AT skew were negative from -0.475 (nad6) to -0.111 (NCRs), and inversely, the values of GC-skew were positive with scope 0.226 (nad4) to 0.674 (nad3), suggesting Ts and Gs were more frequently used in the genome.

Figure 1

TTG was the most common initial codon in this study, followed by ATT. TTG was used as the start codon for nine genes (cox1-3, cytb, nad1-4, and nad6) (Table 1). The rest three PCGs (atp6, nad4L, and nad5) used ATT as the initial codon. Generally, TAG and TAA were common stop codons in metazoans (Hu et al., 2004). In this study, TAA was the most frequent termination among nad6, nad4L, nad4, cytb, and nad2. The genes nad1 and nad3 used TAG as their stop codon. The rest genes used incomplete stop codons T (cox1 and cox3) or TA (atp6, cox2 and nad5), respectively.

A total of 3,422 amino acids were translated by 12 PCGs. TTT (480) was the most common codon used in encoding Phe, followed by GTT (219, Val), TTG (216, Leu), and ATT (214, Ile). Leu (519) and Phe (499) were the most frequent amino acids, while Arg (34) was the least. There was a tendency of Gs and Ts in the same amino acid by comparing the relative synonymous codon usage (RSCU) (Table 2). The AT content of 12 protein genes ranged from 66.7% (cox1) to 78.9% (nad6) (Table 3). The values of AT skew ranged from –0.475 (nad6) to –0.373 (cox2), while the values of GC skew were 0.226 (nad4) to 0.674 (nad3), suggesting the bias of T and G bases.

The length of 22 tRNAs ranged from 51 bp (tRNA-Ser1) to 65 bp (tRNA-His). The typical structure of tRNA consisted of one acceptor stem, a dihydrouridine loop (D-loop), an anticodon loop, TYC loop, and related arms fixing with them (Su et al., 2020). However, the TψC loop was always replaced by a TV replacement loop in nematodes. In our study, 16 of 20 tRNAs (excluding tRNA-Ser1 and tRNA-Ser2) lacked a TYC loop, replaced by several nucleotide residues, which compromised the TV replacement loop (Hu et al., 2004). The tRNA-His, tRNA-Ile, and tRNA-Met were observed in a relatively standard cloverleaf structure with a TYC loop, although the latter two (tRNA-Ile and tRNA-Met) lacked DHU stem. The tRNA-Ser1 and tRNA-Ser2 were similar to previous reports with one TYC-loop but lacked D-loop (Su et al., 2020).

Ribosomal RNAs of Contracaecum sp. were fixed as a GA3 pattern. The rrnL was located between tRNA-His and nad3 with a size of 959 bp, and the rrnS gene was located between tRNA-Glu and tRNA-Ser2 with a size of 711 bp (Table 1). The content of A + T for rrnL and rrnS was 75.6% and 70.6%, respectively. There were two NCRs among the mt genome of Contracaecum sp. One short region was placed in nad4 and cox1 with a length of 122 bp, and the long region was situated between tRNA-Ser2 and tRNA-Asn with a length of 691 bp.

Based on aligned nucleotide sequences among species C. osculatum, C. rudolphii, C. ogmorhini, and Contracaecum sp., nucleotide diversities (Pi) were calculated based on the sliding window. The values of Pi ranged from 0.124 to 0.181 by analyzing a window of 300 bp and a default step of 25 bp (Fig. 2). The most variable genes were cytb (0.178), nad2 (0.181), nad4 (0.179), and nad6 (0.172), and the most conserved genes were cox1 (0.124) and cox2 (0.130) in Contracaecum (Fig. 2). Protein genes cox1 and cox2 seemed to be the most stable genes in Contracaecum nematodes with the least variation, which could be used as molecular markers to identify species from Contracaecum. Results also supported that nad2 and nad4 could act as alternative markers among nematodes isolated from different distributions.

Figure 2

The present phylogenetic trees were constructed based on the 12 PCGs of 41 available mt genome sequences from the superfamilies Ascaridoidea and Heterakoidea (Table S2 in Supplementary Material). Two phylogenetic trees, both Bayesian inference (BI) and ML, had similar topologies, excluding species within the superfamily Heterakoidea. The topologies of ML and BI phylogenetic trees were highly similar to those of previous studies (Liu et al., 2016; Zhang et al., 2021; Zhao et al., 2021). Contracaecum sp. formed a branch with Contracaecum nematodes, indicating a closer relationship within the genus with strong support (Fig. 3); however, the long distance between Contracaecum sp. and the other three Contracaecum species (C. osculatum, C. rudolphii, and C. ogmorhini) was longer than the branch distance within other anisakid nematodes, which further indicated Contracaecum sp. was a novel species and verified the hypothesis of Zhang et al. (2021) proposed. According to the structure of phylogenetic trees, results supported that the superfamilies Ascaridoidea and Heterakoidea were monophyletic and evidenced families, including Ascarididae, Anisakidae, Heterocheiidae, Toxocaridae, and Cucullanidae, were monophyletic, consistent with previous studies (Li et al., 2018; Zhao et al., 2021).

Figure 3

Within the superfamily Ascaridoidea, both BI and ML showed identical topologies. Among the family Ascarididae, the genera Ascaris, Baylisascaris, Toxascaris, and Parascaris had a closer relationship than Ophidascaris, which was similar to Zhou et al. (2021) reported. Based on morphological descriptions, the genus Ophidascaris was classified as a member of the superfamily Ascaridoidea (Pinto et al., 2010), and phylogenetic analyses suggested the genus Ophidascaris was more related to the family Ascaridae. However, the distance between the genus Ophidascaris and other Ascaridae genera was longer, suggesting there was systematic controversy in the Ascaridae. In the present study, the family Ascarididae was closely related to the family Anisakidae (Fig. 3), different from the previous study where the family Ascarididae was closely related to Toxocaridae (Zhou et al., 2021). In addition, results also supported the monophyly of all 5 families and all 11 genera within the superfamily Ascaridoidea with strong support (BPP = 1, Bf >70, Fig. 3), consistent with records (Liu et al., 2016; Zhao et al., 2018).

Liu et al. (2013) confirmed that Ascaridia columbae was more related to Ascaridia sp. than A. galli. The phylogenetic analyses in the present study also confirmed this. In ML analysis, the topology showed that A. galli was more related to Heterakis species with high statistical support (Bf = 87), in line with Liu et al. (2016) studied. However, BI analysis presented a totally different topology from that of ML tree. A. galli formed a distinct branch from genera Heterakis and Ascaridia with strong support (BPP = 1), hypothesizing A. galli might be another genus. Results also showed the family Heterakidae was a sister taxon to Ascaridiidae, and phylogenetic analyses (BI and ML) suggested the family Ascaridiidae might be paraphyly.