Internal transcribed spacer (ITS) barcoding reveals species differentiation and evolutionary lineages of Trypanorhyncha (Cestoda) in teleosts imported to Jordan from the Mediterranean and Red Seas

General Data and Parasitological Analysis

The Jordanian Slaughterhouses Department of the Amman Municipality provided specimens of parasitized fish (n=36) between September 2023 and July 2024. The examined samples were obtained from the Egyptian Mediterranean Sea Fisheries' Thunnus or tuna fish, Xiphias gladius (Abu Saif fish), and Solea aegyptiaca (Mousa fish). The Red Sea species under study include Plectropomus areolatus (Najel or the red grouper) and Epinephelus tauvina (commonly known as Hammor or Arabian grouper). Whole fish were eviscerated and subjected to macroscopical examination to assess the presence of larval cestodes and associated focal tissue necrosis within the coelomic cavity and somatic musculature (Fig. 1). Following initial verification of infection, specimens were transported on ice to the Laboratory of Parasitology at the Department of Applied Biology, Al-Balqa Applied University. The fish were eviscerated to be macroscopically inspected for the presence of tapeworms and associated necrosis in the body cavity and muscles (Fig 1). Inspections revealed that the fish had a tapeworm infection. Nevertheless, the precise quantity of larvae has not been determined. The frozen fish were subsequently sent to the laboratory of parasitology at Al-Balqa Applied University's department of applied biology. The visible encapsulated larvae (plerocerci and merocercoids), were removed to collect the tapeworm specimens, which were found in the muscle and viscera, as previously described by Palm and coworkers (Palm et al., 2009). Following their placement in Petri dishes containing a 0.65 % NaCl solution, the larval cestodes were examined under an optical microscope. The infection status was recorded as encapsulated larvae embedded in fish flesh. Larval cestodes were frozen at −20°C for the molecular identification.

Fig. 1.

The Isolated plerocercoids were thoroughly rinsed in sterile phosphate-buffered saline (PBS), in order to remove host contaminants. Individual larvae were mechanically homogenized using sterile pestles within 1.5 mL microcentrifuge tubes; containing a tissue lysis buffer. Homogenization was enhanced by the addition of Proteinase K, followed by an extended digestion phase; exceeding two hours at (56°C) in order to guarantee total tissue breakdown (Al-Khlifeh et al., 2021, 2023a). The total genomic DNA was extracted using the G-spin Total DNA Extraction Kit (iNTRON Biotechnology, Seongnam, Korea), according to modified manufacturer protocols. To optimize yields from highly fibrous helminth tissue, the recommended volumes of lysis solution, absolute ethanol, and Proteinase K were doubled. Extracted DNA was incubated on the silica membrane for 15 minutes before final elution. DNA yield, purity, and (260/280nm) absorbance ratios were quantified using a NanoDrop spectrophotometer.

PCR Amplification and Sequencing

The complete internal transcribed spacer (ITS) region was amplified using the universal primer pair BD1/DB1 (forward: 5′-GTCG-TAACAAGGTTTCCGTA-3′) and S4/4S (reverse: 5′-TCTAGATG-CGTTCGAA(G/A)TGTCGATG-3′), which target the conserved 3′ terminus of the 18S rRNA gene and the 5′ terminus of the 28S rRNA gene, respectively (Luton et al., 1992). Since these primers target the deeply conserved structural flanking regions, the resulting amplicons encompass variable lengths of flanking 18S and 28S rRNA genes on either side of the core ITS1–5.8S–ITS2 cassette, a layout standard for eucestode barcoding (Scholz & Kuchta, 2022).

Amplification protocol includes four minutes initial denaturation at 94°C, 40 cycles. To separate the double-stranded DNA, the template DNA was denatured for 30 seconds at 94°C. The ITS1 gene was annealed for 25 seconds at 49°C, whereas the 28S rRNA gene was annealed for 25 seconds at 54°C. The extension stage lasted one minute and twenty seconds at 74°C. To guarantee full amplification of all the DNA fragments, a final extension was carried out at 74°C for 4 minutes after the cycle steps. The target genes were effectively amplified for subsequent investigation via our modified PCR technique. A nucleic acid dye and 1.5 % (w/v) agarose gel electrophoresis were used to separate the PCR products. Digital pictures of the PCR product bands under UV light were taken via gel documentation equipment. Zymo Research DNA Clean & Concentrator® reagents were used to purify the PCR amplicons and remove excess primers and nucleotides. PCR products of appropriate quality and purity were sequenced via the Sanger method (ABI sequencer at Microgen, Korea), and the sequences were analyzed to confirm the species and construct a phylogenetic tree.

Quality Control and Sequence Filtering

All sequences underwent rigorous quality control, which included verification of a minimum length (300 bp), ambiguity filtering (sequences with >5 % Ns were excluded), BLAST validation against the NCBI nr database (E value < 1×10⁻⁵), taxonomic verification via the NCBI Taxonomy Database, redundancy removal, and manual screening for chimeras. Chimera screening was carried out by inspecting BLAST hit-coverage profiles for each sequence; sequences with non-overlapping high-identity hits to multiple distinct lineages would be flagged for exclusion. No sequences from this study were excluded on these grounds. The PV189195 and PV189196 sequences were excluded because of quality deficiencies, and EU343734 (Gyrocotylidea outgroup) was removed following initial rooting.

Initial Sequence Identity Assessment

Our preliminary sequence identity assessment via BLASTN with the megablast algorithm indicated that the putative new species isolates shared unusually low identity values with their nearest GenBank matches: PV917196 (81.5 %), PV917197 (81.5 %), and PV917198 (81.4 %). These identities are well below the conventional species-level threshold of ~97 % for ITS sequences in helminths, providing initial evidence of distinctiveness and motivating a rigorous phylogenetic investigation.

Phylogenetic and Compositional Analyses

Two complementary analytical frameworks were applied to a combined dataset of 131 ITS sequences (comprising 22 newly generated sequences and 109 GenBank references): (i) an alignment-based Neighbor-Joining strategy calculated from Tamura-Nei (TN93) (Tamura & Nei, 1993) distances, and (ii) an alignment-free k-mer composition strategy. This dual-track approach accounted for the pronounced length heterogeneity inherent to eucestode ITS amplicons (ranging from 318 to 10,312 bp in this dataset), which stems from variable tracking of flanking structural rRNA genes across diverse historical submissions.

Multiple Sequence Alignment and Column Trimming

A total of 131 sequences were aligned using MAFFT v7.490 (Katoh & Standley, 2013) with the ‘--auto’ option, which selects the FFT-NS-2 strategy for datasets of this size. Because of the length heterogeneity described above, the initial alignment spanned 12,999 positions, the majority of which were long gap-only regions contributed by sequences that contained only partial flanking rRNA segments. A custom Python script removed alignment columns in which more than 50 % of sequences carried a gap, thereby retaining only the conserved core ITS region in which all sequences carry informative bases. The trimmed alignment comprised 728 positions (94.4 % of columns removed), corresponding closely to the expected 700 – 800 bp ITS-1–5.8S length for trypanorhynch cestodes (Palm et al. 2009). The 50 % threshold was selected after testing 30 %, 50 % and 70 %: 30 % was too permissive and retained informative-site-poor regions, while 70 % was too aggressive and eroded the ITS-2 signal. Both the raw MAFFT alignment (Supplementary File S1) and the trimmed alignment used for downstream analyses (Supplementary File S2) are deposited with this article.

Clarification regarding the alignment length: The initial alignment length of 12,999 positions is a direct consequence of length-het-erogeneous input sequences, not an artefact of poor alignment quality: MAFFT must accommodate the longest input sequence by inserting gaps in shorter ones. The reduction to 728 bp after gap-trimming is therefore the standard procedure for recovering the homologous core for phylogenetic analysis and should not be interpreted as loss of information.

Evolutionary distance estimation

The pairwise genetic distances were calculated under the Tamura-Nei (1993) model (TN93), which indicates that base frequencies are unequal (π_A ≈ 0.199, π_C ≈ 0.220, π_G ≈ 0.313, π_T ≈ 0.268; estimated from the trimmed alignment), and separate transition rates for A↔G and C↔T pairs. The TN93 was chosen in preference to simpler models (JC69, K2P) because our sequences exhibit a clear GC bias and transition/transversion rate asymmetry, both of which these models assume to be absent. All the distances were computed in python using a custom implementation following the closed-form analytical solution (Tamura & Nei, 1993).

Tree construction and support

A Neighbour-Joining (NJ) tree (Saitou & Nei, 1987) was inferred from the 128 × 128 TN93 distance matrix; three duplicated reference accessions were collapsed prior to inference. Branch support was estimated using bootstrap resampling of the trimmed alignment, with 100 pseudoreplicates (Felsenstein, 1985). For each pseudoreplicate, columns were resampled with replacement, TN93 distances were recalculated, and a new NJ tree was constructed. Bootstrap support ≥70 % was considered strong, 50 – 69 % moderate, and <50 % weak (Hillis & Bull, 1993). To enable consistency checking, we also constructed a UPGMA dendrogram (average linkage) from the same TN93 matrix, which recovered the same three major groups without assuming a molecular clock at individual branches.

Two study isolates (PV920021 and PV920022, C. gracilis from Epinephelus tauvina), were included in the dataset, but they were excluded from the Neighbour-Joining tree, because their sequenced region (≈340 bp) overlapped only partially with the trimmed alignment, leading to TN93 distance saturation (d = 2.0) against many references. These sequences are retained in Supplementary Files S5 and S6, and their species identity (>97.99 % BLAST identity to C. gracilis references) is shown in Table 1. Therefore, he final phylogeny tree presents 121 sequences.

Alignment-free k-mer composition analysis

To provide an alignment-independent confirmation of phylogenetic signal, the frequency of all 256 possible tetranucleotides (4-mers) was computed independently for each sequence by sliding a window of length (k = 4) along the full-length sequence and recording only those windows composed exclusively of A/T/G/C (Al-Khlifeh et al., 2026). Counts were normalised to frequencies, yielding a 128 × 256 matrix. The value k = 4 was selected after evaluating k ⋲ {3, 4, 5, 6} on the basis of: (i) balance between descriptive power and overfitting risk at 128 sequences; (ii) biological relevance of tetranucleotides to rRNA secondary-structure motifs; and (iii) perfect discriminability of the putative new species (F1 = 1.0 at k = 4). Because this analysis does not require sequence alignment, it is robust to the length heterogeneity of the input dataset.

Compositional dissimilarity and hierarchical clustering

Pairwise compositional dissimilarity between k-mer vectors was quantified with the Bray–Curtis dissimilarity index, which is appropriate for normalized frequency data and is bounded in [0, 1] (Al-Khlifeh et al., 2026). The Bray–Curtis matrix was subjected to hierarchical clustering with Ward's linkage (Bakker, 2026), producing a dendrogram complementary to the distance-based phylogeny. Principal Component Analysis (PCA) of the standardised 256-dimensional k-mer matrix was used to visualise compositional relationships.

Classification and feature importance

Binary classification (putative new species vs. all others) was performed with Random Forest (100 trees; Breiman 2001), XGBoost (100 boosting rounds, learning rate 0.1, max depth 6; (Chen & Guestrin, 2016)) and L2-regularised logistic regression. All classifiers were evaluated with 10-fold stratified cross-validation (random_ state = 42). Given the small sample size of the positive class (n = 3), classification scores were treated as descriptive rather than generalisable. Feature importance was quantified with SHAP values (TreeExplainer, (Lundberg & Lee, 2017)) on the XGBoost model, and the top 20 discriminative k-mers were presented in the results..

Software and reproducibility

All analyses were performed in Python 3.9 using NumPy 1.23, Pandas 1.5, SciPy 1.10, Biopython 1.81, scikit-learn 1.2, XGBoost 1.7, SHAP 0.41, matplotlib 3.7 and seaborn 0.12. Alignment was performed with MAFFT v7.490. All stochastic analyses used ‘random_state = 42’ for reproducibility. Analysis scripts and intermediate files are deposited as Supplementary Material.

Morphological and Parasitological Observations

Tapeworm larvae were collected from the 36 examined fish. These larvae were first morphologically referred to the order Trypanorhyncha. Merocercoid larvae excised from the somatic musculature were morphologically assigned to the genus Molicola Dollfus, 1935. Conversely, plerocercoid larvae were localized predominantly within the coelomic cavity and visceral mesenteries, often presenting as encapsulated forms closely to the host musculature (Fig. 1).

Molecular Profile and Sequence Analysis

The PCR amplification of the complete ITS rRNA cassette, yielded an amplicon footprint of approximately 750 bp across the viable isolates. Sequence identity screening and annotation against the NCBI GenBank database indicated that the majority of the new sequences obtained showed clear structural homology to C. gracilis (Rudolphi, 1819) Pintner, 1931. These isolates fulfilled high-stringency sequence similarity criteria, yielding BLASTN identity values exceeding 97.5 %, and high query coverage metrics (see Table 1). Notably, initial macroscopical and superficial morphological screening tentatively referred the somatic merocercoid larvae to the genus Molicola Dollfus, 1935 based on structural larval boundaries. However, the ITS sequence analysis definitively unmasked these specimens as C. gracilis. This distinct mismatch highlights the known limitations of larval-stage morphometrics in trypanorhynch cestodes, where cryptic inter-generic features and host-induced encapsulation tissue artifacts frequently confound classical taxonomic identification, establishing genetic barcoding as an indispensable tool for definitive helminth diagnosis. A singular exception was noted for isolate PV920020, which showed a sequence identity threshold of 96.3 % which is very significant although slightly lower (see Table 1). The match is likely genuine, and not random because the E value is extraordinarily low.

In contrast, three distinct isolates (PV917196, PV917197, and PV917198) recovered exclusively from S. aegyptiaca demonstrated highly divergent molecular profiles. These isolates yielded top BLASTN matches of similarity of 81.5 % across narrow query coverage windows (<70 %), falling drastically below a 97 % species level similarity threshold for helminth ITS sequences (see Table 1). Thus, these three divergent isolates could not be assigned to any known valid species and they are now classified as a putative novel lineage within the Eucestoda.

Dataset Architecture and Comparative Taxonomy

Trypanorhyncha is known for its cryptic biodiversity. We compared local sequencing samples with global data from the GenBank. a comprehensive master matrix comprising 128 distinct sequences was assembled. This dataset architecture integrates 20 newly generated local isolates, 103 global reference sequences spanning major valid Eucestoda orders (including Bothriocephalidea, Diphyllobothriidea, Lecanicephalidea, and Proteocephalidea), and 5 highly divergent outgroup sequences (Gyrocotyle urna [n = 4]) and Brachycladium atlanticum [n= 1]) utilized to anchor basal compositional relationships. Our approach provides significant advantages in understanding the diversity of the isolates under investigation in addition to balancing the phylogenetic and geographic data (Fig. 2). To optimize taxonomic resolution—that is, the capacity of the ITS marker to differentiate between closely related taxa—reference selection was restricted to sequences with high similarity (>94 %) to sequences generated in this study.

Fig. 2.

The comparative matrix displayed in (Fig. 2) highlights a well-represented dataset for Trypanorhyncha in Mediterranean marine teleosts. Direct sequence identity evaluations between global reference entries and our study isolates confirmed that 17 of the recovered sequences—spanning both Red Sea and Mediterranean hosts—exhibited high levels of structural homology, consistently exceeding the standard species-delineation threshold (>97.5 %). Nevertheless, the sequence of the IDs, PV917196, PV917197, and PV917198, is below the normal threshold for species-level identity in comparison with other sequences of Trypanorhyncha. They most likely belong to a novel species or a distinct variant (Fig. 3). Following the exclusion of the 5 outgroup sequences and the removal of 2 short local isolates (PV920021 and PV920022) due to sequence length limitations, the final alignment-based phylogenetic framework presented a core matrix of 121 sequences (Fig. 4).

Fig. 3.

Fig. 4.

Model-Based Phylogenetic Topology

Both the NJ and MP methods converged on the same topology (log-likelihood: −13,359.68; parsimony score: 34,805 steps). Their agreement suggests a robust evolutionary relationship (Fig. 4). Isolates under study were divided into three significant clades. All six Eucestoda orders formed distinct, well-supported clades, confirming strong phylogenetic signals at quite high taxonomic levels. Isolates under study were divided into three significant clades: First, the three specific novel isolates (PV917196, PV917197, and PV917198) clustered as a monophyletic group positioned basally near the order Lecanicephalidea (Fig. 4). This novel cluster was characterized by short internal branch lengths (0.014) substitutions per site), a divergence footprint approximately 7- to 10-fold greater than the typical baseline intra-specific variation observed within valid Trypanorhyncha (<0.002) substitutions per site), thus approaching recognized inter-specific evolutionary distances.

The second clade comprised 12 study isolates displaying high intra-clade consistency (mean distance <0.003 substitutions per site) that clustered in close proximity to reference Lecanicephal-idea sequences. The unexpected topological attraction between these valid trypanorhynch isolates, and the lecanicephalidean references, indicates a Long-Branch Attraction (LBA) artifact within the NJ framework, likely driven by accelerated nucleotide substitution rates at hypervariable ITS positions across these lineages. The third clade comprised three specific sequences (PV918672–PV918674), originating from Red Sea sampling localities, that grouped with absolute structural affinity alongside validated global C. gracilis references.

Hierarchical Clustering and Host-Driven Structure

The distance-based hierarchical clustering separated successfully isolates of this study into well-separated lineages (see Fig. 5), validating the presence of highly structured genetic diversity over a uniform single-population model. Seventeen of the examined isolates exhibited clear sequence identity matching global C. gracilis reference sequences based on BLASTN parameters. In the distance-based Ward dendrogram (Fig. 5), eleven of these isolates demonstrated tight within-group clustering with baseline C. gracilis entries. Interestingly, the remaining six verified C. gracilis isolates (IDs PV920011–PV920014, PV920017, and PV920022), clustered distinctly in an intermediate position between global Pro-teocephalidae and Bothriocephalidea reference entries, aligning close to Indian Ocean Molicola reference sequences harvested from X. gladius. While this intermediate positioning stems from high sequence variance and marker constraints across deep tax-onomic scales, this cohesive grouping highlights a powerful parasite population sub-structure dictated by host species (X. gladius) rather than strict geographic constraints, since isolates PV920013, PV920014, and PV920017 were collected from Mediterranean X. gladius. Finally, the three putative novel species isolates (IDs PV917196–PV917198) separated into an independent sub-cluster closely allied with global reference sequences labeled as Molicola sp. (GenBank Accessions: KX712333.1 and ON197571.1), providing secondary validation that these larvae represent a unique taxonomic lineage distinct from the core C. gracilis group.

Fig. 5.

Differential analysis (k-diff/split k-mer) for dissecting the molecular differences between potential new isolates and reference sequences

The research isolates (PV917196, PV917197, and PV917198) are distinguished by nucleotide motifs in the ITS region, which are represented by the identified k-mer patterns. These are referred to as “putative new species”. We focused on the k-mer present in the isolate but absent in the reference sequences (or vice versa). The summary of these findings is presented in (Fig. 6).

Fig. 6.

Binary Classification

All of the three classifiers, which we employed achieved perfect discrimination between isolates (PV917196, PV917197 and PV917198), and the 125 other Eucestoda sequences, as shown in (Table 2). Such a notable agreement across all algorithms with fundamentally different assumptions, like ensemble trees, gradient boosting, and linear models, indicates that these three isolates possess genuinely unique compositional signatures. Particularly, logistic regression, which achieved 100 %, confirming the linear separability in the k-mer space.

Multiclass Classification

The 7-class random forest classifier obtained an overall cross-validation accuracy of 82.87 % ± 9.31 %, this represnts a 5.8-fold improvement over the random baseline of 14.29 % (p < 0.001, permutation test). The performance of each class is summarized in (Table 2).

Discriminative Compositional Signatures (SHAP Interpretability)

The Shapley Additive exPlanations (SHAP) analysis of the XG-Boost binary classifier identified 15 top discriminative k-mers, as presented in (Table 3). These useful metrics showed three main differences in the composition of the novel isolates in relation to reference databases:

• 1. Enrichment of TG-rich Motifs: The frequencies of TGTG (1.87-fold), GTGT (1.80-fold), and CTGT (2.72-fold), were significantly elevated within the novel isolates (see Table 3). Such systematic enrichment suggests that the presence of tandem TG dinucleotide repeat expansions, which are likely linked to specific replication slips, or secondary structural constraints within the ITS cassette (Fig. 6A).

• Lineage-Specific Motif Absences: The motifs TGGT and GCGT that maintain mean baselines of (1.05 %) and (0.52 %) across reference Eucestoda respectively, were completely absent within all three novel isolates. This absolute drop indicates targeted lineage-specific loss, or strong negative evolutionary selection pressures (Fig. 6B).

• 2. Depletion of GC-rich Elements: The frequencies of GTGG (0.44-fold), GTTG (0.45-fold), and CTGC (0.25-fold), were significantly reduced, matching the overall lower GC-content bias, which is observed across these novel genotypes (Fig. 6B).

Our validation analysis revealed the initial and validation k-mer patterns to be totally consistent with each other. This proved that these patterns are not the results of bad analysis.

Compositional Space Visualization

The results of the PCA analysis of the 256-dimensional k-mer space are shown in (Fig. 7), and explained 37.8 % of the total variance in two components (PC1: 23.2 %, PC2: 14.6 %). The putative new species sequences clustered separately from all other taxonomic groups, in a unique region of compositional space that did not overlap with any order cluster. The Hierarchical clustering of k-mer Bray-Curtis distances confirmed tight internal grouping of putative new species (within-group distance: 0.036), with increased separation from Eucestoda (~0.33). A Bray-Curtis distance of 0.33 between groups indicates that the putative new species and the reference Eucestoda diverge in about 1/3 of the compositional space. This is a significant divergence, since in Eucestoda inter-order distances are generally between 0.05 and 0.20 for our dataset. The three novel isolates mapped in a distinct, independent region of coordinate space, with no spatial overlap or proximity with established clusters of eucestode orders, (see Fig. 7A).

Fig. 7.