Prevalence and Molecular Characterization of Ascaridoid Parasites of Philippine Decapterus species

Ethical approval for animal use was obtained from the Institutional Animal Care and Use Committee (IACUC) University of the Philippines, Los Baños (Protocol No: CAS-2019-001). Approval and sampling were coordinated to the National Stock Assessment Program (NSAP) of the Bureau of Fisheries and Aquatic Resources Region IV-A (BFAR-4A) in Lucena City Regional Office.

Balayan Bay (13.8447° N, 120.8499° E) and Tayabas Bay (13.5999° N, 121.7195° E) are among the productive fishing grounds of the Philippines. These bays are adjacent to the West Philippine Sea and are part of the Verde Island Passage, which is considered to be rich in marine biodiversity. These two bays are included in Fisheries Management Area 12 (FMA 12) designated by the Bureau of Fisheries and Aquatic Resources (BFAR) for conservation and management of fisheries in the Philippines.

Fish samples were obtained from two fish landing areas in Batangas and Lucena in May 2018. The Batangas fish landing area located in Lemery and Calatagan has most of their fish catch from the adjacent Balayan Bay. The fish landing in Dalahican Port in Lucena acquires most of the fishes and other marine resources from Tayabas Bay (Fig. 1).

Figure 1

A total of 371 fish samples were collected from Tayabas Bay and Balayan Bay and were identified based on key morphological characteristics found in FishBase (Froese and Pauly, 2019). Three species of Decapterus were identified in the two sampling sites, namely, D. maruadsi (n = 120), D. tabl (n = 130), and D. macrosoma (n = 121).

The fish samples were maintained in cold temperature and immediately transported to the Parasitology Research Laboratory, Institute of Biological Sciences, University of the Philippines Los Baños (UPLB) for dissection and processing. Before dissection, the total body length and weight of individual fish samples were measured. Sex of fish samples was also recorded based on observation of the gonads following dissection (Ohshimo et al., 2014).

The body cavity, stomach, and intestine of each dissected fish were carefully examined with the naked eye and under a stereomicroscope (Leica EZ4; Leica Microsystems Pte Ltd, Singapore) for the presence of helminth larvae (third-stage larvae of nematodes). The helminths were removed, counted, and preserved in 70% ethanol in separate Eppendorf tubes. These larvae were then cleared in a mixture of 5% glycerin in 70% ethanol for further examination under a compound microscope (Nikon, Tokyo, Japan). Parasites were identified based on morphological characteristics (Mattiucci et al., 2005; Borges et al., 2012; Hien et al., 2021). The specimens were photographed using digital microscope cameras OptixCam Summit Series SK2-5.2Xx USB (Microscope LLC, VA, USA) and Optika version 2.1 (Optika Microscopes, Bergamo, Italy). All parasites were removed and preserved in 70% ethanol, but only representative individuals of those initially identified as marine ascaridoid parasites through light microscopy were subjected to molecular processing for species identification.

Fifteen (five per fish species) morphologically similar ascaridoid larvae were randomly selected for molecular identification. Genomic DNA was extracted from ethanol-fixed midsections of individual larvae using a DNeasy^TM Tissue Kit (Qiagen, Hilden, Germany). Identification of the parasites was carried out through analysis of the small subunit ribosomal ribonucleic acid sequences. Primers were designed using the Primer-BLAST (Primer3web v. 4.1.0) tool (Ye et al., 2012) to target the 18S rRNA gene and the internal transcribed spacer 2 (ITS2) region, flanked by partial sequences of 5.8S and 28S genes. The 18S gene marker was amplified using the primer pair SSU-F (5'-ATGAGAGGGCAAGTCTGGTG-3') and SSU-R (5'-CTGTCAATCCTCACGGTGTC-3'). The amplification reactions for this gene were performed in a 25 ml volume mastermix containing 2.5 ml 10X PCR buffer, 1.875 ml 10 mM dNTPs, 2 ml 25mM MgCl₂, 1 ml 10 mM each of forward and reverse primers, 0.125 U Taq polymerase, and 2 ml 40ng/ml DNA template. The PCR conditions included an initial denaturation at 94°C for 4 min, followed by 40 cycles of denaturation at 94°C for 1 min, annealing at 55°C for 30 s, extension at 65°C for 1 min, and a final extension step at 72°C for 10 min. The 5.8S-ITS2-28S region was amplified using forward primer AscF (5'-GATCGATGAAGAACGCAGCC-3') and reverse primer AscR (5'-TTTGCAACTTTCCCTCACGG-3') with the following PCR conditions: initial denaturation of 94°C for 4 min, followed by 40 cycles of denaturation at 94°C for 1 min, annealing at 60°C for 30 s, and extension at 65°C for 1 min. A final extension step was added at 72°C for 10 min. The PCR was performed in a 25 ml volume with a similar composition of the PCR mastermix for the 18S gene amplification. The PCR products were resolved using a 1% agarose gel stained with 1% EtBr run in a horizontal gel electrophoresis system at 100 V for 30 min. Bands were visualized using a UV transilluminator. Distinct bands were sliced out from the gel and were stored in a properly labeled 2-ml microcentrifuge tube ready for gel extraction.

From the PCR amplification, 14 amplicons with a molecular weight of about 1 kb were purified using a QIAquick Gel Extraction Kit (Qiagen, Germany) and were sent to Macrogen, South Korea, for standard Sanger sequencing. Sequence assembly was performed using the Staden Package 4.0 (Staden et al., 2000). The generated consensus sequences were subjected to local alignment using the Basic Local Alignment Search Tool (https://blast.ncbi.nlm.nih.gov/Blast.cgi?PAGE_TYPE=BlastSearch) to identify the closest species matches (Altschul et al., 1990). Sequences from families Anisakidae and Raphidascaridae with closest species matches from BLAST were obtained from GenBank (Table 1). The aligned sequences of 18S and ITS2 were edited and trimmed in BioEdit Sequence Alignment Editor (Hall, 1999). The final sequence alignment for each marker was converted to a PAUP/NEXUS file in Data Analysis in Molecular Biology and Evolution (DAMBE) (Xia, 2013). The optimal model of evolution was chosen according to the Akaike Information Criterion (AIC) using jModelTest 1 (Darriba et al., 2012). In the 18S final sequence dataset, the TrNeF + G model (Tamura and Nei, 1993) was implemented involving 34 nucleotide sequences, with a total of 628 positions. And in the ITS2 final sequence dataset, the HKY + G model (Hasegawa et al., 1985) was implemented involving 33 nucleotide sequences, with a total of 460 positions. The phylogenetic relationships between nematodes were reconstructed using neighbor-joining (NJ) trees (Saitou and Nei, 1987) in PAUP version 4.0b10 (Swofford, 2002) and maximum-likelihood (ML) (Felsenstein, 1981) trees in PhyML version 3.0 (Guindon and Gascuel, 2003). The outgroup used in the reconstructed 18S phylogenetic tree was Ascaris lumbricoides (Linnaeus, 1758) from the family Ascarididae (Baird, 1853). For ITS2, Raphidascaris (Sprentascaris) andersoni (Malta et al., 2018) from the family Raphidascarididae (Hartwich, 1954) was used as an outgroup. Bootstrap resampling for both NJ and ML trees computed for 1,000 replicates and all bootstrap values >50 are shown.

The prevalence of parasite infection was calculated by dividing the number of infected fish by the total number of fish examined, multiplied by 100, and expressed as percentage. Mean intensity was calculated by dividing the total number of parasites by the total number of infected fish hosts examined. On the other hand, chi-square tests were performed to determine variations in helminth prevalence and mean intensities among fish host species and between sampling sites. This was performed using IBM SPSS Statistics v. 25 with 95% confidence interval and 0.05 level of significance. Parasite aggregation in the fish host population was measured using Poulin’s discrepancy index (D), considering the v variance/mean ratio (s²/m), and k value of negative binomial distribution, with values ranging from 0 (no aggregation) to 1 (theoretical maximum) (Estaño et al., 2021). Aggregation analysis to generate frequency distribution of the helminths was done in Quantitative Parasitology (QPweb) (Reiczigel et al., 2019).

Three hundred seventy-one Decapterus spp. (mean weight: 68.59 ± 33.72 g; mean length: 18.56 ± 2.97 cm) morphologically identified as D. maruadsi, D. tabl, and D. macrosoma were collected from Tayabas Bay and Balayan Bay adjacent to the West Philippine Sea. Of the 371 fishes examined, 80 (22%) were infected with helminth parasites initially identified as marine ascaridoid larvae through visual examination. Macroscopically, the parasites were small, thin, whitish nematodes. Under the microscope, the body of the larvae was cylindrical in shape and attenuated at both ends. Shared characteristics among individual larvae are based on general morphology of the third-stage larvae of ascaridoid nematodes, which include a smooth cuticle with transverse striations, inconspicuous lips, presence of boring tooth at the anterior end, straight gut structure composed of a muscular esophagus, irregularly shaped ventriculus, an esophagus–intestinal junction which was poorly defined in some specimens, and a tail with or without a terminal mucron (Fig. 2). Among the fish species, D. tabl showed the highest marine ascaridoid larvae infection rate (27.69%), followed by D. macrosoma (19%), and D. maruadsi (17.50%) showed the lowest, as shown in Figure 3. Each of the 80 infected Decapterus spp. hosted from one to seven anisakid larvae, except for one specimen of D. maruadsi with 39 larvae recovered. Consequently, the highest intensity was observed in D. maruadsi (4 ± 8.24), followed by D. tabl (2 ± 1.42) and then by D. macrosoma (2 ± 0.73). Both prevalence (x² [10] = 12.0, P = 0.285) and intensity (x² [6] = 6.0, P = 0.423) of infection were not significantly different across fish host species.

Figure 2

Figure 3

The variation in the prevalence and mean intensity of marine ascaridoid larvae with respect to the sampling areas is shown in Figure 4. A total of 185 individuals were collected from Tayabas Bay fish landing area in Quezon Province, while 186 individuals were from Balayan Bay in Batangas. Of the 371 fish samples, 51 (27.57%) were observed as infected from Tayabas Bay and 29 (15.59%) from Balayan Bay. While a higher prevalence is observed in Tayabas Bay than in Balayan Bay, there was a difference in the mean intensity of infection of fish samples collected from the two sites, 3 ± 5.33 in Tayabas Bay and 2 ± 1.49 in Balayan Bay. Statistical analysis revealed that prevalence (x² [5] = 6.0, P = 0.306} and mean intensity (x² [3] = 6.0, P = 0.112) did not vary significantly between the sites.

Figure 4

The marine ascaridoid parasites showed an aggregated distribution in Decapterus spp., with the value of variance significantly greater than the mean (4.79 > 0.46). The resulting variance/mean ratio showed a significant departure from Poisson (random) distribution. A negative binomial distribution parameter was also used to estimate the degree of aggregation. The resulting measurement of aggregation (k) of the negative binomial revealed that the parasites are highly aggregated in distribution, where the value of k is <1 (k = 0.19) (Fig. 5).

Figure 5

Evaluation of the identity of the marine ascaridoid larvae was performed using molecular analysis based on the 18S rRNA and the 5.8S-ITS2-28S markers. Twelve of the 15 selected larval samples generated DNA sequences that were subjected to subsequent molecular analyses. Single-gene alignment of both 18S rRNA and 5.8S-ITS2-28S regions identified two genera, Anisakis (n = 8) and Raphidascaris (n = 4). The closest species matches with the highest identity score for the ascaridoid larvae in the present study include A. typica (Deising, 1860) Baylis 1920, Anisakis simplex (Rudolphi, 1809), Raphidascaris (Ichthyascaris) lophii (Wu, 1949) Hartwich 1975, and Raphidascaris (Ichthyascaris) longispicula (Table 2). The partial sequences of 18S rRNA (628 bp) generated alignments, which presented genetic matches of the eight samples with zoonotic Anisakis simplex complex (99.20%–99.68%) and four others with Raphidascaris longispicula (99.52%–99.84%). Aligning the sequences using the 5.8S-ITS2-28S region (460 bp), the ascaridoid larvae were found to be genetically similar with A. typica (98.61%–99.88%) and Raphidascaris lophii (99.18%–99.79%).

Pairwise distance comparison (Table 3) indicates low genetic variation within species for both 18S and 5.8S-ITS2-28S markers. However, sequence variation between species and between genera is at lowest for 18S rRNA gene, with P-distance ranging from 0 to 0.0425. Hence, an overlap for the pairwise comparison is observed between the paired entities for both intra- and interspecific categories (Fig. 6). A more pronounced overlap, however, is found among the 18S sequences of the GenBank and sample isolates. Meanwhile, a higher variability in the sequences of the different species and genera is observed for the 5.8S-ITS2-28S region (P-distance 0.004 to 1.877).

Figure 6

Phylogenetic reconstruction using the ML and NJ methods based on the 18S rRNA gene showed two major clusters of the currently studied ascaridoid larvae. Four samples were resolved in close association with Raphidascaris longispicula and Hysterothylacium species, with a bootstrap percentage of 90/95.6% ML/NJ bootstraps (Fig. 7). Meanwhile, the other eight samples were grouped together, forming an external cluster with the zoonotic A. simplex species complex (85/87.7% ML/NJ bootstraps). In this reconstruction, A. typica species was not represented due to the unavailability of GenBank accessions.

Figure 7

The phylogeny based on the partial sequence of the 5.8S-ITS2-28S rRNA region (Fig. 8) showed similar grouping of the four Raphidascaris sp. larval samples as an external sister clade of R. lophii with strong nodal support for both ML and NJ trees (96/100% bootstraps). The inclusion of A. typica reference accessions generated similar clustering of Anisakis sp. samples with full support (100/100% ML/NJ bootstraps). It can be inferred that all Anisakis sp. samples are genetically similar to A. typica (P-distance = 0.0077422) of cosmopolitan distribution.

Figure 8