Differentiating paramphistome species in cattle using DNA barcoding coupled with high-resolution melting analysis (Bar-HRM)

Paramphistomosis is caused by paramphistome or amphistome parasites, including Fischoederius elongatus, Gastrothylax crumenifer, Orthocoelium parvipapillatum, and Paramphistomum epiclitum. The control and prevention of these parasite outbreaks are difficult because of the wide occurrence of these species. Besides, the clinical manifestations and their egg characteristics are similar to those of other intestinal flukes in the paramphistome group, leading to misdiagnosis. Here, we employed DNA barcoding using NADH dehydrogenase (ubiquinone, alpha 1) (ND1) and cytochrome c oxidase subunit I (COI), coupled with high-resolution melting analysis (Bar-HRM), for species differentiation. As a result, ParND1_3 and ParCOI4 resulted in positive amplification in the paramphistomes and Fasciola gigantica, with significantly different melting curves for each species. The melting temperatures of each species obtained clearly differed. Regarding sensitivity, the limit of detection (LoD) for all species of paramphistomes was 1 pg/µl. Our findings suggest that Bar-HRM using ParND1_3 is highly suitable for the differentiation of paramphistome species. This approach can be used in parasite detection and epidemiological studies in cattle.

Traditionally, parasite diagnosis in ruminant livestock is carried out by egg observation in feces under light microscopy. Furthermore, the differentiation of species of the families Paramphistomidae and Gastrothylacidae based on egg morphology is difficult because of similarities in morphological characteristics such as size, color, transparency of the eggshell, and inside components (Chryssafidis et al. 2015; Nak-on and Chontananarth 2020). The adult stages of both paramphistome families are required for species identification under a scanning electron microscope (SEM) to investigate the inner and outer characteristics to perform species taxonomy according to the keys of Yamaguti (1958), Eduardo (1982aEduardo ( , b-1985, Sey (1991), and Jones et al. (2005). Detection of fluke species from cattle might affect the choice of parasite treatment. For example, albendazole is more in the elimination of liver flukes than of rumen flukes (Arias et al. 2013). Hence, species taxonomy using fluke eggs, based on morphological characteristics, is still insufficient. Molecular approaches have been used to facilitate species differentiation/identification for unidentified specimens due to incomplete or immature stages of the organisms.
For species differentiation/identification, the DNA barcodes, the unique sequence of the short, standardized barcode (200-800 bp) which are capable of representative of the specific species, has recently received attention (Hebert et al. 2004;Ward et al. 2009). Currently, DNA barcodes are used extensively in parasites on the mitochondrial DNA, including NADH dehydrogenase subunit 1 (ND1), cytochrome c oxidase subunit 1 (COI), cytochrome b (cytb), and ribosomal RNA (Wood 2018) as well as nuclear DNA including 5.8S rRNA and internal transcribed spacer (ITS) (Khan et al. 2020). For example, in our previous study, COI was successfully applied to distinguish the species of rumen flukes including Gastrothylax, Carmyerius, Fischoederius, Paramphistomum, Explanatum, and Fasciola from egg specimens in feces (Anucherngchai et al. 2020). Although DNA barcoding enables species identification, it is complex, timeconsuming, and costly. To solve this issue, high-resolution melting analysis (HRM) has been adopted to discriminate the differences in sequences obtained from different species based on the melting behavior (DNA denaturation with increasing temperature). DNA barcodes contain the variation in sequences among different species. When HRM is used for DNA barcodes, it would provide the unique behavior of DNA melting shape. The melting shape and melting temperature obtained from DNA barcodes coupled with HRM (Bar-HRM) can be use as representative of the certain species. Bar-HRM has been employed to accelerate the species differentiation in various parasite genera, such as Schitosoma species (Li et al. 2015), liver flukes, Fasciolosis hepatiga and Fasciolosis gigantiga (Hosseini-safa et al. 2019), and Echninostoma species (Buddhachat and Chontananarth 2019). The use of Bar-HRM can eliminate the sequencing process, thereby shortening the process of species differentiation and making it more economically viable.
In this study, we aimed to develop a molecular tool for accommodating species differentiation among different species of rumen flukes using mini-DNA barcodes, ND1 and COI, combined with high-resolution melting analysis, the so-called Bar-HRM, which offers accuracy, robustness, and speed in fluke diagnosis and epidemiology.

Specimens
Rumen flukes of four paramphistome species, namely, Fischoderirus elongatus, Gastrothylax crumenifer, Orthocoelium parviparpillatum, and Paramphistomum epiclitum, were collected from cattle rumens (Bubalus bubalis for gastrothylacids and Bos taurus for paramphistomids) obtained from slaughterhouses. Also, other flukes which are both closely and distantly related species were included, such as Fasciola gigantica, Echinostoma revolutum, Hypoderaum conoideum, Thapariella sp., Clinostomum sp., Raillietina echinobothrida, and Ascaridia galli. All adult flukes of each species in this study were previously used to confirm the species based on morphological features under a light microscope and a scanning electron microscope, according to taxonomic keys of Eduardo (1982aEduardo ( , b-1985 and Sey (1991) for paramphistome species and Jones et al. (2005) and Yamaguti (1958) for all specimens. In addition, paramphistome eggs were collected from incubation of adult paramphistomes in PBS solution. Egg specimens that fell to the bottom of the plate overnight were observed under a stereo microscope and collected into a microcentrifuge tube which contained absolute ethanol. We used the DNA barcodes ND1 (Bowles and McManus 1993) and COI (Morgan and Blair 1998) loci to confirm species (Supplementary file 1), and the sequences obtained from Fi. elongatus, G. crumenifer, O. parviparpillatum, and P. epiclitum were deposited in GenBank for ND1 as OP939413-OP939415 except for Fi. elongatus, respectively, and COI as OP928154-OP928157, respectively (Supplemental file 2). The animal ethic of this experiment was approved by Srinakharinwirot University under SWE-A-022-2563.

DNA extraction
Genomic DNA was extracted from fresh rumen fluke tissue samples using a Genomic DNA Isolation Kit (Tissue) (PureDireX, Taiwan), according to the manufacturer's protocol. All tissues were incubated at 60 °C overnight (16-18 h) for complete digestion. The extracted DNA was quantitatively measured by a NanoDrop LITE spectrophotometer (Thermo Scientific, USA) and subjected to 1.5% agarose gel electrophoresis for quality and quantity checks. The DNA concentration was adjusted to 10 ng/µl, and the DNA was stored at − 20 °C until further analysis.

Primers design for HRM
The mitochondrial gene sequences, ND1 and the COI or cox, of the studied flukes were retrieved from GenBank as shown in Supplementary file 3. These sequences of each gene were obtained for multiple alignment in MUTALIN (https:// www. multa lin. toula ouse. intra. fr/ multa lin), and subsequently, we evaluated the nucleotide diversity ( ) values of both regions by spiting window at 20 bp, which separated two groups: (i) the paramphistome group and (ii) all fluke species under BioEdit (Hall 1999). Subsequently, the primer sets were selected at low for the paramphistome only group, with a large difference of values between the paramphistome only group and all fluke species and with an amplicon length below 300 bp. Finally, we obtained seven primer sets, including three primer sets for ND1 and four primer sets for COI (Table 1).

Primer selection for DNA amplification
To find the suitable primer to amplify the DNA of paramphistome group, PCR was performed with four paramphistome species, namely, Fi. elongatus, G. crumenifer, O. parvipapillatum, and P. epiclitum, in a final volume of 25 µl, containing 12.5 µl of 2 × DreamTaq PCR Master Mix (Thermo Fisher Scientific, USA), 0.3 µM of each primer, the addition of 1 mM of MgCl 2 , and 1 µl of DNA (10 ng/ µl), with the addition of sterile water to obtain a volume of 25 µl. The reaction mixture was heated for pre-denaturation at 95 °C for 5 min, followed by 40 cycles of amplification at 95 °C for 30 s for denaturation, 55 °C for 30 s for annealing, 72 °C for 30 s, with a final extension step at 72 °C for 7 min after 40 cycles using a MiniAmp Plus Thermal Cycler (Thermo Fisher Scientific, USA). The PCR products obtained were analyzed by 1.5% agarose gel electrophoresis at 100 V for 25 min and stained by EtBr to visualize the DNA amplicon under an ultraviolet (UV) transilluminator.

Specificity determination of Bar-HRM
After attaining the suitable primers, ParND1_3 and Par-COI4, they were used for performing HRM analysis to test their specificity. The four paramphistome and other fluke species were investigated via HRM in a total volume of 25 µl, using a LightCycler®480 Instrument II (Roche Dianostics, Japan). Each reaction mixture consisted of 2.5 µl of 1 × PCR buffer, 1 unit of Platinum® Taq DNA polymerase, 0.2 mM dNTP mixture, 2 mM SYTO®9 1 µl (Invitrogen, USA), 1.5 mM MgCl 2 , 0.2 µM of each primer, and 1 µl of DNA (10 ng/µl), with the addition of sterile water to obtain a volume of 25 µl. The reaction mixture was heated for pre-denaturation at 95 °C for 5 min, followed by 40 cycles of amplification at 95 °C for 30 s for denaturation, 55 °C for 30 s for annealing, 72 °C for 30 s, with a final extension step at 72 °C for 7 min after 40 cycles. The HRM analysis was performed using an initial hold at 95 °C for 1 min and ramping from 65 to 95 °C, with 0.1 °C per step. The acquired fluorescence data were analyzed for the normalization plot, the difference plot, and the melting peak to obtain the melting profile for each species amplified by Bar-HRM, using the ParND1_3 and COI4 primer sets for species differentiation. The products were subjected to 1.5% agarose gel electrophoresis at 100 V for 25 min and stained by EtBr for visualization under a UV transilluminator to check whether the amplicon size obtained met the expected size. The PCR products of the four paramphistome obtained from Bar-HRM using the ParND1_3 and COI4 were sequenced

Species identification of fluke eggs by Bar-HRM
Fluke eggs were collected from overnight incubation of adult paramphistomes, collected into five tubes, and the DNA was extracted using a Genomic DNA Isolation Kit (Tissue) (PureDireX, Taiwan). Single-blinded species identification was performed via Bar-HRM using ParND1_3; here, we investigated fluke eggs along with four paramphistomes and other samples in a final volume of 25 µl. Each reaction mixture contained 2.

Statistical analysis
The melting temperature of each species derived from Bar-HRM for both primers was tested for significant differences using ANOVA, and multiple ranking tests for comparison among species were performed using Duncan's test at p < 0.05. The detection limit of Bar-HRM using ParND1_3 for species discrimination across paramphistome species was determined as the least concentration showing positive DNA amplification. Besides, to access the precision within a run and among run, intra-assay and inter-assay variability were calculated based on the coefficient of variation (CV). The intra-assay CV was calculated by performing Bar-HRM of each parasite in two or three replicates each, on a single plate. The inter-assay CV was calculated by performing different species in two different plates. The percent of coefficient of variation (%CV) for the intra-assay and inter-assay CV was determined.

Designing and selecting suitable primers
In this study, there were seven primer sets designed from two regions of mitochondrial ND1 and COI genes; the DNA binding sites of these primers were specifically selected for paramphistomes based on multiple alignment and nucleotide diversity ( ). Most of the primers were selected for the DNA region, with low , along with large differences in values among and within paramphistomes and all fluke species (Fig. 1). We obtained three primer sets from ND1 (Fig. 1A) and four primer sets from COI (Fig. 1B). Out of two primer sets from either ND1 and COI, ParND1_3 and ParCOI4 enabled the DNA amplification of four species of paramphistomes (Fi. elongatus, G. crumenifer, O. parvipapillatum, and P. epiclitum) (Fig. 1C, D), although ParCOI4 appeared to produce a fade band in P. epiclitum (Fig. 1D). Both primers were used for fluke species discrimination of paramphistomes by Bar-HRM.

Species discrimination among paramphistomes by Bar-HRM
To evaluate the feasible use of Bar-HRM using either ParND1_3 or ParCOI4 for species differentiation of flukes in the paramphistome group, four species were investigated. As a result, Bar-HRM using ParND1_3 showed that the melting shapes of four paramphistome species differed in the normalization plot ( Fig. 2A) and the melting curve (Fig. 2B). The melting temperatures obtained from Bar-HRM using ParND1_3 of four species followed the order P. epiclitum (82.6 °C) > G. crumenifer (81.7 °C) > Fi.  Table 2); the T m of each species differed by about 1 °C in average. However, the T m of P. epiclitum showed a larger variation (SD = 0.6) than that of other species (Table 2). Besides, the intra-and inter-assay CV displayed less than 0.5% for HRM using both ParND1_3 and ParCOI4, indicating the reproducibility of hands-on setting and plate-to-plate consistency (Table 2). For HRM using ParCOI4, there were three species yielding positive values, namely, G. crumenifer, Fi. elongatus, and O. parvipapillatum (Fig. 2C, D). The normalization plot (Fig. 2C) and melting peaks (Fig. 2D) exhibited a small difference in melting shapes among three species, resulting in little discrepancy in T m . We noted that G. crumenifer showed higher values, followed by Fi. elongatus and O. parvipapillatum (Table 2).
Bar-HRM using ParND1_3 and ParCOI4 was determined for cross-reactions with other fluke species both closely and distantly related (Fig. 3). For Bar-HRM using ParND1_3, there were three spices positive of other flukes, excluding paramphistome species such as F. gigantica, Thapariella sp., and Clinostomum sp., with melting profiles similar to that of P. epiclitum (Fig. 3A, B). In contrast, the HRM using ParCOI4 gave the positive for F. gigantica (Fig. 3C, D).

Sensitivity determination for species identification by HRM using ParND1_3
Based on the specificity test, Bar-HRM using ParND1_4 showed a great ability to distinguish paramphistome species. We used this primer to determine the sensitivity on each species of four paramphistomes in duplicate. The results showed that the sensitivity was different among the studied species. The detection limit (showing positive for duplicate) for G. crumenifer was 10 −4 ng/µl, followed by those for Fi. elongatus (10 −2 ng/µl), O. parvipapillatum (10 −1 ng/µl), and P. epiclitum (10 ng/ µl) ( Fig. 4 and Table 3).

Bar-HRM-based species identification from fluke eggs
In the present study, fluke eggs were collected from overnight incubation of adult paramphistomes; the isolated DNA was single-blind identified by Bar-HRM using ParND1_3. The melting profiles of five samples showed three different patterns: (i) egg 2 showed a pattern similar to that of O. parvipapillatum, (ii) egg 1 was similar to Fi. elongatus, and eggs 3-5 could be grouped into the melting shape of P. epiclitum (Fig. 5). As shown in Table 3, the T m of each egg was considered to fit the T m of paramphistome species; for example, the T m of egg1 was 79.9, matching the T m of O. parvipapillatum, and the T m of eggs 3-5 was approximately 82.0, close to that of P. epiclitum (range of T m : 82.0-88.6) ( Table 4).

Discussion
Cattle infected with paramphistomes, the rumen flukes show paramphistomosis which is the manifestation of symptoms such as diarrhea, green foamy droppings with a strong stench, anorexia, and weight loss (Boray 1969). Therefore, the morbidity and mortality of this economically important livestock lead to economic losses, partly because of insufficient knowledge of their biology and because of inadequate farm management. The establishment of accuracy, sensitive, and speed to facilitate the species differentiation especially in the egg stage is a pivotal task. This study is the first report of the use of DNA barcodes, ND1 and COI, conjugated with HRM (Bar-HRM), for species differentiation in paramphistomes, including the species Fi. elongatus, G. crumenifer, O. parvipapillatum, and P. epiclitum based on their melting profiles (melting and melting temperature). Bar-HRM using ParND1_3 showed a high ability to discriminate four paramphistomes and some other fluke species, such as F. gigantica. The effective amount of DNA to amplify by Bar-HRM using ParND1_3 must be above 10 ng/ µl for all four species. We could also verify this approach to identify paramphistome species from eggs as the species Fi. elongatus, O. parvipapillatum, and P. epiclitum were successfully identified.
The DNA barcodes ND1, ND2, and COI have been widely used for studying species taxonomy and evolution in platyhelminths (Van Steenkiste et al. 2015;Yang et al. 2016;Buddhachat and Chontananarth 2019;Semyenova et al. 2006;Amer et al. 2011;Chaudhry et al. 2016;Anucherngchai et al. 2020). For example, Van Steenkiste et al. (2015) showed that the sequence of the COI gene at 1000 bp Table 2 Melting temperature (T m ) and the coefficient of variation (CV) of intra-assay and the inter-assay obtained from Bar-HRM using ParND1_3 and ParCOI4 a, b, c, d represent the significant difference by ANOVA and multiple comparison by Duncan at p-value < 0.05. Not determined is abbreviated as ND. Standard deviation is abbreviated as SD  F. elongatus, F. cobboldi), significant nucleotide differences of several genes were observed, ranging from 6.0-13.2%; the ND gene family seemed to have a higher nucleotide variability than the COI (cox) gene family (Watthanasiri et al. 2021).
Both ND1 and COI have been used for facilitating species discrimination in various parasites such as Echinostoma (Buddhachat and Chontananarth 2019), Fasciola (Semyenova et al. 2006;Amer et al. 2011;Chaudhry et al. 2016), and rumen flukes (Anucherngchai et al. 2020). Based on multiple alignment and nucleotide diversity, the DNA binding sites of ParND1_3 and ParCOI4 primer have different nucleotides among and within the paramphistome group and all flukes, possibly contributing to the specificity of primer sets to paramphistome species and closely related species. The nucleotides within the amplicon amplified of both ParND1_3 and ParCOI4 showed a high diversity ( > 0.1), interspersed in few regions, contributing to the melting behavior of double-stranded DNA and to species differentiability among the paramphistome group for the four species. Thus, it was not surprising that these two primer sets are promising for the identification of paramphistome species.
Furthermore, we tested the sensitivity of this approach for DNA amplification, resulting in different detection limits for different species; G. crumenifer showed the highest sensitivity, followed by Fi. elongatus, O. parvipapillatum, and P. epiclitum. We postulate that there were several reasons for this phenomenon: (i) the variation in nucleotides on the primer binding site due to some portion of DNA sequence on primers contained the inter-specific nucleotide variation for different paramphistome species, (ii) the different amounts of mitochondrial copy numbers in each species for even cell type and stage (Furda et al. 2012;Rooney et al. 2015), and (iii) the quality of the DNA isolated from different species  Table 3 The sensitivity of Bar-HRM using ParND1_3 for species identification among four species of paramphistome species Sample DNA concentration (ng/µl) 10 10 0 10 −1 10 −2 10 −3 10 −4 10 −5 10 −6 10 −7 Fi. elongatus because of the chemical composition and difficulties in tissue breakdown, resulting from the tightly connective tissues of some parasites. However, we found that the DNA concentration for Bar-HRM using ParND1_3 should be higher than 10 ng/µl for successful species discrimination. In addition to Bar-HRM using ParND1_3 and ParCOI4, DNA amplification for four species of the paramphistome group enabled us to yield amplicons for F. gigantica. Yang et al. (2016) demonstrated that the genus Fasciola (F. gigantica and F. hepatica) is a sister clade of paramphistomes, based on the phylogenetic tree constructed from concatenated amino acid sequence data representing 12 proteincoding genes. When the oligonucleotide sequences of both ParND1_4 and ParCOI3 were checked for cross-binding with ND1 and COI of F. gigantica, the primers could bind to sequences of F. gigantica with more than 90%, completely binding at the 3'end with the first 10 positions. With respect to this, it was not surprising that these primer sets could produce the amplicon from the genus F. gigantica. However, although the melting profile of F. gigantica obtained from Bar-HRM using ParND1_3 was similar to that of P. epiclitum, Bar-HRM using ParCOI4 did not seem to amplify P. epiclitum, and the melting profile of F. gigantica was close to those of Fi. elongatus and G. crumenifer. Thereby, combining the melting profiles obtained from the two primer sets (ParND1_3 and ParCOI4) or multilocus detection can be used to correctly identify F. gigantica. Similarly, differentiating Echinostoma revolutum from other species of the genus Echinostoma, relying on the melting profile obtained from HRM, required two loci of ND1 and COI for correct identification (Buddhachat and Chontananarth 2019). Buddhachat et al. (2020) developed multiplex HRM (mHRM) for the simultaneous detection of canine hemoparasites (Babesia sp., Hepatozoon canis, Ehrlichia canis, and Anaplasma platys) at a single reaction.
To evaluate the performance of Bar-HRM for cattle rumen fluke detection, the eggs were identified by singleblinded Bar-HRM. The melting profiles obtained for five fluke eggs appeared to be similar to those of three paramphistome species, namely, Fi. elongatus (n = 1), O. parvipapillatum (n = 1), and P. epiclitum (n = 3). This indicates that Bar-HRM can be adopted for multiple species identification in the paramphistome group with a single reaction. Traditionally, parasite diagnosis in ruminant livestock is carried out by egg observation in feces under light microscopy, using different methods such as direct wet smear (Thienpont et al. 1986), the floating technique (Thienpont et al. 1986), and concentration techniques, which consist of two approaches: the sediment technique and the formalin-ether technique (Anuar et al. 2013). The latter is a suitable method for detecting paramphistome operculated eggs because the egg is heavier than the solution of formalin-ether, resulting in sedimentation (Truant et al. 1981). However, parasite determination under the microscope is time-consuming and requires expertise to discriminate the fluke eggs based on morphological features (Chryssafidis et al. 2015). The use of Bar-HRM using ParND1_3 for species discrimination among paramphistomes can address the restriction of microscopic examinations, contributing to multiplex detection of several species with specificity, sensitivity, and speed.
Although Bar-HRM using ParND1_3 exhibited a great alternative tool for aiding species differentiation in paramphistomes, the restriction of this method might occur to practical use, for example, the DNA variation in the amplicon causing the shift of melting shape and melting temperature especially the substitution of G or C to A or T and vice versa, resulting in 0.5-1 °C shift (Venter et al. 2001). This may cause the incorrect species prediction. We  suggest that the sample with shifted melting shape should be confirmed by sequencing or performing Bar-HRM using other barcodes like COI. Besides, the quantitation of egg by Bar-HRM should be done for further study since the number of parasite egg found in feces can indicate the severity of the disease and suitable treatments. In summary, Bar-HRM using ParND1_3 is a promising approach to accelerate species identification in paramphistomes, with high specificity and multiple species discriminations at a single reaction. The amount of DNA used for this method should not be less than 10 ng to ensure positive amplification. For more accuracy and reliability, we suggest that HRM using multiple DNA loci or mHRM should be used for species identification in a single test. We believe that this method is a powerful tool to determine species taxonomy and epidemiology and for detection in groups comprising various species, in particular in the field of parasite epidemiology.