Introduction
The phylum Nematoda is one of the most successful groups of organisms on Earth. Although their morphology may appear simplistic, these animals are extremely diverse and abundant, being present in nearly every habitat (Dorris et al. Reference Dorris, De Ley and Blaxter1999; Poinar Reference Poinar2009). Many lineages are parasitic, affecting a wide range of species including those of agricultural, human, and veterinary interest (Brooker Reference Brooker2010; Taylor et al. Reference Taylor, Coop and Wall2007) leading to severe societal and monetary implications. Despite this importance, much of the evolutionary history of nematodes remains unresolved. This limitation partially stems from the availability of few taxonomically informative morphological characters, due in part to their small size and relative structural simplicity, and the scarcity of fossil records (Dorris et al., Reference Dorris, De Ley and Blaxter1999; Poinar, Reference Poinar2003; Blaxter et al., Reference Blaxter, Koutsovoulos, Jones, Kumar, Elsworth, Olson, Hughes and Cotton2016), but also from the uneven representation of nematode diversity in molecular studies. Most sequencing efforts have focused on species of medical or economic relevance, leaving the vast diversity of free-living and wildlife-infecting nematodes underrepresented. This taxonomic imbalance restricts our ability to reconstruct a comprehensive picture of nematode evolution.
The first phylum-wide phylogenetic assessments were performed using the 18S ribosomal DNA genes, and led to substantial revisions of nematode taxonomy (De Ley and Blaxter Reference De Ley, Blaxter and Lee2002, Reference De Ley and Blaxter2004), including the proposal of new clades. More recently, additional mitochondrial genes have also been used, both for broader phylogenetics and population assessments, some studies even using complete mitogenomes (Gilabert and Wasmuth Reference Gilabert and Wasmuth2013; Kang et al. Reference Kang, Sultana, Eom, Park, Soonthornpong, Nadler and Park2009; Kim et al. Reference Kim, Eom and Park2006). Mitochondrial sequences are widely used to infer phylogenetic relationships at different evolutionary levels due to various characteristics of mitogenomes, including their high evolutionary rate, reduced recombination, relatively constant gene content and order, haploidy, maternal transmission (including in nematodes, Anderson et al. Reference Anderson, Komuniecki, Komuniecki and Jaenike1995; Tsang and Lemire Reference Tsang and Lemire2003), and general ease of sequencing (Jex et al. Reference Jex, Hall, Littlewood and Gasser2010). To date, over 600 complete mitochondrial genome entries are available on GenBank for the phylum Nematoda (https://www.ncbi.nlm.nih.gov), which still represents a very small fraction of the huge diversity of these organisms. Moreover, sequencing efforts have been biased toward nematode species of medical, veterinary, or agricultural importance or linked to model organisms. Our understanding of nematode evolutionary history would benefit from sequencing a wider variety of organisms (Otto et al. Reference Otto, Gilabert, Crellen, Böhme, Arnathau, Sanders, Oyola, Okouga, Boundenga and Willaume2018; Woo et al. Reference Woo, Ansari, Otto, Klinger, Kolisko, Michálek, Saxena, Shanmugam, Tayyrov and Veluchamy2015), for example, allowing for more effective comparative studies, including on adaptation, genome structure, gene expression, and protein function, host-switching, and origins of parasitism. For this, the development of high-throughput sequencing (HTS) techniques has become an essential asset. HTS has substantially decreased the cost of DNA sequencing, and obtaining new genomes has become a much less arduous and expensive task (Goodwin et al. Reference Goodwin, McPherson and McCombie2016). As a result, the number of mitogenomes being sequenced has grown exponentially (Jex et al. Reference Jex, Hall, Littlewood and Gasser2010), which in turn is allowing for an increased resolution in deciphering their diversification. Mitochondrial genome assessments have identified niche-specific selection dynamics within gastrointestinal parasites of the superfamily Strongyloidea (Li et al. Reference Li, Zhao, Xu, Lai, Gong, Wang, Zhang, Chen, Guo, Wang and Liu2025), but data for representatives of many orders are still unavailable, limiting comparative studies to determine aspects such as host-switching or mitochondrial gene rearrangements.
Pinworms of the superfamily Oxyuroidea are one of the most notable groups of nematodes, and probably among the first parasites known to man. They are haplodiploid intestinal parasites with a direct life-cycle and oral–faecal transmission (Adamson Reference Adamson1989). The family Pharyngodonidae is ubiquitous and diverse, comprising 25 genera (Adamson Reference Adamson1989; Bursey et al. Reference Bursey, Goldberg and Kraus2009) that infect all classes of vertebrates except for birds (Anderson et al. Reference Anderson, Chabaud and Willmott2009). However, the evolutionary history of the Pharyngodonidae remains unclear. Currently, the availability of genetic data on this family is restricted to a few lizard-infecting species belonging to a handful of genera (Falk and Perkins Reference Falk and Perkins2013; Jorge et al. Reference Jorge, Perera, Roca, Carretero, Harris and Poulin2014; Mockett et al. Reference Mockett, Bell, Poulin and Jorge2017). This limited data has revealed cases of taxonomic incongruences, cryptic species, and undescribed diversity (Falk and Perkins Reference Falk and Perkins2013; Jorge et al. Reference Jorge, Perera, Carretero, Harris and Roca2013, Reference Jorge, Perera, Poulin, Roca and Carretero2018; De Sousa et al. Reference De Sousa, Jorge, Carretero, Harris, Roca and Perera2018) suggesting that increasing data generation and availability would be beneficial for further research on the biology of these organisms. Currently, only the 18S rRNA, 28S rRNA, and Cox1 genes for a small number of species have been sequenced. In this study, with the objective of developing new molecular markers, we have reconstructed and annotated four mitochondrial genomes from species of the genera Spauligodon, Pharyngodon, Parapharyngodon, and Thelandros. We combined the new data from representatives of the family Pharyngodonidae with available mitogenomes from two other families, Oxyuridae and Heteroxynematidae, to assess phylogenetic relationships across the superfamily Oxyuroidea. These mitogenomes provide the groundwork for future studies to determine patterns of evolution and adaptations to hosts across a broader taxonomic range of these important parasites.
Methodology
Nematodes were collected from faecal samples of wild lizards. Thelandros galloti from host Gallotia caesaris (50 pooled individuals), and Spauligodon atlanticus and Parapharyngodon echinatus, both from Gallotia atlantica (113 and 28 individuals respectively), were all collected in the Canary Islands. Nematodes from the genus Pharyngodon were collected from Chalcides ocellatus (80 individuals) in Morocco. Nematode individuals were isolated from faecal material under an Olympus SZX2 stereomicroscope. Given their small size and resulting low amount of total DNA, pools of nematode individuals were produced in order to obtain at least 1000 ng of DNA. In order to minimize the effect of this pooling approach, pools consisted of nematodes from the same host individual or, when this was not possible, infecting the same host species. A 600 bp fragment of the Cox1 gene was sequenced for all pools in order to confirm the morphological identification and ensure limited genetic variability across individuals for subsequent HTS sequencing. Metazoan invertebrate universal primers were used (Folmer et al. Reference Folmer, Black, Hoeh, Lutz and Vrijenhoek1994), PCR conditions were the same as Jorge et al. (Reference Jorge, Roca, Perera, Harris and Carretero2011), and Sanger sequencing was performed by an external company (Beckman, UK).
Total DNA was extracted from the pools using the Qiagen QIAamp DNA Micro kit (Valencia, California, USA), following the manufacturer’s instructions. Up to 2000 ng of gDNA was sheared to an average size of ~400 bp using the Q800R sonicator (QSonica, Newton, Connecticut, USA) and subsequently prepared for Illumina sequencing using Nextera-style adapters with dual-indexed 8 bp barcodes. Samples were pooled together for paired-end sequencing (read length 2 × 150 bp) on an Illumina HiSeq 2500 platform, performed by an external company (Sanford Burnham Prebys, USA).
General sequence quality was checked on FastQC (Andrews Reference Andrews2010). Adapter removal and quality trimming were performed using fastp 0.24 (Chen et al. Reference Chen, Zhou, Chen and Gu2018), using default settings. For the mitogenome assemblies, two different approaches were explored. First, a reference-based approach, using the Bowtie v2.3 read mapper (Langmead and Salzberg Reference Langmead and Salzberg2012) on loose stringency settings, was applied where reads were mapped to six mitochondrial genomes of other Oxyurida species from the family Oxyuridae (Aspiculuris tetraptera (KT764937), Enterobius vermicularis (EU281143), Oxyuris equi (KP404095), Passalurus ambiguus (KT879302), Syphacia obvelata (KT900946), and Wellcomia siamensis (GQ332427)). Mapped sequences were then assembled into contigs using SOAPdenovo v2.04 (Luo et al. Reference Luo, Liu, Xie, Li, Huang, Yuan, He, Chen, Pan and Liu2012). Secondly, a baiting approach using the programs NOVOPlasty v4.3.5 (Dierckxsens et al. Reference Dierckxsens, Mardulyn and Smits2016) and MITObim v1.9.1 (Hahn et al. Reference Hahn, Bachmann and Chevreux2013) was used. The Cox1 sequences previously obtained for the samples were used as starting points (baits) for the consecutive runs of iterative mapping employed by both programs. Additionally, for MITObim, forward and reverse reads were paired using PEAR v0.9.6 (Zhang et al. Reference Zhang, Kobert, Flouri and Stamatakis2013), keeping only reads longer than 250 bp to fix repetitive parts that caused contig extension to halt. Afterwards, the obtained assemblies were annotated for protein-coding, tRNA and rRNA genes using a combination of MitoZ 3.6 (Meng et al. Reference Meng, Li, Yang and Liu2019) and MITOS2 (Bernt et al. Reference Bernt, Donath, Jühling, Externbrink, Florentz, Fritzsch, Pütz, Middendorf and Stadler2013), with remaining tRNAs found using either MiTFi (Jühling et al. Reference Jühling, Pütz, Bernt, Donath, Middendorf, Florentz and Stadler2012) or ARWEN (Laslett et al. Reference Laslett and Canbäck2008).
To estimate phylogenetic relationships, 13 additional mitogenomes were downloaded from GenBank, including all available representatives from the families Oxyuridae (8 mitogenomes) and Heteroxynematidae (2 mitogenomes). Three additional mitogenomes, from Thelazia calipaeda, Metathelazia capsulata, and Toxocara canis, were included as outgroups. CDS of 12 protein-coding genes were analysed. These were aligned using MAFFT (Katoh and Standley Reference Katoh and Standley2013), with ambiguously aligned regions removed using Gblocks (Castresana Reference Castresana2000). Our final dataset consisted of 17 taxa, with 8937 base pairs. Bayesian inference phylogenies were inferred using MrBayes v3.2.7 (Ronquist et al. Reference Ronquist, Teslenk, Mark, Van Der Ayers, Darling, Höhna, Larget, Liu, Suchard and Huelsenbeck2012), under a partition model for 10 million generations, and with the initial 25% of sampled data discarded as burnin.
Results and discussion
Mapping approaches were ineffective in obtaining the mitogenomes of Pharyngodonidae nematodes, as they produced fragmented and incomplete reconstructions. Baiting approaches were more successful, producing a linear complete sequence or long fragments that could be aligned together into a complete molecule. The ease and efficiency of reconstruction varied, with samples containing repetitive regions proving more difficult to form into a circle. For instance, in the case of Thelandros, multiple TATA and poly-T repeats present would repeatedly halt contig extension at the same location and several runs with different starting points were needed to produce a circular contig. For all samples, a manually curated consensus of the assemblies produced by MITObim and NOVOPlasty is shown, with remapping data used to solve small discrepancies. Regarding annotation, MitoZ performed well with protein-coding genes but was missing several tRNA annotations, which were subsequently added using MITOS2, MiTFi, and ARWEN. Therefore, the results shown are a combination of all annotations.
The four resulting molecules (GenBank accession numbers PX608396-PX608399) had genome sizes of 13,692 bp for Spauligodon, 13,706 bp for Pharyngodon, 14,672 bp for Thelandros, and 16,700 bp for Parapharyngodon (Fig. 1). The disparate size of the mitogenome from Parapharyngodon was due to the existence of two large non-coding regions (with 1,202 and 1,394 bp). The overall nucleotide composition was very AT rich: 68% for Pharyngodon, 70.9% for Spauligodon, 71.4% for Parapharyngodon, and 70.5% for Thelandros. All four mitochondrial genomes contained 12 protein-coding genes (Cox1 to Cox3, Cyt b, NAD1 to NAD6, NAD4L, and ATP6), 2 rRNA genes (12S and 16S), and 22 tRNA genes (Fig. 1). All these genes were encoded on the same strand. Similar to most other nematode genomes, the ATP8 gene was not present. The mitogenomes of Spauligodon and Pharyngodon were not only similar in size, but also presented the same gene order. In contrast, Thelandros and Parapharyngodon showed a broadly similar gene arrangement to each other with exception of the position of the tRNA-E, yet both genera were highly divergent in gene order from Spauligodon and Pharyngodon. Several non-coding regions were present across the mitogenomes, ranging in size from just a couple of base pairs to 1,394 bp. Notably, the large non-coding regions present in Thelandros and Parapharyngodon occupy the same genomic location where Spauligodon and Pharyngodon contain the 16S, NAD6, and Cytb genes; however, alignment and BLAST analyses indicate that these non-coding regions bear no significant sequence similarity to that gene-block.
Mitochondrial genomes of the four pharyngodonid nematodes: Spauligodon, Pharyngodon, Thelandros, and Parapharyngodon. The 12 protein-coding genes are shown in green, the 2 rRNA genes in red, and the 22 tRNA genes in pink. Sequence length is also given under each genus name.
In this work, we present the first mitochondrial genomes of the nematode family Pharyngodonidae, from the genera Thelandros, Spauligodon, Pharyngodon, and Parapharyngodon. Methodologically, two very different approaches were used to assemble these mitogenomes: mapping to reference genomes and baiting followed by iterative mapping. The latter was found to be more successful, which was not unexpected given the high variability of mitogenome structure in nematodes (Yatawara et al. Reference Yatawara, Wickramasinghe, Rajapakse and Agatsuma2010). Together with the short read size, capturing mitochondrial sequences from whole genome data can be quite difficult, even at low stringency settings (Yuan et al. Reference Yuan, Xia, Zheng and Zeng2016). When available reference sequences are limited or from distantly related organisms, baiting approaches for mitogenome assembly provide better results, as was the case in this study. Their only limitation here is the requirement of a closely related ‘bait’ sequence to start the reconstruction. In our case, we used the Cox1 gene initially obtained by PCR and Sanger sequencing for all four samples, although initial mapping to distantly related mitochondrial genomes can also identify conserved areas to be used as initial bait sequence. Baiting and iterative mapping provide an easy, fast, and computationally light way of reconstructing mitochondrial and other circular genomes from complex whole-genome sequencing data. Additionally, we recommend testing different algorithms and settings as in our study the two softwares used, MITObim and NOVOPlasty, produced slightly different results.
Metazoan mitochondrial genomes are considered to be relatively uniform in size and gene content across animals (Lang et al. Reference Lang, Gray and Burger1999): a 14–18 kb circular maternally transmitted double-stranded DNA molecule encoding 37 genes. The arrangement of these genes is usually relatively conserved, with some blocks of genes being shared even among different phyla (Boore Reference Boore1999). However, nematode mitogenomes deviate from several of these norms (reviewed in Kern et al. Reference Kern, Kim and Park2020). Molecule size varies considerably from 13 to 24 kb (Hu and Gasser Reference Hu and Gasser2006); the rRNA and tRNA genes are smaller, the ATP8 protein-coding gene is not present in most cases, among other differences and unique organizations (Gibson et al. Reference Gibson, Blok, Phillips, Hong, Kumarasinghe, Riley and Dowton2007; Lavrov and Brown Reference Lavrov and Brown2001). Considering this, the characteristics of the four mitogenomes here reconstructed were in accordance with those found in other families of the nematode class Secernentea. The four mitogenomes were within the normal size range, included 2 rRNA, 22 tRNA, and 12 protein-coding genes, and the ATP8 was not present. Likewise, all genes were located in the same strand and transcribed in the same direction which is the case for almost all chromadorean nematodes (Kim et al. Reference Kim, Kern, Kim, Sim, Kim and Kim2017), and base content was highly AT biased (Yatawara et al. Reference Yatawara, Wickramasinghe, Rajapakse and Agatsuma2010). Similar to other nematode mitogenomes, we also found extensive and multiple non-coding regions (e.g., Lavrov and Brown Reference Lavrov and Brown2001; Besnard et al. Reference Besnard, Jühling, Chapuis, Zedane, Lhuillier, Mateille and Bellafiore2014; Kikuchi et al. Reference Kikuchi, Afrin and Yoshida2016, Zhao et al. Reference Zhao, Tu, Wu and Li2018; Palomares-Rius et al. Reference Palomares-Rius, Cantalapiedra-Navarrete, Archidona-Yuste, Blok and Castillo2017; Huo et al. Reference Huo, Bao, Long, Qin and Zhang2023), which explains most of the size variation between the four mitochondrial genomes. Most non-coding regions were small in size and characterized as intergenic ‘spacers’, but for the longer non-coding regions, there were also sites of TATA tandem repeats. Such repeats are commonly found in the putative control region of metazoan mitogenomes and, as they potentially form a stem-loop secondary structure, they may function as a replication and/or transcription origin (Okimoto et al. Reference Okimoto, Chamberlin, Macfarlane and Wolstenholme1991, Reference Okimoto, Macfarlane, Clary and Wolstenholme1992; Wolstenholme Reference Wolstenholme1992).
In our estimate of relationships, representatives of the same species and genera all grouped as expected with strong support, as did Spauligodon with Pharyngodon, and Thelandros with Parapharyngodon (Fig. 2). However, neither of the families Pharyngodonidae nor Oxyuridae was monophyletic. Previous estimates of phylogenetic relationships based on mitogenomes within the superfamily Oxyuroidea have differed not only between them but also estimates based on maximum likelihood and Bayesian approaches have been incongruent (e.g. Huo et al. Reference Huo, Bao, Long, Qin and Zhang2023). Within the data there is notable base composition heterogeneity, as well as high genetic distances – maximum uncorrected distances reached 41.6%. Mitogenome-derived phylogenies have already been shown to be inconsistent with phylogenies derived from genomic data in certain groups of nematodes (Ahmed et al. Reference Ahmed, Roberts, Adediran, Smythe, Kocot and Holovachov2022). Additional assessments of the influence of mitochondrial gene rearrangements and taxon sampling will be needed to shed light on the causes of these topological discrepancies.
Estimate of phylogenetic relationships based on a Bayesian analysis of aligned protein-coding genes (8937bp) of the mitogenome. Bayesian Posterior Probability (BPP) support values are indicated beside each node. Inset are images of the four newly sequenced species, A) Parapharyngodon echinatus, B) Thelandros galloti, C) Pharyngodon sp., and D) Spauligodon atlanticus (all photographs by FJ).
Among the four reconstructed mitogenomes, there were clear similarities and major differences in the gene organization. Spauligodon and Pharyngodon presented the same gene order, while for Thelandros and Parapharyngodon only one tRNA varied in its location. Between these two pairs, however, there were extensive rearrangements of rRNA, tRNA, and protein-coding genes. This pattern likely reflects the taxonomic affinities of these genera, as each pair belongs to distinct subfamilies: Pharyngodoninae and Thelandroinae (Adamson Reference Adamson1989; De Sousa et al. Reference De Sousa, Jorge, Carretero, Harris, Roca and Perera2018; Pereira et al. Reference Pereira, Luque and Tavares2018). Such major gene rearrangements are also common in nematodes (Hu and Gasser Reference Hu and Gasser2006), even among closely related species (Kikuchi et al. Reference Kikuchi, Afrin and Yoshida2016; Palomares-Rius et al. Reference Palomares-Rius, Cantalapiedra-Navarrete, Archidona-Yuste, Blok and Castillo2017), and may constitute an additional mechanism of diversification. It is still not clear what regulates gene rearrangements in animal mitogenomes, but several models have been suggested, such as the ‘duplication and random loss model’ (Macey et al. Reference Macey, Larson, Ananjeva, Fang and Papenfuss1997) and the ‘intra-mitochondrial recombination model’ (Rokas et al. Reference Rokas, Ladoukakis and Zouros2003). Additionally, heteroplasmy (which has been reported in nematodes; Lunt and Hyman Reference Lunt and Hyman1997; Tsang and Lemire Reference Tsang and Lemire2003), the same direction of transcription for all genes, and the presence of multiple non-coding regions might further facilitate rearrangements. Other animal groups are notable for mitochondrial gene rearrangements. These include amphibians, for which gene rearrangements linked to the origin of light-strand replication were identified (Fonseca and Harris Reference Fonseca and Harris2008). Furthermore, these rearrangements help explain the difficulty in development of new genetic tools for non-model nematodes and why (along with sequence divergence) universal primers may not work for some species. Increasing taxon sampling will also be necessary to further compare mitochondrial genomes to help identify evolutionary patterns for the most varied processes, and to better understand the origin of the impressive biodiversity of nematodes. For parasitic nematodes, this can be of particular importance for evolutionary and health research, as it will shed light on the emergence of life strategies, mechanisms of infection, immunology, metabolism, and development of vaccines and drug treatments (Coghlan et al. Reference Coghlan, Tyagi, Cotton, Holroyd, Rosa, Tsai, Laetsch, Beech, Day and Hallsworth-Pepin2019).
Acknowledgements
BT, KPM, AP, and JA were funded by the Fundação para a Ciência e Tecnologia (FCT) through contracts PD/BD/52601/2014, PD/BD/52604/2014, IF/01257/2012, and PTDC/BIA-EVL/29115/2017, respectively. KPM was supported by FWO post-doctoral fellowships under grant numbers 1224223N and 12AFD26N. The study was funded by the project IF01257/2012/CP0159/CT0005 (to AP) and ‘Genomics and Evolutionary Biology’, the latter cofinanced by North Portugal Regional Operational Programme 2007/2013 (ON.2 – O Novo Norte) under the national strategic Reference Framework (NSRF) and the European Regional Development Fund (ERDF). Special thanks to F. Oxelfelt for his help in organizing some of the funding. The authors declare no conflicts of interest in the publication of this manuscript.
