Phylogeographic analysis of Fasciola hepatica populations across the Indian subcontinent using mitochondrial and nuclear markers

Introduction

The earliest written descriptions of liver rot in sheep-later attributed to Fasciola hepatica infection-appear in European veterinary texts from the 17th century, and for much of the subsequent three hundred years, fasciolosis was regarded primarily as a disease of temperate climates [1] . This Eurocentric view began to shift in the 1990s when molecular tools revealed F. hepatica populations coexisting with F. gigantica across large swathes of Asia and Africa, sometimes within the same host animal and even hybridising to produce intermediate forms [2] . The Indian subcontinent presents a particularly complex epidemiological landscape. F. hepatica thrives in the cooler highland regions of northern India, Nepal, and northern Thailand, while F. gigantica dominates the tropical lowlands. Between these extremes lies a broad overlap zone where both species-and their hybrids-circulate among cattle, buffalo, goats, and sheep [3] . Understanding how F. hepatica populations are structured across this heterogeneous terrain is not merely an academic exercise: anthelmintic resistance alleles can spread through parasite gene flow, and the efficacy of any future vaccine will depend on matching circulating genotypes to the immunogenic targets used in vaccine design [4] . Mitochondrial markers, particularly the NADH dehydrogenase subunit 1 (nad1) and cytochrome c oxidase subunit 1 (cox1) genes, have become standard tools for Fasciola phylogeography because they evolve relatively rapidly, are maternally inherited (avoiding recombination), and have large reference databases in GenBank [5] . Nuclear markers such as the internal transcribed spacer 2 (ITS-2) of ribosomal DNA provide a complementary, biparentally inherited perspective, though their lower mutation rate means they capture deeper phylogenetic structure rather than recent population dynamics [6] . Despite the recognised importance of the Indian subcontinent as a contact zone between Fasciola species, no published research has simultaneously genotyped F. hepatica populations across Thailand, India, Nepal, Sri Lanka, and Bangladesh using both mitochondrial and nuclear markers in a single comparative framework. Existing data are fragmented into single-country or single-marker reports that cannot be directly compared owing to differences in primers, sequencing platforms, and analytical approaches [7] . This research aimed to fill that gap by conducting a multi-marker phylogeographic analysis of F. hepatica across eight geographic regions spanning five countries, with three objectives: to quantify haplotype diversity and population differentiation, to reconstruct demographic history and dispersal patterns, and to identify implications for parasite control strategies.

Reagents and Chemicals

Genomic DNA was extracted using the DNeasy Blood and Tissue Kit (Qiagen, Hilden, Germany). PCR amplification was performed using GoTaq Green Master Mix (Promega, Madison, WI, USA). Primers for nad1 (forward: 5'-TGTTTTGATTTTGCCGAAGG-3'; reverse: 5'-AAACCAACCACTAACTAATCACT-3'), cox1 (forward: 5'-TGTTTAGTTTTTGGGCATCC-3'; reverse: 5'-AGGCCAATGAAAGAGCCAAC-3'), and ITS-2 (forward: 5'-GGTACCGGTGGATCACTCGG-3'; reverse: 5'-GGGATCCTGGTTAGTTTCTTTTCCTCCG-3') were synthesised by Macrogen Inc. (Seoul, South Korea). Agarose (molecular biology grade), ethidium bromide, and 100 bp DNA ladder were obtained from Bio-Rad Laboratories (Hercules, CA, USA). Proteinase K, RNase A, phosphatebuffered saline tablets, and molecular-grade ethanol were purchased from Sigma-Aldrich (St. Louis, MO, USA). All PCR reagents were stored at-20 °C until use.

Materials and Methods

Materials

Adult F. hepatica specimens were collected from the bile ducts of naturally infected cattle (n = 187) and water buffalo (n = 108) at abattoirs and slaughter facilities across eight geographic regions: northern Thailand (Chiang Mai and Chiang Rai provinces), central Thailand (Nakhon Pathom), southern Thailand (Nakhon Si Thammarat), eastern India (West Bengal and Odisha), western India (Maharashtra and Gujarat), Nepal Terai (Chitwan and Rupandehi districts), Sri Lanka (Central and North Central provinces), and Bangladesh (Mymensingh and Sylhet divisions). Collection took place between March 2022 and November 2023. Ethical approval was obtained from the Chiang Mai University Animal Use Committee (Protocol CMU-AUC-2022-016) and from institutional review boards at collaborating institutions in each country. At each collection site, 30-50 flukes were collected, washed in phosphate-buffered saline, morphologically identified as F. hepatica based on body length-to-width ratio, cephalic cone morphology, and shoulder prominence [2] , and preserved individually in 95% ethanol at 4 °C. A total of 295 specimens confirmed as F. hepatica by morphological and molecular criteria were included in the final analysis. Specimens showing intermediate morphology suggestive of F. hepatica × F. gigantica hybrids (n = 23) were excluded.

Methods

Genomic DNA was extracted from approximately 20 mg of fluke tissue (posterior end to avoid reproductive organs) using the DNeasy Blood and Tissue Kit following the manufacturer's protocol. Partial fragments of nad1 (approximately 535 bp), cox1 (approximately 440 bp), and ITS-2 (approximately 361 bp) were amplified by PCR in 25 µL reactions containing 12.5 µL GoTaq Green Master Mix, 1 µL each primer (10 µM), 2 µL template DNA, and 8.5 µL nuclease-free water. Thermocycling conditions were: initial denaturation at 95 °C for 5 min; 35 cycles of 95 °C for 30 s, 55 °C (nad1 and cox1) or 58 °C (ITS-2) for 30 s, and 72 °C for 45 s; and final extension at 72 °C for 7 min. PCR products were visualised on 1.5% agarose gels, purified using ExoSAP-IT (Applied Biosystems), and sequenced bidirectionally on an ABI 3730xl DNA Analyser (Macrogen Inc., Seoul). Sequences were assembled and edited in Geneious Prime v2023.2, aligned with ClustalW, and trimmed to uniform lengths. Haplotype diversity (Hd), nucleotide diversity (π), and neutrality tests (Tajima's D, Fu's Fs) were calculated in DnaSP v6.12. Pairwise FST values and analysis of molecular variance (AMOVA) were computed in Arlequin v3.5. Median-joining haplotype networks were constructed in PopART v1.7. Maximum likelihood phylogenetic trees were built in IQ-TREE v2.2 with 1,000 ultrafast bootstrap replicates, using the best-fit substitution model selected by ModelFinder. All novel sequences were deposited in GenBank.

Extraction Protocol

Each fluke specimen was removed from ethanol storage and rinsed three times in sterile phosphate-buffered saline to remove residual fixative. Approximately 20 mg of tissue was excised from the posterior body region using a sterile scalpel and placed in a 1.5 mL microcentrifuge tube. Tissue was lysed overnight at 56 °C in 180 µL ATL buffer with 20 µL Proteinase K (20 mg/mL). The following morning, 4 µL RNase A (100 mg/mL) was added and samples were incubated at room temperature for 2 minutes before proceeding with the DNeasy spin-column protocol. DNA was eluted in 100 µL AE buffer (prewarmed to 70 °C) and quantified using a NanoDrop 2000 spectrophotometer. DNA purity was assessed by A260/A280 ratio (accepted range 1.7-2.0) and integrity was confirmed by electrophoresis on 0.8% agarose gels. Extracted DNA was stored at-20 °C in 20 µL aliquots to minimise freeze-thaw cycles. Extraction blanks (no tissue) were included in every batch of 20 samples to monitor contamination.

Quality Control

Multiple quality control measures were implemented throughout the molecular workflow. Each PCR batch included a positive control (DNA from a morphologically confirmed F. hepatica specimen with a known GenBank-deposited sequence), a negative control (nuclease-free water), and an extraction blank. Bidirectional sequencing was performed for all specimens, and forward and reverse chromatograms were manually inspected for ambiguous base calls. Sequences with more than 2% ambiguous positions were re-amplified and resequenced. Species identity was confirmed molecularly by BLAST comparison against GenBank, requiring greater than 99% identity to published F. hepatica reference sequences (GenBank accession AF216697 for nad1, AF216698 for cox1). Specimens showing less than 97% identity or intermediate identity between F. hepatica and F. gigantica references were classified as potential hybrids and excluded. Repeatability of haplotype assignments was verified by independently re-extracting and re-sequencing 10% of samples (n = 30); all repeat sequences matched the original haplotype assignment.

Results

Sequencing yielded clean, aligned datasets of 521 bp for nad1 (295 sequences), 428 bp for cox1 (295 sequences), and 348 bp for ITS-2 (289 sequences; 6 failed amplification). Across the three markers, 200 unique haplotypes were identified in total: 87 for nad1, 72 for cox1, and 41 for ITS-2. No insertions, deletions, or premature stop codons were detected in the protein-coding mitochondrial genes, confirming amplification of functional mitochondrial DNA rather than nuclear pseudogenes. Eastern India (West Bengal and Odisha) consistently harboured the highest haplotype diversity: Hd = 0.921 (nad1), 0.873 (cox1), and 0.714 (ITS-2). Sri Lanka showed the lowest diversity across all markers (Hd = 0.768, 0.713, and 0.561 respectively), consistent with an isolated island population subject to founder effects and genetic drift. Nucleotide diversity followed the same geographic gradient: eastern India showed the highest π values (0.0187 for nad1) and Sri Lanka the lowest (0.0086). Tajima's D was significantly negative for eastern India (D =-1.847, p<0.05) and Bangladesh (D =-1.692, p<0.05), suggesting historical population expansion in these regions.

Comprehensive Interpretation

Taken together, the results reveal a clear east-to-west and south-to-east gradient in F. hepatica genetic diversity across the Indian subcontinent, with eastern India as the apparent diversity hotspot. The star-like topology of the haplotype network, centred on eastern Indian haplotypes, combined with significantly negative Tajima's D and Fu's Fs values in this region, supports a scenario of historical demographic expansion from an eastern Indian origin. Thai populations, while genetically less diverse, showed low internal differentiation (FST = 0.028-0.051), suggesting recent or ongoing gene flow among Thai provinces-likely mediated by domestic cattle trade. Sri Lanka's genetic isolation (FST = 0.127-0.238 against all mainland populations) is consistent with limited livestock importation and the island's geographic separation.

Discussion

The central finding of this research is that F. hepatica populations across the Indian subcontinent are genetically structured but not fully isolated, with eastern India serving as the likely centre of diversity and a source for secondary dispersal into Southeast Asia. This finding has three layers of significance: for understanding parasite biogeography, for predicting anthelmintic resistance spread, and for designing regionally appropriate control interventions.

The high haplotype diversity in eastern India (Hd = 0.921 for nad1) agrees with the proposal by Mas-Coma et al. [3] that the Ganges-Brahmaputra floodplain, with its dense cattle populations, year-round surface water, and abundant Lymnaea snail habitats, provides ideal conditions for maintaining large, genetically diverse fluke populations. Our FST values between Thai populations (0.028-0.051) are comparable to the 0.02-0.06 range reported by Itagaki et al. [8] for F. hepatica in Japanese prefectures, reinforcing the pattern that withincountry livestock movement homogenises parasite populations over moderate geographic distances.

A notable disagreement with existing literature concerns Sri Lanka. Rajapakse et al. [9] suggested that Sri Lankan F. hepatica populations were recently introduced from India based on morphological similarity, but our molecular data show moderate to high FST values (0.127-0.238) and reduced haplotype diversity, more consistent with an older colonisation event followed by genetic drift in relative isolation. The positive Tajima's D value for Sri Lanka (0.342) further supports a stable or contracting population rather than a recent expansion. The discrepancy between mitochondrial and nuclear marker resolution merits comment. ITS-2 showed consistently lower FST values and fewer haplotypes than nad1 or cox1, which is expected given its slower mutation rate and biparental inheritance. However, the general geographic patterns were concordant across all three markers, strengthening confidence that the observed structure reflects genuine demographic history rather than marker-specific artefacts [6] .

For parasite control, these findings carry two practical messages. First, the moderate gene flow among mainland populations means that anthelmintic resistance alleles arising in one region-such as the triclabendazole resistance increasingly reported in Indian livestock [10] -could spread to neighbouring countries within years to decades via livestock trade. Second, any vaccine targeting F. hepatica in this region should account for the genetic variation documented here; a vaccine designed against European genotypes may show reduced efficacy against South Asian variants with divergent epitope sequences [4] .

Conclusion

The key takeaway from this phylogeographic analysis is that Fasciola hepatica populations across the Indian subcontinent are genetically structured along a diversity gradient centred on eastern India, with Thai populations showing recent internal homogenisation and Sri Lankan populations displaying genetic isolation consistent with founder effects. These patterns reflect the interplay of livestock trade, geographic barriers, and demographic history in shaping parasite population structure. The implications are both scientific and practical. Scientifically, eastern India emerges as a priority region for deeper genomic investigation, including whole-genome sequencing of F. hepatica to resolve fine-scale population dynamics and identify signatures of selection associated with drug resistance. Practically, livestock trade corridors between India, Nepal, and Bangladesh represent conduits for resistance gene dispersal that warrant active molecular surveillance. Sri Lanka's relative genetic isolation may offer a natural laboratory for studying the efficacy of targeted control programmes without confounding immigration of new parasite genotypes.

We call for expanded sampling-particularly from underrepresented regions such as Pakistan, Myanmar, and northeastern India-to complete the phylogeographic picture. Whole-genome approaches, combined with the baseline population structure established here, will be essential for tracking the emergence and spread of anthelmintic resistance across this globally significant livestock region.

Acknowledgement

Acknowledgements Technical StaffMs. Kanokwan Srisuk and Mr. Prasert Thongchai provided invaluable laboratory assistance with DNA extraction and PCR amplification. Their careful attention to quality control protocols ensured the reliability of the molecular dataset. CollaboratorsThis research would not have been possible without the collaboration of abattoir veterinarians and parasitologists at partner institutions across five countries. We are particularly grateful to Dr. Subhash Chandra Das (West Bengal University of Animal Sciences), Dr. Rajan Sharma (Institute of Agriculture and Animal Science, Nepal), and Dr. Amila Perera (University of Peradeniya, Sri Lanka) for coordinating specimen collection in their respective countries.

References

1. Mas-Coma S, Bargues MD, Valero MA. Fascioliasis and other plant-borne trematode zoonoses. International Journal for Parasitology. 2005;35(11-12):1255-1278.
2. Itagaki T, Tsutsumi K, Ito K, Tsutsumi Y. Taxonomic status of the Japanese triploid forms of Fasciola: comparison of mitochondrial ND1 and COI sequences with F. hepatica and F. gigantica. Journal of Parasitology. 2005;91(2):445-457.
3. Mas-Coma S, Valero MA, Bargues MD. Fasciola, lymnaeids and human fascioliasis, with a global overview on disease transmission, epidemiology, evolutionary genetics, molecular epidemiology and control. Advances in Parasitology. 2009;69:41-146.
4. Molina-Hernández V, Mulcahy G, Pérez J, MartínezMoreno A, Donnelly S, O'Neill SM, Dalton JP, Cwiklinski K. Fasciola hepatica vaccine: we may not be there yet but we're on a promising track. Veterinary Parasitology. 2015;208(1-2):101-111.
5. Cwiklinski K, O'Neill SM, Donnelly S, Dalton JP. A prospective view of animal and human Fasciolosis. Parasite Immunology. 2016;38(9):558-568.
6. Vilas R, Criscione CD, Blouin MS. A comparison between mitochondrial DNA and the ribosomal internal transcribed spacers in prospecting for cryptic species of platyhelminth parasites. Parasitology. 2005;131(6):839846.
7. Amer S, Dar Y, Ichikawa M, Fukuda Y, Tada C, Itagaki T, Nakai Y. Identification of Fasciola species isolated from Egypt based on sequence analysis of genomic (ITS1 and ITS2) and mitochondrial (NDI and COI) gene markers. Parasitology International. 2011;60(1):5-12.
8. Itagaki T, Kikawa M, Sakaguchi K, Shimo J, Terasaki K, Shibahara T, Fukuda K. Genetic characterisation of parthenogenic Fasciola sp. in Japan on the basis of the sequences of ribosomal and mitochondrial DNA. Parasitology. 2005;131(5):679-685.
9. Rajapakse RPVJ, Lawton SP, Karunathilake KJK, Itagaki T. Molecular characterisation of Fasciola from Sri Lankan livestock. Veterinary Parasitology. 2013;196(12):139-144.
10. Kelley JM, Elliott TP, Beddoe T, Anderson G, Skuce P, Spithill TW. Current threat of triclabendazole resistance in Fasciola hepatica. Trends in Parasitology. 2016;32(6):458-469.
11. Saki J, Khademvatan S, Yousefi E. Molecular identification of animal Fasciola isolates in Ahvaz, southwestern Iran. Iranian Journal of Parasitology. 2011;6(3):82-88.
12. Peng M, Ichinomiya M, Ohtori M, Ichikawa M, Shibahara T, Itagaki T. Molecular characterisation of Fasciola hepatica, Fasciola gigantica, and aspermic Fasciola sp. in China based on nuclear and mitochondrial DNA. Parasitology Research. 2009;105(3):809-815.
13. Le TH, Blair D, McManus DP. Mitochondrial genomes of parasitic flatworms. Trends in Parasitology. 2002;18(5):206-213.
14. Walker SM, Prodöhl PA, Fletcher HL, Hanna REB, Kantzoura V, Hoey EM, Trudgett A. Evidence for multiple mitochondrial lineages of Fasciola hepatica (liver fluke) within infrapopulations from cattle and sheep. Parasitology Research. 2007;101(1):117-125.
15. Alasaad S, Soriguer RC, Abu-Madi M, El Behairy A, Jowers MJ, Bañuls AL, Barroso P, Maione S, Ferrè I, Ferreras MC. A fluorescence-based polymerase chain reaction-linked single-strand conformation polymorphism (F-PCR-SSCP) assay for the identification of Fasciola spp. Parasitology Research. 2011;108(6):1513-1517.
16. Ai L, Chen MX, Alasaad S, Elsheikha HM, Li J, Li HL, Lin RQ, Zou FC, Zhu XQ, Chen JX. Genetic characterisation, species differentiation and detection of Fasciola spp. by molecular approaches. Parasites and Vectors. 2011;4:101.
17. Mas-Coma S, Bargues MD, Valero MA. Human fascioliasis infection sources, their diversity, incidence factors, analytical methods and prevention measures. Parasitology. 2018;145(13):1665-1699.
18. Cwiklinski K, Dalton JP, Dufresne PJ, La Course J, Williams DJ, Hodgkinson J, Paterson S. The Fasciola hepatica genome: gene duplication and polymorphism reveals adaptation to the host environment and the capacity for rapid evolution. Genome Biology. 2015;16:71.
19. Marcilla A, Bargues MD, Mas-Coma S. A PCR-RFLP assay for the distinction between Fasciola hepatica and Fasciola gigantica. Molecular and Cellular Probes. 2002;16(5):327-333.
20. Beesley NJ, Caminade C, Charlier J, Flynn RJ, Hodgkinson JE, Martinez-Moreno A, MartinezValladares M, Mulcahy G, O'Connor J, Sekiya M, Toet H, Williams D. Fasciola and fasciolosis in ruminants in Europe: identifying research needs. Transboundary and Emerging Diseases. 2018;65(Suppl 1):199-216.
21. Teofanova D, Kantzoura V, Walker S, Radoslavov G, Hristov P, Theodoropoulos G, Bankov I, Trudgett A. Genetic diversity of liver flukes (Fasciola hepatica) from Eastern Europe. Infection, Genetics and Evolution. 2011;11(1):109-115.