Phylomitogenomic Analyses Provided Further Evidence for the Resurrection of the Family Pseudoacanthocephalidae (Acanthocephala: Echinorhynchida)

Simple Summary Acanthocephalans, commonly known as spiny-headed or thorny-headed worms, are a small group of endoparasites with veterinary, medical and economic importance due to their ability to cause disease in domestic animals, wildlife, and humans. In recent decades, great progress has been made using mitochondrial genome data to clarify the phylogenetic relationships of acanthocephalans. However, the current mitochondrial genome database for acanthocephalans remains very limited. Herein, the characterization of the mitochondrial genome of Pseudoacanthocephalus bufonis (Shipley, 1903), the first representative of the family Pseudoacanthocephalidae, is reported. Phylogenetic analyses using the amino acid sequences of 12 protein-coding genes supported the validity of the family Pseudoacanthocephalidae and suggested a close affinity between Pseudoacanthocephalidae and Cavisomatidae. Our phylogenetic results also showed that the families Polymorphidae and Centrorhynchidae have a closer relationship than Plagiorhynchidae in the Polymorphida. These findings contribute to revealing the patterns of mitogenomic evolution in this group and represent a substantial step towards reconstructing the classification of the phylum Acanthocephala. Abstract The phylum Acanthocephala is an important monophyletic group of parasites, with adults parasitic in the digestive tracts of all major vertebrate groups. Acanthocephalans are of veterinary, medical, and economic importance due to their ability to cause disease in domestic animals, wildlife, and humans. However, the current genetic data for acanthocephalans are sparse, both in terms of the proportion of taxa surveyed and the number of genes sequenced. Consequently, the basic molecular phylogenetic framework for the phylum is still incomplete. In the present study, we reported the first complete mitochondrial genome from a representative of the family Pseudoacanthocephalidae Petrochenko, 1956. The mitogenome of Pseudoacanthocephalus bufonis (Shipley, 1903) is 14,056 bp in length, contains 36 genes (12 protein-coding genes (PCGs) (lacking atp8), 22 tRNA genes, and 2 rRNA genes (rrnL and rrnS)) and two non-coding regions (NCR1 and NCR2), and displayed the highest GC-skew in the order Echinorhynchida. Phylogenetic results of maximum likelihood (ML) and Bayesian inference (BI) using the amino acid sequences of 12 protein-coding genes in different models provided further evidence for the resurrection of the family Pseudoacanthocephalidae and also supported that the order Echinorhynchida is paraphyletic. A monophyletic clade comprising P. bufonis and Cavisoma magnum suggests a close affinity between Pseudoacanthocephalidae and Cavisomatidae. Our phylogenetic analyses also showed that Polymorphidae has a closer relationship with Centrorhynchidae than Plagiorhynchidae in the monophyletic order Polymorphida.


Introduction
Acanthocephala is an important group of obligate endoparasites, with more than 1300 species parasitizing the digestive tracts of all major lineages of vertebrates and their larvae developing in arthropods [1][2][3][4]. According to the current classifications based on a combination of morphological and ecological traits, the phylum is divided into three classes, Archiacanthocephala, Eoacanthocephala, and Palaeacanthocephala, which include 10 orders, 26 families, and over 160 genera [1,[5][6][7].
In order to further test the monophyly of the orders Echinorhynchida and Polymorphida, assess the validity of the recently resurrected family Pseudoacanthocephalidae, and clarify the evolutionary relationships of the Pseudoacanthocephalidae and the other families in Palaeacanthocephala using mitogenome data, the complete mitochondrial genome of Pseudoacanthocephalus bufonis, the first mitogenome from the Pseudoacanthocephalidae, was sequenced and annotated for the first time. Moreover, phylogenetic analyses of the protein-coding genes of all available acanthocephalan mitogenomes were performed using maximum likelihood (ML) and Bayesian inference (BI) in different models.

Parasite Collection and Species Identification
A total of 14 spot-legged tree frogs, Polypedates megacephalus Hallowell (Anura: Rhacophoridae), were caught by hand at night in the Diaoluo Mountains, Hainan Island, China, and euthanized by injection of an overdose of pentobarbitone sodium solution. The acanthocephalan specimens were collected from the intestine of the host. For light microscopical studies, acanthocephalans were cleared in glycerine. Photomicrographs were recorded using a Nikon ® digital camera coupled to a Nikon ® optical microscopy (Nikon ECLIPSE Ni-U, Nikon Corporation, Tokyo, Japan). The specimens were identified as P. bufonis based on morphological features reported in previous studies [31][32][33][34][35][36]. The terminology is according to the previous study [37]. Voucher specimens were deposited in the College of Life Sciences, Hebei Normal University, Hebei Province, Shijiazhuang, China (HBNU-A-2022A001L).

Molecular Procedures
For molecular analysis, the genomic DNA was extracted using a modified CTAB (pH 8.0)-based DNA extraction protocol as described in Zhao et al. [38]. The genomic DNA library was constructed, and a total of 20 GB of clean data were generated using the The complete mitochondrial genome was assembled using GetOrganelle v1.7.2a [39]. Protein coding genes (PCGs), rRNAs, and tRNAs were annotated using the MitoS web server (http://mitos2.bioinf.uni-leipzig.de/index.py, accessed on 20 January 2022) and MitoZ v2.4 [40]. The open reading frame (ORF) of each PCG was confirmed manually by the web version of ORF finder (https://www.ncbi.nlm.nih.gov/orffinder/, accessed on 10 March 2022). The "lost" tRNA genes ignored by both MitoS and MitoZ were identified using BLAST based on a database of the existing tRNA sequences of Acanthocephala. The secondary structures of tRNAs were predicted by the ViennaRNA module [41], building on MitoS2 [42] and RNAstructure v6.3 [43], followed by a manual correction. MitoZ v2.4 was used to visualize and depict gene element features [40]. The base composition, amino acid usage, and relative synonymous codon usage (RSCU) were calculated by a Python script, which refers to codon adaptation index (CAI) [44]. The total length of the base composition included ambiguous bases. The base skew analysis was used to describe the base composition of nucleotide sequences. The relative values were calculated using the formulas: ATskew = A−T A+T and GCskew = G−C G+C . The complete mitochondrial genome sequence of P. bufonis obtained herein was deposited in the GenBank database (http://www.ncbi.nlm.nih.gov, accessed on 5 October 2022).

Phylogenetic Analyses
Phylogenetic analyses were performed based on concatenated amino acid sequences of 12 PCGs using maximum likelihood (ML) and Bayesian inference (BI). Gnathostomula armata and G. paradoxa (Gnathostomulida) were chosen as the out-group. The in-group included 7 species of rotifers and 24 species of acanthocephalans. Detailed information on representatives included in the present phylogeny was provided in Table 1. The phylogenetic trees were re-rooted on Gnathostomulida. Genes were aligned separately using MAFFT v7.313 under the iterative refinement method of E-INS-I [45]. Ambiguous sites and poorly aligned positions were pruned using BMGE v1.12 (m = BLOSUM90, h = 0.5) [46]. The aligned and pruned sequences were concatenated into a matrix by PhyloSuite v1.2.2 [47]. The pruned alignments were then concatenated into the "AA" matrix with the amino acid sequences of PCGs (2087 sites). Bayesian inference (BI) was implemented under the CAT + GTR + G4 model, using PhyloBayes-MPI 1.8 [45,[48][49][50][51]. Two independent Markov Chain Monte Carlo (MCMC) runs of 8000 generations each were executed. A consensus tree was simultaneously built by pooling the remaining MCMC trees from both runs. Convergence was evaluated with the "bpcomp" and "tracecomp" procedures in the PhyloBayes package with a burn-in of the first 1000 generations. The maximum discrepancy in the convergence result is 0.017. The maximum likelihood (ML) inference was conducted in IQTREE v2.1.2 [52]. Substitution models were compared and selected according to the Bayesian Information Criterion (BIC) by using ModelFinder [53]. Additionally, a profile mixture model (C60) was used based on the best-fit substitution model of the NP datasets of amino acid datasets [54]. An edge-unlinked model was specified for both the full partition and the merged partition schemes. The best run is selected from the four independent runs based on log-likelihood. A total of 1000 Ultrafast bootstraps were used to evaluate the nodal support of the ML tree [55], and to estimate the consensus tree. For the "AA" matrix, three partition schemes were applied for ML ( Table 2): (1) no partition (NP); (2) full partition (FP) that provides the best-fitting model for each individual gene; and (3) merged partition (MP) that implements a greedy strategy starting with the full partition model and subsequently merging pairs of genes until the model fit does not improve any further. The phylogenetic-terrace aware (PTA) data structure was used to facilitate the efficient analysis of the "AA" matrix under each partition model [56]. We selected the best final maximum likelihood and consensus trees according to the Akaike Information Criterion (AIC). Phylogenetic analyses ranked nodes with posterior probabilities (PP) and bootstrap support values (BS) = 1/100 as fully supported, 0.98-0.99/95-99 as strongly 3. Results and Discussion 3.1. Morphology of Pseudoacanthocephalus bufonis ( Figure 1, Table 3) Trunks are medium-sized, smooth, and cylindrical. Females are much larger than males. Proboscis is nearly cylindrical, armed with 16-20 longitudinal rows of 3-5 rooted hooks each. Proboscis receptacle is double-walled with the cerebral ganglion at the posterior of the proboscis receptacle. The neck is short. Lemnisci are more or less equal, slightly longer or shorter than the proboscis receptacle. Morphometric data of the present specimens and morphometric comparisons of P. bufonis between our specimens and previous studies are shown in Table 3.   (Figure 1, Table 3)

Morphology of Pseudoacanthocephalus bufonis
Trunks are medium-sized, smooth, and cylindrical. Females are much larger than males. Proboscis is nearly cylindrical, armed with 16-20 longitudinal rows of 3-5 rooted hooks each. Proboscis receptacle is double-walled with the cerebral ganglion at the posterior of the proboscis receptacle. The neck is short. Lemnisci are more or less equal, slightly longer or shorter than the proboscis receptacle. Morphometric data of the present specimens and morphometric comparisons of P. bufonis between our specimens and previous studies are shown in Table 3.  Abbreviations: ac, acanthor; bu, copulatory bursa; cg, cement glands; ep, epidermis; h, hooks; le, lemniscs; me, membrane; p, proboscis; pr, proboscis receptacle; sh, shell of egg, te, testis; ub, uterine bell; va, vagina. The morphology and measurements of the present material are more or less identical to the previous descriptions of P. bufonis [31][32][33][34][35][36], including the morphology and size of trunk and proboscis, the number of the longitudinal rows of proboscis hooks and the hooks per longitudinal row, the number and length of testes and cement-glands, and the morphology and size of eggs. However, the lengths of the proboscis receptacle and lemnisci are slightly ser than those of the previous studies. Additionally, we also sequenced the ITS region (OQ550505, OQ550506) of our specimens. Pairwise comparison of the ITS sequences of our specimens with the available ITS data (KC491878-KC491883) of P. bufonis reported in the previous study [70], displayed only 0.17% nucleotide divergence. Thus, we confirmed our specimens to be P. bufonis.

Gene Content and Organization of the Mitogenome
The complete mitogenome of P. bufonis is 14,056 bp in length and includes 36 genes, containing 12 PCGs (cox1-3, nad1-6, nad4l, cytb, and atp6), 22 tRNA genes, 2 rRNA genes (rrnS and rrnL), and two non-coding regions (NCR1 and NCR2) ( Figure 2, Table 4). The lack of atp8 in the mitogenome of P. bufonis is typical for most of the available mitogenomes of acanthocephalans, except for Leptorhynchoides thecatus, which has two putative atp8 genes [55]. All genes in the mitogenome of P. bufonis are encoded on the same strand and in the same direction. Furthermore, the highest GC-skew (0.53) and the second lowest AT-skew (−0.28) of the mitogenome of P. bufonis in the order Echinorhynchida show its preference for G and T nucleotides (Figure 3), which was possibly a result of the propensity for low use of A-rich codons in their PCGs (Table 5). A similar situation also occurred in Polyacanthorhynchus caballeroi and some species of Polymorphida [19,21,25]. The morphology and measurements of the present material are more or less identical to the previous descriptions of P. bufonis [31][32][33][34][35][36], including the morphology and size of trunk and proboscis, the number of the longitudinal rows of proboscis hooks and the hooks per longitudinal row, the number and length of testes and cement-glands, and the morphology and size of eggs. However, the lengths of the proboscis receptacle and lemnisci are slightly ser than those of the previous studies. Additionally, we also sequenced the ITS region (OQ550505, OQ550506) of our specimens. Pairwise comparison of the ITS sequences of our specimens with the available ITS data (KC491878-KC491883) of P. bufonis reported in the previous study [70], displayed only 0.17% nucleotide divergence. Thus, we confirmed our specimens to be P. bufonis.

Gene Content and Organization of the Mitogenome
The complete mitogenome of P. bufonis is 14,056 bp in length and includes 36 genes, containing 12 PCGs (cox1-3, nad1-6, nad4l, cytb, and atp6), 22 tRNA genes, 2 rRNA genes (rrnS and rrnL), and two non-coding regions (NCR1 and NCR2) (Figure 2, Table 4). The lack of atp8 in the mitogenome of P. bufonis is typical for most of the available mitogenomes of acanthocephalans, except for Leptorhynchoides thecatus, which has two putative atp8 genes [55]. All genes in the mitogenome of P. bufonis are encoded on the same strand and in the same direction. Furthermore, the highest GC-skew (0.53) and the second lowest AT-skew (−0.28) of the mitogenome of P. bufonis in the order Echinorhynchida show its preference for G and T nucleotides (Figure 3), which was possibly a result of the propensity for low use of A-rich codons in their PCGs (Table 5). A similar situation also occurred in Polyacanthorhynchus caballeroi and some species of Polymorphida [19,21,25].
The composition and usage of codons in the mitogenome of P. bufonis were shown in Figure 4 and Table 6. ATN (i.e., ATA, ATG, and ATT), GTG, and TTG are used as start codons for the 12 PCGs in the mitogenome of P. bufonis, whereas TAA, TAG, and incomplete codons of T or TA are used as termination codons, in accordance with those of other acanthocephalans [15,[18][19][20][21][22][23][24][25][58][59][60][61][62][63][64]. GTG is the most common start codon, being used for six PCGs (cox1, cox3, nad2, nad4l, nad5, and nad6), followed by ATN for four PCGs (ATA: atp6; ATG: nad3 and nad4; ATT: nad1). Two genes (cytb and cox2) were inferred to use TTG as the start codon. Among the 12 PCGs, six genes (atp6, cytb, nad1, nad2, nad3, and nad4l) are terminated with complete stop codon TAA, while three genes (cox1, cox2, and nad5) were inferred to terminate with complete stop codon TAG. The incomplete stop codons T and TA are used for cox3 and nad4, and nad6, respectively.    In the PCGs of P. bufonis, the codon with the highest RSCU value is AGG (Ser), while the rarest codon is CTC (Leu). Val is the most frequently used amino acid (16.66%) in 12 PCGs of P. bufonis. Gln is the least commonly used amino acid (0.59%). The high frequency of Val (encoded by GTN) is associated with the high proportions of G and T in their protein-coding sequences (Figure 4 and Table 6).
In the mitogenome of P. bufonis, two rRNAs, rrnL, located between trnY and trnL1, and rrnS, located between trnM and trnP, were identified. The rrnL is 908 bp in length, with 62.86% A + T content, whereas the rrnS is 572 bp in length with 64.37% A + T content ( Figure 2 and Table 5). Animals 2023, 13, x FOR PEER REVIEW 14 of 21 In the mitogenome of P. bufonis, two rRNAs, rrnL, located between trnY and trnL1, and rrnS, located between trnM and trnP, were identified. The rrnL is 908 bp in length, with 62.86% A + T content, whereas the rrnS is 572 bp in length with 64.37% A + T content ( Figure 2 and Table 5).

Non-Coding Regions
In the mitogenome of P. bufonis, there are two non-coding regions (NCR1 and NCR2). NCR1 is located between trnI and trnM, is 610 bp in length. NCR2, located between trnW and trnV, is 503 bp. Their A + T contents are 61.97% and 56.86%, respectively (Table 5).

Molecular Phylogeny
Phylogenetic trees generated from BI and ML methods under different models have similar topologies and indicate that the Acanthocephala are monophyletic, which was widely accepted in previous studies. However, the evolutionary relationships of the The gene order of trnS2 in the mitogenome of P. bufonis is type A (trnS2, atp6, nad3, trnW, trnV, trnK, trnE, trnT). The gene arrangement of trnK is of type D (trnK, trnV). The trnS1 of P. bufonis is of type F (trnM, rrnS, trnF, cox2, trnC, cox3, trnA, trnR, trnN, trnS1) ( Figure 6).

Non-Coding Regions
In the mitogenome of P. bufonis, there are two non-coding regions (NCR1 and NCR2). NCR1 is located between trnI and trnM, is 610 bp in length. NCR2, located between trnW and trnV, is 503 bp. Their A + T contents are 61.97% and 56.86%, respectively (Table 5).

Molecular Phylogeny
Phylogenetic trees generated from BI and ML methods under different models have similar topologies and indicate that the Acanthocephala are monophyletic, which was widely accepted in previous studies. However, the evolutionary relationships of the Acanthocephala and the three subtaxa of Rotifera (Monogononta, Bdelloidea, and Seisonidea) have been under debate for a long time [19,[71][72][73]. The present phylogenetic results showed that the Acanthocephala is sister to Bdelloidea (Rotaria rotatoria, Philodina citrina) and rejected the monophyly of Eurotatoria (Monogononta + Bdelloidea), which are identical to the previous phylogenetic results using EST libraries [72] and mitogenomic data [19], but conflicted with some other phylogenies based on 18S rDNA and transcriptomic data [71,73].
Our phylogeny also supported the division of the phylum Acanthocephala into three large clades (Clade I, Clade II, and Clade III) (Figure 7).
Clade I, including Macracanthorhynchus hirudinaceus and Oncicola luehei (Oligacanthorhynchida: Oligacanthorhynchidae), represents Archiacanthocephala, a monophyletic group located at the base of the phylogenetic trees of Acanthocephala (Figure 7). The present results agree well with some previous phylogenetic studies [8,9,[11][12][13]18,20,[22][23][24][25]64,66]. The representative of Polyacanthocephala (Polyacanthorhynchus caballeroi) nested with species of Eoacanthocephala (Pallisentis celatus + Acanthogyrus bilaspurensis + Neoechinorhynchus violentum + Paratenuisentis ambiguus), forming Clade II. The present phylogenetic results challenged the validity of Polyacanthocephala, as have some previous molecular phylogenetic studies [19,22,23,25,64,66]. Acanthocephala and the three subtaxa of Rotifera (Monogononta, Bdelloidea, and Seisonidea) have been under debate for a long time [19,[71][72][73]. The present phylogenetic results showed that the Acanthocephala is sister to Bdelloidea (Rotaria rotatoria, Philodina citrina) and rejected the monophyly of Eurotatoria (Monogononta + Bdelloidea), which are identical to the previous phylogenetic results using EST libraries [72] and mitogenomic data [19], but conflicted with some other phylogenies based on 18S rDNA and transcriptomic data [71,73]. Our phylogeny also supported the division of the phylum Acanthocephala into three large clades (Clade I, Clade II, and Clade III) (Figure 7). Clade I, including Macracanthorhynchus hirudinaceus and Oncicola luehei (Oligacanthorhynchida: Oligacanthorhynchidae), represents Archiacanthocephala, a monophyletic group located at the base of the phylogenetic trees of Acanthocephala (Figure 7). The present results agree well with some previous phylogenetic studies [8,9,[11][12][13]18,20,[22][23][24][25]64,66]. The representative of Polyacanthocephala (Polyacanthorhynchus caballeroi) nested with species of Eoacanthocephala (Pallisentis celatus + Acanthogyrus bilaspurensis + Neoechinorhynchus violentum + Paratenuisentis ambiguus), forming Clade II. The present phylogenetic results challenged the validity of Polyacanthocephala, as have some previous molecular phylogenetic studies [19,22,23,25,64,66]. The representatives of Palaeacanthocephala formed Clade III. The monophyly of the order Polymorphida, including the representatives of Plagiorhynchus transversus, Polymorphus minutus, Southwellina hispida, Centrorhynchus clitorideus, C. milvus, C. aluconis, Sphaerirostris lanceoides, and S. picae, is strongly supported in our phylogenetic results. In Polymorphida, the Polymorphidae are more closely related to the Centrorhynchidae than the Plagiorhynchidae, in accordance with other recent mitogenomic phylogenies, but inconsistent with some previous phylogenetic studies based on nuclear and mitochondrial genetic markers [17,68,69,74]. Our phylogenetic results showed that the order Echinorhynchida is paraphyletic, which is consistent with previous molecular phylogenetic studies The representatives of Palaeacanthocephala formed Clade III. The monophyly of the order Polymorphida, including the representatives of Plagiorhynchus transversus, Polymorphus minutus, Southwellina hispida, Centrorhynchus clitorideus, C. milvus, C. aluconis, Sphaerirostris lanceoides, and S. picae, is strongly supported in our phylogenetic results. In Polymorphida, the Polymorphidae are more closely related to the Centrorhynchidae than the Plagiorhynchidae, in accordance with other recent mitogenomic phylogenies, but inconsistent with some previous phylogenetic studies based on nuclear and mitochondrial genetic markers [17,68,69,74]. Our phylogenetic results showed that the order Echinorhynchida is paraphyletic, which is consistent with previous molecular phylogenetic studies [12,14,28,30,69]. Furthermore, they supported the resurrection of Pseudoacanthocephalidae [16]. In the present mitogenomic phylogeny, P. bufonis clustered together with Cavisoma magnum, suggesting an affinity between Pseudoacanthocephalidae and Cavisomatidae. Our results agreed well with some recent phylogenetic studies based on nuclear gene sequences [16,75]. These indicated that the current classification of Echinorhynchida is based on unique combinations of characteristics, not shared derived features [13]. The systematics of Echinorhynchida needs to be revised so that its constituent families, subfamilies, and genera reflect the underlying lineages. This will require phylogenetic analysis of both nuclear and mitochondrial DNA datasets from representatives of a more diverse range of taxa than are currently available. It is essential to sequence mitogenomes from yet unrepresented taxa for constructing the molecular phylogenetic framework of Acanthocephala and further exploring the unusual patterns of mitogenomic evolution in this group. The complete mitogenome of P. bufonis obtained herein represents a valuable building block for future work.

Conclusions
In the present study, the complete mitochondrial genome of P. bufonis, the first representative of the family Pseudoacanthocephalidae, was characterized. Phylogenetic analyses based on the amino acid sequences of 12 protein-coding genes further confirmed the sister relationship of the Acanthocephala and Bdelloidea and rejected the monophyly of Eurotatoria (Monogononta + Bdelloidea) and Pararotatoria (Seisonidea + Acanthocephala). Our phylogeny also revealed that the order Echinorhynchida and the family Echinorhynchidae are both paraphyletic in the Acanthocephala. The current systematic status of Pseudoacanthocephalus in the Echinorhynchidae is challenged. The present phylogenetic results supported the recent resurrection of Pseudoacanthocephalidae and showed a close affinity between Pseudoacanthocephalidae and Cavisomatidae. Phylogenetic analyses also strongly supported the monophyly of the order Polymorphida and indicated that the Polymorphidae and Centrorhynchidae have a closer relationship than the Plagiorhynchidae. The present phylogenetic studies provided a new insight into the evolutionary relationships of higher taxa within Acanthocephala.