Phenotypic plasticity, genetic structure and systematic position of Neoechinorhynchus emyditoides Fisher, 1960 (Acanthocephala: Neoechinorhynchidae): a parasite of emydid turtles from the Nearctic and Neotropical regions

Abstract Abstract The taxonomy of the 10 recognized Neoechinorhynchus species associated with emydid turtles is complex due to the morphological conservatism. In the present study, specimens of N. emyditoides from northern and southeastern Mexico exhibit great phenotypic plasticity on its diagnostic characteristics. We sequenced three molecular markers: the internal transcribed spacers ITS1, ITS2 and 5.8S gene, the D2 + D3 domains of the large subunit from nuclear DNA and cytochrome c oxidase subunit I (cox1) from mitochondrial DNA. Sequences of the nuclear molecular markers were aligned and compared with other congeneric species associated with emydids available in GenBank. Phylogenetic analyses supported the polyphyly of Neoechinorhynchus. The species from emydids formed a clade, which was subdivided into five subclades that correspond with each species analysed (N. pseudemydis, N. chrysemydis, N. emydis, N. schmidti and N. emyditoides). To understand better the genetic structure of N. emyditoides a haplotype network was inferred with 29 cox1 sequences, revealing the presence of 13 haplotypes, two of which were shared and 11 were unique. The high values of fixation index, Fst (0.4227–0.8925) detected between the two populations from southeastern and the two from northern Mexico indicated low genetic flow among the populations. Our data suggest that the Neoechinorhynchus species associated with emydid turtles diversified in the eastern USA and that of N. emyditoides expanded its distribution range reached southeastern Mexico.


Introduction
Phenotypic plasticity is defined as the ability of a genotype to produce multiple phenotypes in response to environmental conditions (Roff, 2002;Miner et al., 2005). In organisms with complex and diverse life cycles such as parasites, different environmental conditions include (1) the host's immune system, (2) different host species and (3) the geographical distribution of the definitive hosts (Poulin, 2007). The recent application of molecular markers has helped to define, recognize, delineate and better understand intraspecific variation that can be attributed to differences in the development and phenotypic plasticity of acanthocephalans (Steinauer et al., 2007;Rosas-Valdez et al., 2012, 2020Alcántar-Escalera et al., 2013;Perrot-Minnot et al., 2018;Pinacho-Pinacho et al., 2018a;Lisitsyna et al., 2019).
The first morphological and molecular phylogenies inferred with a few species of Neoechinorhynchus revealed that the genus is paraphyletic (Pinacho-Pinacho et al., 2017). However, one of the most emblematic and enigmatic groups of the genus Neoechinorhynchus from North America associated with emydid turtles (Barger, 2004;Barger and Nickol, 2004) formed a monophyletic assemblage (García-Varela and Pinacho-Pinacho, 2018;Koch et al., 2021). The taxonomy of the 10 currently described species of Neoechinorhynchus associated with emydid turtles is complex due to their morphological homogeneity. Species delimitation relies entirely on characteristics of the female worm, including the contour of the posterior end and the shape and membrane structure of the fully formed egg, whereas male worms are nearly identical among species (Cable and Hopp, 1954;Barger, 2004;Barger and Nickol, 2004). In addition, Dezfuli and Tinti (1998) and Barger and Nickol (2004) mentioned that species of Neoechinorhynchus from turtles exhibit great phenotypic plasticity, which may lead to the misclassification of female specimens and leave males of different species completely indistinguishable. Barger (2004) performed one of the most comprehensive taxonomic reviews of Neoechinorhynchus species that infect emydid turtles from North America, acknowledging that eastern USA represents a biodiversity hot spot for this group of acanthocephalans. Barger reported that Neoechinorhynchus emyditoides Fisher, 1960 has the widest distribution range extending from the eastern USA to southeastern Mexico through the Gulf of Mexico.
As part of our long-term studies on the biodiversity of helminth parasites of turtles, acanthocephalans belonging to the Neoechinorhynchidae were recovered from the intestines of emydid turtles from five localities in northern and southeastern Mexico. Our extensive sampling allowed us recognize morphologically two species of acanthocephalans, i.e. Neoechinorhynchus schmidti Barger, Thatcher and Nickol, 2004 and N. emyditoides.
The objective of the current study was to combine morphological and molecular characteristics to investigate the phenotypic plasticity of N. emyditoides along its distribution range and explore the genetic structure of populations by using sequences of the cytochrome c oxidase subunit 1 (cox1) from mitochondrial DNA. In addition, we briefly discuss the systematics of Neoechinorhynchus species associated with emydid turtles by using sequences of two molecular markers: the D2 + D3 domain of the large subunit (LSU) and the internal transcribed spacer region, including 5.8S (ITS) from nuclear ribosomal DNA.

Morphological study
For taxonomic identification, specimens were stained with Mayer's paracarmine, dehydrated in a graded ethanol series, cleared with methyl salicylate and mounted on permanent slides with Canada balsam. Mounted specimens were examined under a bright-field Leica DM 1000 LED microscope (Leica, Wetzlar, Germany), and drawings were made using a drawing tube attached to the microscope. Measurements were taken using Leica Application Suite microscope software (Leica) and are given in micrometres (μm).
Vouchers of N. emyditoides from Tlacotalpan, Veracruz (No. 6695) and the new vouchers (Nos. 11647-11653) were examined and deposited in the Colección Nacional de Helmintos (CNHE), Instituto de Biología, Universidad Nacional Autónoma de México (UNAM), Ciudad de México, Mexico. The species identification was conducted following the revision of the genus Neoechinorhynchus by Barger (2004) and the original description (Fisher, 1960). For scanning electron microscopy (SEM), two specimens of N. emyditoides of each locality sampled were dehydrated with an ethanol series, critical point dried, sputter coated with gold and examined with a Hitachi Stereoscan Model S-2469N scanning electron microscope operating at 15 kV from the Instituto de Biología, UNAM.

Amplification and sequencing of DNA
A total of seven specimens identified as N. emyditoides were analysed. The posterior end of four specimens, two male and two female (hologenophores; Pleijel et al., 2008), were used for DNA extraction, whereas the rest of the body was stained with Mayer's paracarmine and mounted on permanent slides with Canada balsam. In addition, one female of N. emyditoides from Purification River was cut-off into three sections. The anterior and posterior sections were processed for SEM, whereas the middle section was processed for DNA extraction. Finally, other two specimens identified as N. emyditoides were placed individually in tubes and digested overnight at 56°C in a solution containing 10 mM Tris-HCl (pH 7.6), 20 mM NaCl, 100 mM Na 2 EDTA (pH 8.0), 1% sarkosyl and 0.1 mg mL −1 proteinase K. Following digestion, DNA was extracted from the supernatant using the DNAzol reagent (Molecular Research Center, Cincinnati, OH, USA) according to the manufacturer's instructions. The domains D2 + D3 and ITS from nuclear ribosomal DNA and the cox1 from mitochondrial DNA were amplified using polymerase chain reaction (PCR). The domains D2 + D3 from LSU of approximately 700 bp were amplified using forward primer 502, 5 ′ -CAA GTA CCG TGA GGG AAA GTT GC-3 ′ and reverse primer 536, 5 ′ -CAG CTA TCC TGA GGG AAAC-3 ′ (García- Varela and Nadler, 2005). A fragment of approximately 750 bp from ITS was amplified using the forward primer BD1, 5 ′ -GTC GTA ACA AGG TTT CCG TA-3 ′ and the reverse primer BD2, 5 ′ -ATC TAG ACC GGA CTA GGC TGT G-3 ′ (Luton et al., 1992). The cox1 of approximately 550 bp was amplified using the forward primer 516, 5 ′ -ATT TTT TAG TTT GAG TGT GAG GAG-3 ′ and the reverse primer 517, 5 ′ -ATGA CGA ATT AAT ATT ACG ATC CA-3 ′ (see Pinacho-Pinacho et al., 2018a, 2018b. The amplification reactions (25 μL) consisted of 1 μL of each primer (10 μM), 2.5 μL of 10× buffer, 1.5 μL of 2 mM MgCl 2 , 0.5 μL of dNTPs (10 mM), 16.37 μL of water, 2 μL of genomic DNA and 1 U of Taq DNA polymerase (Platinum Taq, Invitrogen Corporation, São Paulo, Brazil). The PCR cycling conditions for amplification included denaturation at 94°C for 3 min, followed by 35 cycles of 94°C for 1 min, annealing at 50°C for LSU, ITS and 40°C for cox1 for 1 min and extension at 72°C for 1 min, with a final postamplification incubation at 72°C for 10 min. The sequencing reactions were performed using the initial primers for LSU, ITS and cox1, plus two internal primers: 503, 5 ′ -CCT TGG TCC GTG TTT CAA GAC G-3 ′ and 504, 5 ′ -CGT CTT GAA ACA CGG ACT AAGG-3 ′ for LSU (García-Varela and Nadler, 2005). Sequencing reactions were performed using ABI Big Dye (Applied Biosystems, Boston, MA, USA) terminator sequencing chemistry, and the reaction products were separated and detected using an ABI 3730 capillary DNA sequencer. Contigs were assembled, and base-calling differences were resolved using CodonCode Aligner 9.0.1 (CodonCode Corporation, Dedham, MA, USA). Newly generated sequences were deposited in the GenBank database under the accession numbers; Hologenophore, male (OM892136 and OM892138 for LSU; OM892142 and OM892144 for ITS; OM891889 for cox1), female (OM892137 and OM892139 for LSU; OM892143 and OM892145 for ITS; OM891888 and OM891890 for cox1). The accession numbers for the female of N. emyditoides cut-off into three sections were OM892146 for ITS and OM892140 for LSU. Finally, the accession numbers of the other specimens were OM892141 for LSU and OM892147-148 for ITS.

Alignments and phylogenetic analyses
Sequences obtained in the current research for LSU and ITS were aligned with other congeneric species downloaded from GenBank (see Table 1). In addition, sequences of the LSU and ITS of other genera, such as Polyacanthorhynchus Travassos 1920, Acanthosentis Verma and Datta 1929, Floridosentis Ward 1953, Mayarhynchus Pinacho-Pinacho, Hernández-Orts, Sereno-Uribe, Pérez-Ponce de León and García-Varela, 2017, and Atactorhynchus Chandler 1935, were used as outgroup. Sequences of each nuclear molecular marker were aligned separately using the software Clustal W with default parameters implemented in MEGA version 7.0 (Kumar et al., 2016). The best-fitting nucleotide substitution models for LSU and ITS dataset were GTR + G + I and were estimated with the Akaike Information Criterion (AIC) implemented in MEGA version 7.0 (Kumar et al., 2016). The phylogenetic analyses were inferred through maximum likelihood (ML) with the program RAxML v7.0.4 (Silvestro and Michalak, 2012) and Bayesian inference (BI) criteria, employing the nucleotide substitution model identified for AIC. BI trees were generated using MrBayes v3.2 (Ronquist et al., 2012), running two independent MC3 runs of four chains for 5 million generations and sampling tree topologies every 1000 generations. 'Burn-in' periods were set to 1 million of generations according to the standard deviation of split frequencies values (<0.01). To support each node, 10 000 bootstrap replicates were run with the ML method. Posterior probabilities of clades were obtained from 50% majority rule consensus of sample trees after excluding the initial 20% as 'burn-in'. The genetic divergence among taxa was estimated using uncorrected 'p' distances with the program MEGA version 7.0 (Kumar et al., 2016).
To examine cox1 haplotype frequency among the populations of N. emyditoides, an unrooted statistical network was constructed using the PopART program with the median-joining algorithm (Bandelt et al., 1999). The degree of genetic differentiation among the populations was estimated using the fixation indices F st (Hudson et al., 1992), with the program Arlequin v.3.5 (Excoffier and Lischer, 2010). To investigate the population history and demography, Tajima's D (Tajima, 1989) and Fu's F (Fu, 1997) tests were calculated using DnaSP v. 5. 10 (Rozas et al., 2003). Values less than 0.05 were considered statistically significant.

Morphological identification
The acanthocephalans were collected from several previous field expeditions in four localities ( Fig. 1; Table 1), two from the northern Mexico and two from the southeastern Mexico (Figs 2 and 3). The acanthocephalans exhibited the typical morphological characteristics of N. emyditoides. For example, trunk elongated cylindrical, slender, curved ventrally, body wall with reticular lacunar system, usually five (occasionally three) dorsal giant hypodermic nuclei and one ventral ( Fig. 2C-D). Proboscis nearly cylindrical, longer than wide (Figs 2A and 3A-D). Anterior hooks markedly larger than middle hooks not in a perfect circle but at different levels, posterior hooks smaller than middle hooks (Figs 2A and 3A-H). Neck short, wider than longer. Lemnisci appreciably different in length; longer than the proboscis receptacle. Proboscis receptacle single-walled (Fig. 2B). Male with two oblong testes, contiguous or slightly overlapping, equatorial, just behind anterior trunk; anterior testis relatively longer than the posterior testis. Prominent cement gland   ) elongated, longer than anterior testis. Rounded-ovoid cement gland reservoir with two lateral ducts. Saefftingen's pouch beginning at level of anterior end of the cement reservoir duct (Fig. 2E).
Posterior end of the female usually rounded with two lobes (Figs 2F and 3I-L), short vagina, slim uterus and uterine bell attached to the anterior body wall (Fig. 2F). Elliptical eggs, slightly inflated at the poles (Fig. 2G). These characteristics of our new morphometric data correspond to those reported previously (see Table 2) (Bravo-Hollis, 1946;Fisher, 1960;Barger, 2004).

Remarks
Fisher (1960) -Hollis, 1946;Fisher, 1960). Barger (2004) and Barger and Nickol (2004) performed one of the most comprehensive taxonomic reviews of the species of Neoechinorhynchus that infects emydid turtles, mentioning that N. emyditoides is the species with the most extensive distribution range including 11 states from eastern USA and single state (Veracruz) from southeastern Mexico. In the current research, we analysed the morphology of specimens from four populations of N. emyditoides from Mexico, exhibiting a wide phenotypic plasticity along its distribution. For instance, newly collected material showed some level of morphological intraspecific variation with respect to previous descriptions performed by Fisher (1960), Bravo-Hollis (1946) and García-Varela et al. (2011) (see Table 2). Interestingly, the specimens from Mexico possess lower limits for the following characteristics: body size, receptacle proboscis and proboscis hooks length with respect to original description (see Table 2). Cable and Hopp (1954) mentioned that the size and shape of fully developed eggs and the posterior end of the female are morphological traits keys on the differentiation at species level. For instance, the female of N. emyditoides is characterized by possessing a posterior end rounded with two lobes. Our specimens showed phenotypic plasticity on the shape of the posterior end (see Fig. 3I-L) (see Barger, 2004).

Phylogenetic analyses
The LSU dataset was conformed with 830 characters and 82 sequences, including 23 sequences of 10 Neoechinorhynchus species parasitizing teleost fishes and 46 sequences of five Neoechinorhynchus species associated to emydid turtles plus other sequences from diverse genera that were used as outgroups (see Table 1). The phylogenetic analyses inferred with ML and BI yielded that the genus Neoechinorhynchus is polyphyletic, due to the fact that several genera that were used as outgroups such as Mayarhynchus and Atactorhynchus from the Neoechinorhynchidae were nested inside of Neoechinorhynchus (Fig. 4A). The five species of Neoechinorhynchus analysed associated with emydid turtles formed a clade, which was subdivided into five subclades that correspond with each species analysed as N. pseudemydis, Neoechinorhynchus chrysemydis Cable and Hopp, 1954, N. emydis, N. schmidti and N. emyditoides with high bootstrap values and Bayesian posterior probabilities (Fig. 4A). In addition, all the samples (including six newly generated sequences) identified as N. emyditoides recovered in the four localities from northern and southeastern Mexico formed a subclade together with other sequence identified previously as N. emyditoides available in the GenBank (MK238158) from Oklahoma, USA (Fig. 4A). The  Song et al. (2014) Sequences in bold were generated in the current study. a Source for sequences misidentified as Octospiniferoides sp.
intraspecific genetic divergence among the isolates of N. emyditoides was low that ranged from 0 to 1.6%. Two juvenile worms previously identified as N. emyditoides from Monterrey, Mexico (KR086331-332) were nested in a subclade together with other two sequences available in GenBank (MK238159-160) identified as N. pseudemydis, a parasite of the red-eared slider turtle from Oklahoma, USA. This subclade was supported with high bootstrap values and Bayesian posterior probabilities (Fig. 4A). A second dataset with the ITS was conformed with 886 characters and 120 sequences. This alignment included 27 sequences of 12 Neoechinorhynchus species parasitizing teleost fishes and 84 sequences of five Neoechinorhynchus species associated with emydid turtles plus other sequences from diverse genera that were used as outgroups (see Table 1). The phylogenetic analyses inferred with ITS dataset agreed with the LSU tree with respect to polyphyly of Neoechinorhynchus due to the fact that the several genera that were used as outgroup taxa namely Mayarhynchus, Atactorhynchus and Floridosentis from the Neoechinorhynchidae were nested inside Neoechinorhynchus (Fig. 4B). The five species of Neoechinorhynchus associated with   emydid turtles formed a single clade, which was subdivided into five subclades that correspond with each species analysed as N. pseudemydis, N. chrysemydis, N. emydis, N. schmidti and N. emyditoides with high bootstrap values and Bayesian posterior probabilities (Fig. 4B). All samples (including seven newly generated sequences) identified as N. emyditoides recovered in the four localities from northern and southeastern Mexico formed a subclade together, with other three sequences identified previously as N. emyditoides available in the GenBank dataset (MW520495-497) from Oklahoma and Texas, USA (Fig. 4B). The intraspecific genetic divergence among the isolates of N. emyditoides was low that ranged from 0 to 0.5%. Two juvenile worms previously identified as N. emyditoides from Monterrey, Mexico (MG870860-861) were nested in a subclade together with other seven sequences available in GenBank (MK238090-901, MW520486-491) identified as N. pseudemydis, plus other nine sequences of Neoechinorhynchus sp. (MW520472-480), recovered from the ostracods of the genus Physocypra from Oklahoma, USA (Fig. 4B).

Haplotype network
A haplotype network was built with the cox1 alignment of 546 bp and contained sequences of 29 specimens (including three newly generated sequences) identified as N. emyditoides from four populations two from northern Mexico (localities; Monterrey, Nuevo León and Purificación river, Tamaulipas) and other two from southeastern Mexico (localities; Catemaco, Veracruz and Tlacotalpan, Veracruz). The haplotype network revealed the presence of 13 haplotypes, two of them (H1, H3) were shared between the localities of Catemaco, Veracruz and Monterrey, Nuevo León and 11 were unique haplotypes (Fig. 5). The level of haplotype diversity (Hd = 0.913) was very high, and nucleotide diversity was low (pi = 0.03728) among the populations (  Table 4).

Discussion
To the best of our knowledge, two species of the genus Neoechinorhynchus (N. schmidti and N. emyditoides) associated with emydid turtles have been recorded in Mexico (Bravo-Hollis, 1946;García-Varela et al., 2011;Pinacho-Pinacho et al., 2018a, 2018b. The current records suggest that the species do not have a sympatric distribution; for example, N. emyditoides has been recorded from northern Mexico to the Papaloapan River basin, which is considered the second largest hydrobiological system in southeastern Mexico, whereas N. schmidti has been recorded in the Centla Swamps, which is one of the largest wetlands in southeastern Mexico and is formed by the delta of the Grijalva and Usumacinta rivers (González-Ramírez and Parés-Sierra, 2019). The species N. emyditoides was described in red-eared slider turtles (T. scripta elegans) from the St. Francis River, Arkansas, USA (Fisher, 1960). This acanthocephalan is considered one of the biogeographical core parasite fauna of emydid turtles in the eastern USA (Barger, 2004;Barger and Nickol, 2004). In the current study, adults recovered from several emydid turtles from northern and southeastern Mexico were identified as N. emyditoides. Our observations and morphometric data found remarkable morphological differences among the populations from Mexico (see Table 2), including the contour posterior end of the females (Fig. 3). The inclusion of molecular characteristics in this study was key to achieving a better understanding of the Longer lemiscus L (mean)
The haplotype network analysis of cox1 sequences inferred with 29 sequences revealed the presence of 13 haplotypes, two of which (H1, H3) were shared between the localities of Catemaco, Veracruz and Monterrey, Nuevo León, and 11 were unique haplotypes (Fig. 5). The two populations from southeastern Mexico (Catemaco and Tlacotalpan, Veracruz) had low F st value (0.1603), suggesting genetic flow between both populations, which can be explained by the fact that both localities belong to the same hydrological basin (González-Ramírez and Parés-Sierra, 2019). In contrast, the two localities from northern Mexico (Monterrey, Nuevo León and Purificación river, Tamaulipas) did not belong to the same freshwater system. However, the value of F st was also low (0.1848), which suggests that definitive or intermediate hosts have the capacity to disperse and maintain the connection between the two localities.
Our phylogenetic analyses inferred with the LSU and ITS datasets confirmed that Neoechinorhynchus is polyphyletic, with most of its members nested in several independent clades that did not share a common ancestor. However, a main clade was formed with the five species of Neoechinorhynchus associated with emydid turtles. This clade contained three species distributed in the eastern USA (N. emydis, N. pseudemydis and N. chrysemydis), one in southeastern Mexico (N. schmidti) and one species (N. emyditoides) distributed from the eastern USA to southeastern Mexico ( Fig. 4A and B). All the samples identified as N. emyditoides recovered from the four locations from northern and southeastern Mexico formed a subclade together with other sequences previously identified as N. emyditoides available in GenBank for LSU (MK238158) and ITS (MW520495-497) from Texas and Oklahoma, USA ( Fig. 4A and B). Two juvenile worms from red-eared sliders in Monterrey, Mexico (KR086331-332 and MG870860-861 for LSU and ITS, respectively) previously identified as N. emyditoides were nested in a subclade together with other sequences identified as N. pseudemydis. Koch et al. (2021) performed one of the most extensive studies of morphological, ecological and molecular data of the snail, ostracod and turtle hosts of Neoechinorhynchus species from the eastern USA. Their analysis found that the sequences from northern Mexico originally identified as N. emyditoides (KR086331-332 and MG870860-861) corresponded to N. pseudemydis. The intraspecific genetic divergence estimated on the current study among the isolates of N. emyditoides ranged from 0.3 to 1.6% and from 0 to 0.5% for LSU and ITS, respectively. These values of intraspecific genetic divergence are similar to those previously reported for isolates of N. chrysemydis, N. schmidti and N. emydis, which ranged from 0.3 to 1.6% for ITS and from 0 to 0.01% for LSU (García-Varela et al., 2011;Koch et al., 2021).
Members of the Emydidae family have been divided into two monophyletic lineages recognized as the subfamilies or complexes, Deirochelyinae and Emydinae (Gaffney and Meylan, 1988;Spinks et al., 2016). Deirochelyinae includes six genera broadly distributed in North America but also contains a few taxa that extend across the Greater Antilles, Mexico, Central America and South America (Parham et al., 2013(Parham et al., , 2015. Four of the six genera from Deirochelyinae (Trachemys, Graptemys, Pseudemys and Chrysemys), which share a common ancestor (Thomson et al., 2021), have been found to harbour Neoechinorhynchus species ( Barger and Nickol, 2004). The fossil record suggests that Deirochelyine originated in North America approximately 20.91 Ma prior to the Miocene (Spinks et al., 2016). Barger (2004) suggested that the eastern USA represents a biodiversity hot spot for Neoechinorhynchus species (and the same for emydids; see Lindeman, 2013) that infect emydid turtles. The presence of N. emyditoides in northern and southern Mexico could have resulted from ancestral populations of emydid turtles that inhabited the eastern USA and colonized new habitats along the coasts to the south due to the relatively warm and wet temperatures during the mid-Miocene (Spinks et al., 2016), which is mirrored in the diversity of both the hosts and their parasites.
Supplementary material. The supplementary material for this article can be found at https://doi.org/10.1017/S003118202200049X.