Limited genetic structure and high gene flow in Fasciola hepatica populations infecting ruminants in different geographic areas in the UK

28 The liver fluke, Fasciola hepatica, is a major parasitic threat to ruminant health and 29 productivity worldwide, with important implications for food security, animal welfare, and 30 zoonotic risk. This study developed and validated a multiplex deep amplicon sequencing 31 assay targeting the mitochondrial NADH dehydrogenase 1 (mt-ND1) and cytochrome c 32 oxidase subunit 1 (mt-COX1) loci for high-throughput genotyping of F. hepatica. DNA was 33 extracted from eggs sedimented from sheep and cattle faeces (n = 78) received from farms 34 and from adult worm pools (n = 12) isolated at abattoirs from diverse regions across the UK. 35 Following high-throughput sequencing, bioinformatics analysis was performed to 36 demultiplex Illumina sequence reads and extract amplicon sequence variants (ASVs). A total 37 of 11 ASVs were identified at each locus (mt-ND1: 264–279 bp; mt-COX1: 312–319 bp), with 38 two or three predominant ASVs per locus, along with rare variants. Network and PCA 39 analyses revealed two distinct clusters at the mt-ND1 locus: one primarily associated with 40 sheep and another shared between sheep and cattle. In contrast, mt-COX1 sequence reads 41 formed a single dominant cluster. Population analyses revealed extensive ASV sharing 42 across regions, indicating high gene flow, likely facilitated by livestock movement and 43 parasite adaptation. 44

F. hepatica, mt-ND1, mt-COX1, population genetics, ruminants, UK 45 46 47 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 1, 2026. ; https://doi.org/10.64898/2026.04.01.715781doi: bioRxiv preprint 3

48 The liver fluke genus Fasciola comprises two species: Fasciola gigantica and Fasciola hepatica. 49 F. hepatica is prevalent in temperate zones , including the UK, Europe, parts of Oceania and 50 the Americas [1], however, both species can be found in tropical and subtropical regions of 51 Asia and Africa, and their hybrids are also prevalent in Asia and Africa [1,2]. The lifecycle of 52 Fasciola is complex, with various definitive mammalian hosts, including sheep and cattle. The 53 intermediate hosts in this lifecycle are mud snails: Galba truncatula [3], previously known as 54 Lymnaea truncatula, for F. hepatica [4] and L. natalensis for F. gigantica [5]. This parasite is 55 transmitted through food plants and herbage contaminated with metacercariae, infecting 56 mainly small and large ruminants, though other mammalian species including humans can be 57 infected [6]. 58 Various factors can influence the occurrence of the parasite and resulting disease, including 59 environmental conditions such as rainfall [7], moisture levels, and temperature ( with 60 temperatures from 10°C to 25°C being optimal), as well as the geography of grazing areas 61 (e.g., topography and soil type) [8–11] and animal movement [12]. A feature of the lifecycle 62 is the clonal expansion of Fasciola spp. within its intermediate snail host, which contributes 63 to pasture contamination and to the subsequent infection of hosts by metacercariae of the 64 same genetic origin. This clonal expansion may lead to a genetic bottleneck effect in the 65 parasite, particularly when infection levels in snail populations are low [13]. 66 Understanding the population genetics of F. hepatica infection provides crucial insights to 67 aid the design of effective control strategies [14]. Over the past few decades, high levels of 68 animal movement have been reported in domestic ruminants in several European countries 69 [15,16]; hence, analysing the population genetic structure of F. hepatica can assist in 70 understanding the corresponding spread of parasites infections. Determining parasite 71 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 1, 2026. ; https://doi.org/10.64898/2026.04.01.715781doi: bioRxiv preprint 4 population genetics can further inform on transmission dynamics and infection rates, 72 thereby helping identify interventions to reduce disease burden [12]. For example, if a 73 parasite population is characterised by a single dominant amplicon sequence variant (ASV) 74 at high frequency, this suggests a single infection source with clonal expansion of the 75 parasite in the intermediate host, with low metacercariae mixing in pasture settings [12,13]. 76 On the other hand, multiple ASVs at varying frequencies in parasite populations might 77 indicate multiple infection sources on the farm and high mixing of metacercariae [13,17]. 78 Few studies of large and diverse fluke populations examine whether infection has emerged 79 recently in the host at a single time point, or whether burdens have been established 80 repeatedly at different times before spreading. Recently, we have used these methods to 81 study the multiplicity of Calicophoron daubneyi infection in the United Kingdom [17] and 82 Fasciola gigantica infection in Pakistan [12]. Our findings were consistent with multiple 83 independent emergences of C. daubneyi infection, while the identification of common 84 variants across several populations spanning a range of geographic locations highlights the 85 role of animal movements in the parasite’s spread [17]. Moreover, our findings also suggest 86 that most of the hosts were predominantly infected with the emergence of F. gigantica 87 infection, while the identification of identical variant, consistent with clonal multiplication 88 within the snails. The most common variants was identified across several populations 89 spanning a range of geographic locations, again highlighting the role of animal movements 90 in the spread of F. gigantica infections [12]. 91 Population genetics can be determined using mitochondrial DNA (mtDNA) markers [18–22], 92 as well as nuclear microsatellite loci [23]. Microsatellites are usually highly polymorphic 1-6 93 bp sequences that can be used as markers to investigate genetic diversity and genetic 94 differentiation using genomic DNA [24]. Recently, a panel of 15 highly polymorphic nuclear 95 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 1, 2026. ; https://doi.org/10.64898/2026.04.01.715781doi: bioRxiv preprint 5 microsatellite loci tested on DNA extracted from different life cycle stages, such as eggs, 96 adult worm, miracidia and metacercariae has been reported for F. hepatica [23]. High levels 97 of genetic diversity and clonal expansion of parasite and panmictic (randomly mating) 98 populations has been reported by using microsatellite markers on an abattoir-based 99 population genetics study of F. hepatica in cattle in the central England and Wales [13]. 100 There are several advantages of using microsatellite markers include their distribution over 101 eukaryotic nuclear genomes [25], and mutation information can be useful for calculating 102 Hardy-Weinberg Equilibrium (HWE) to study homozygosity and heterozygosity [23,24] 103 because genomic DNA contains information on inheritance from both parents [26]. 104 However, limitations include high mutation rates and elevated levels of polymorphism [27]. 105 Furthermore, a high number of alleles per locus in microsatellites can inflate F-statistic 106 values [28] and confound interpretation. Thus, this issue can sometimes lead to over- or 107 underestimates of genetic diversity when the most common allele occurs at either very low 108 or very high frequencies [29]. Microsatellite datasets can be also prone to genotyping 109 errors, which can bias downstream population genetic analyses [30]. 110 Mitochondrial DNA is frequently used as a marker to study population genetics because of 111 its haploid nature [31], maternal inheritance [32], high copy number [33], clock-like and 112 neutral evolution rates [34], and a lack of recombination after heteroplasmy [35]. These 113 characteristics also make mitochondrial markers valuable for investigating the genetic 114 diversity and genetic differentiation in fluke populations. Mitochondrial markers have been 115 well documented for studying population genetics in Fasciola spp. from different regions of 116 the world [12,14,22]. 117 As sequencing technology advances, deep amplicon sequencing enables the study of the 118 population genetics of Fasciola spp. in greater depth [12,21,23]. For population genetics 119 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 1, 2026. ; https://doi.org/10.64898/2026.04.01.715781doi: bioRxiv preprint 6 study,the mitochondrial marker mt-ND1 was used for F. gigantica and mt-COX1 for C. 120 daubneyi, along with deep sequencing, utilising DNA isolated from adult worms of naturally 121 infected cattle and sheep [17,36]. Deep sequencing enables the detection of dominant and 122 low-frequency variants in different parasite populations [36], which can be useful to study 123 insights into parasite transmission intensity, infection sources, and metacercarie mixing. 124 However, assays for assessing genetic diversity of F. hepatica using deep sequencing of 125 multiplexed mitochondrial markers yet not exist. Such assays would provide opportunities 126 to shed new light on the population genetics of F. hepatica from naturally infected hosts 127 across various regions of the UK. For example, the extent of overlap between sheep and 128 cattle, and the presence of geographical clustering, could provide insights into the roles of 129 pasture sharing and livestock movement in driving fluke infection and disease. 130 The present study aimed to develop and test an improved deep amplicon sequencing 131 approach to investigate the population genetics of F. hepatica infections across the UK using 132 two mitochondrial markers (mt-ND1 and mt-COX1). This multiplexed sequencing method and 133 high-throughput sequencing approach reduces experimental complexity for studying host-134 level population genetic s in F. hepatica populations using a single Illumina sequencing run. 135 This methodology enabled the examination of transmission dynamics and gene flow in both 136 adult worm and egg DNA obtained from natural F. hepatica infections in the UK. 137

138 Validation of multiplex PCR and demultiplexing of mitochondrial markers 139 The multiplex meta-barcoded PCR targeting mitochondrial positive control F. hepatica DNA 140 successfully amplified distinct bands at 311 bp for mt-ND1 and 359 bp for mt-COX1, with no 141 evidence of non-specific amplification (Supplementary Fig. 1a, b and c). 142 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 1, 2026. ; https://doi.org/10.64898/2026.04.01.715781doi: bioRxiv preprint 7 Out of a total of 90 individual samples processed, 78 (86.66%) samples generated sequence 143 reads for the mt-ND1 marker (n = 12 adult worm DNA, n = 66 egg DNA). These sequences 144 obtained from mt-ND1 loci were categorised into 40 parasite populations, comprising cattle 145 (n = 14) and sheep (n = 26) populations. For the mt-COX1 marker, 84 (93.33%) samples 146 successfully produced sequence reads (n = 12 adult worms, n = 72 egg DNA), grouped into 147 42 parasite populations, including cattle (n = 15) and sheep (n = 27). 148 Geographical distribution of mt-ND1 locus ASVs of F. hepatica 149 Eleven ASVs were identified at the F. hepatica mt-ND1 locus (accession numbers PX902280-150 PX902290), ranging from 264 bp to 270 bp, using a library of mt-ND1 reference sequences 151 downloaded from NCBI database (Supplementary Fig. 2), and their frequencies were 152 recorded across 40 fluke populations in different counties of the UK (Fig. 1a and 153 Supplementary Table 1). Across the 40 parasite populations analysed, a total of 1,403,462 154 sequence reads were extracted, of which the majority (1,167,674 reads; 83.1%) belonged to 155 a small number of predominant ASVs, including ASV1, ASV2, and ASV3. These predominant 156 ASVs reflect the dominant variants circulating in cattle and sheep in different regions. 157 ASV1 was the most abundant variant (39% of total reads), detected in 18 populations across 158 10 counties (Supplementary Table 2). It was predominant in 16 populations across seven 159 regions, exceeding 95% dominance in sheep-derived populations in Southern Scotland and 160 the West of Scotland (Supplementary Table 1). 161 ASV2 was the second most abundant variant (35.9% of total reads), found in 24 populations 162 across 13 counties (Supplementary Table 2) covering nine regions. It was predominant in 15 163 populations, reaching greater than 99% in West Midlands England and South-East England 164 sheep flocks (Supplementary Table 1). ASV3 ranked third (14.7% of total reads), occurring in 165 9 populations in eight counties (Supplementary Table 2) and in seven populations, being the 166 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 1, 2026. ; https://doi.org/10.64898/2026.04.01.715781doi: bioRxiv preprint 8 only ASV present in cattle and sheep fluke in North West England, as well as in sheep fluke 167 in South West England and the East Midlands (Supplementary Table 1). 168 In South-East and North West England, ASV3 and ASV1 predominated in fluke populations in 169 cattle and in South East England in fluke populations in sheep ASV2 predominated, with 170 minor contributions from ASV8. East of England stands out for about equal proportions of 171 ASV1 and ASV3 suggesting within-farm genetic mixing. (Supplementary Table 1). 172 In the Scottish Borders, ASV1 and ASV2 predominated in cattle fluke with ASV6 a rare 173 variant and sheep fluke populations had ASV1 (68.66%) and ASV4 (28.98%) as the major 174 contributors, and ASV2 and ASV8 as minor contributors. Southeastern Scotland cattle fluke 175 exhibited ASV3 dominance (89.19%), with minor proportions of ASV1 and ASV2. In the West 176 of Scotland, most sheep fluke populations ASV2 was most common and one population had 177 ASV10 (39.52%), a variant otherwise absent from the other populations. In Northern 178 Ireland's County Tyrone, fluke in a single sheep population predominantly showed ASV2 179 (78.16%), and other varients were ASV9 (18.15%) and ASV8 (3.69%), (Fig. 1a, Supplementary 180 Table 1). 181 Notably, some ASVs were found in some specific regions and fluke populations. For 182 example, ASV4 was found locally in two Southern Scotland sheep populations (P21S, 77.62% 183 (dominant); P25S, <0.1% (rare)). Similarly, ASV6 was rare but reached 16.1% in a single 184 Scottish Borders cattle population (P14C). ASV7 reached a high frequency (25.41%) in a 185 single Southern Scotland sheep flock (P20S). ASV8 was a widely distributed but rare variant 186 observed at ≤3% abundance but present across multiple regions. Moreover, ASV9 was 187 primarily found in Northern Ireland, while ASV10 was detected in the West of Scotland 188 sheep flock, indicating geographically isolated variants (Supplementary Table 1). 189 Network trees clustering analysis of the mt-ND-1 locus 190 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 1, 2026. ; https://doi.org/10.64898/2026.04.01.715781doi: bioRxiv preprint 9 The Neighbour-Net network tree with pie charts confirmed a highly connected parasite 191 populations dominated by two central, multi-regional ASVs, including ASV1 (the largest 192 node) and ASV2 (the second-largest node), which were present in both cattle and sheep 193 fluke across multiple regions (Fig. 2a and b). ASV1 and ASV2 were connected via ASV3 or 194 low-abundance ASV8, a genetic link between England and Scotland cattle fluke-derived 195 dominant ASVs and Scotland sheep fluke-derived dominant ASVs. Peripheral ASVs for ASV1 196 included ASV5, ASV6, ASV9, ASV10, and ASV11. Peripheral ASVs for ASV2 were ASV4 and 197 ASV7. Notably, all 11 ASVs were found in sheep fluke. In comparison, four ASVs (ASV10, 198 ASV11, ASV5, and ASV7) were not detected in cattle fluke (Fig. 2b). Median Joining Network 199 tree of mt-ND1 (Supplementary Fig. 3a) showed similar linkages among ASVs to the 200 Neighbour-Net tree (Fig. 2). 201 The mt-ND1 PCA plot based on sequence reads of 11 ASVs from 40 populations showed 202 partial clustering of F. hepatica populations by host and region, with the first two principal 203 components explaining PC1 (18.31%) and PC2 (15.18%) of the total sequence read data 204 (33.34%) (Fig. 3a). Populations of F. hepatica in sheep and cattle across the UK showed both 205 substantial overlap and some degree of regional clustering, indicating high gene flow with 206 occasional location-specific patterns. Most sheep derived populations fell into Cluster 1, 207 from 9 geographical regions including South East England, East Midlands England, South 208 West England, North West England, West Midlands England, Scottish Borders, Southern 209 Scotland, West of Scotland, and Northern Ireland. There were two regions not found in 210 cluster 1 including East of England, and Southeastern Scotland. This suggests that sheep 211 across different areas carry genetically similar parasite variants. Sheep fluke populations are 212 more evenly distributed between Cluster 1 and Cluster 2. In Cluster 2, the parasite 213 population in sheep and cattle was widely distributed in 9 regions however, samples from 214 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 1, 2026. ; https://doi.org/10.64898/2026.04.01.715781doi: bioRxiv preprint 10 West Midlands England and Northern Ireland were not represented in this cluster. This 215 indicates that the parasite variants are not geographically restricted and are prevalent in 216 both host species. 217 A few populations demonstrated regional genetic divergence. For example, one sheep fluke 218 population from Southern Scotland (P21S) clustered apart with very high PC1 and PC2 219 values, suggesting unique ASVs in the area. Similarly, Scottish Borders sheep (P9S) and 220 Southern Scotland sheep (P20S) also showed genetic separation from the main clusters. 221 The split topology tree of mt-ND1 sequences showed that one cluster is dominated by ASV1 222 and groups closely with ASVs 9 and 10. The second cluster is defined by ASV2 and ASV3, 223 indicating an evolutionary relationship between these variants. Remaining ASVs (ASV4–224 ASV11) are distributed along shorter branches in between main variants, representing low-225 frequency variants that are genetically closer to one of the two dominant clusters (Fig. 3b). 226 Genetic diversity analysis of the mt-ND-1 locus 227 The analysis of molecular variance (AMOVA) from PopArt revealed that the majority of 228 genetic variation in F. hepatica mt-ND1 populations occurred within populations (137.85%). 229 Variation among groups accounted for only 1.31% of the total, and variation among 230 populations was negative (–39.17%) (Table 1, Supplementary Fig. 3a). The negative variance 231 detected in the AMOVA results showed lack of genetic differentiation and should be 232 interpreted as zero [37,38]. All fixation indices were low and non-significant (Phi ST = –233 0.3785, P = 1.000; Phi SC = –0.3969, P = 0.997; Phi CT = 0.0131, P = 0.544), confirming the 234 absence of significant genetic differentiation between regions or populations. Similar 235 AMOVA results were confirmed by Arlequin (Supplementary File 1). These values suggest 236 high genetic connectivity and gene flow among populations across the UK, consistent with 237 the sharing of common ASVs between regions and hosts. Overall, nucleotide diversity 238 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 1, 2026. ; https://doi.org/10.64898/2026.04.01.715781doi: bioRxiv preprint 11 (π=0.00502) indicated low genetic variation in the population. Tajima's D neutrality test was 239 not significant -0.372, (p = 0.623). 240 Geographical distribution of mt-COX1 locus ASVs of F. hepatica 241 A total of 11 ASVs (312 bp to 319 bp) (accession numbers PX861700-PX861710) were 242 identified at the F. hepatica mt-COX1 locus, using reference sequences of mt-COX1 243 (Supplementary Fig. 4), and their frequencies were recorded in 42 populations in different 244 counties of the UK (Fig. 1b, and Supplementary Table 3). A total of 1.764 million sequence 245 reads were generated across the 42 parasite populations, of which 1.33 million reads 246 (75.2%) were from two predominant ASVs and 437,229 (24.8%) reads corresponded to rare 247 ASVs (Supplementary Table 3). ASV1 and ASV2 were the dominant variants circulating in 248 cattle and sheep fluke in different regions. The rare ASV3 to ASV11 were often 249 geographically restricted. 250 The most abundant variant was ASV1, contributing 45.0% of sequence reads overall and 251 detected in 41 populations across all 17 counties (Supplementary Table 4). It was the 252 predominant variant found both in cattle and sheep fluke in 22 populations, frequently 253 representing 75% of sequence reads. For example, in sheep fluke, ASV1 was observed in the 254 West Midlands England (100% of reads), South East England (>99%), North West England 255 (96%), Southern Scotland (99.46%), in the West of Scotland (>70%), Scottish Borders 256 (86.87%) and a major variant in Northern Ireland (68.23%). In cattle, ASV1 was common in 257 South West England (88.81%) and the Scottish Borders (75.2%). In North West England, 258 ASV1 and ASV2 were found in equal proportions of 50% each (Supplementary Table 3). 259 ASV2 ranked second in abundance, with overall sequence reads of 45.1% in 38 populations 260 and 16 counties. ASV2 was found to be dominant in cattle and sheep fluke across 18 261 populations from 8 regions. It was widespread in cattle fluke from the North West (>95%), 262 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 1, 2026. ; https://doi.org/10.64898/2026.04.01.715781doi: bioRxiv preprint 12 South West (>66.34%), and South East England (86.67%). ASV2 was highly prevalent in 263 sheep fluke across regions of the UK including North West England (81.07%), East of England 264 (75.77%), the West of Scotland (81.69%) and Southern Scotland (>52.06%) (Supplementary 265 Table 3). Of the other variants only ASV6 and ASV10 were common. ASV6 predominated in 266 East Midlands sheep fluke (90.95%), and ASV10 was predominant in Scottish Borders cattle 267 fluke (87.5%). The remaining ASVs were present only as rare variants in different 268 populations (Supplementary Table 3). 269 Although ASV1 and ASV2 were found to be predominant, the analysis of rare variants 270 highlighted their presence in many regions, including Northwest England, Southwest 271 England, Southern Scotland, and South Lanarkshire, ranging from 0.04% to 46.5% of 272 sequence reads. There were certain rare variants noted in sheep and cattle fluke, for 273 example, in sheep fluke ASV3 (24.2%) was in the East of England, ASV4 (34.01%) and ASV5 274 (12.64%) appeared in Scottish Borders, and ASV7 (<0.01%) in East of England. In cattle, ASV3 275 (0.04%) in South East England ASV5 (1.31%) in South West England, ASV7 (5.51%) occurred 276 in Scottish Borders, and ASV8 was noted 13.29% and <0.01% in South East England and 277 Scottish Borders, respectively. There were rare variants found in both hosts including ASV9 278 (< 4%) in South West England, East Midlands England, Scottish Borders, South Lanarkshire, 279 and West of Scotland. ASV10 (0.53% to 29.3%) appeared in Scottish Borders, South Eastern 280 Scotland, and in South West England, while ASV11 (<0.01% to 2.08%) in North West 281 England, Southern Scotland, Northern Ireland, South East England, and East Midlands 282 (Supplementary Table 3). Overall, the rare ASVs demonstrated regional and host-specific 283 patterns and may emerge in the future in both sheep and cattle populations across the UK. 284 Network trees clustering analysis of the mt-COX1 locus 285 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 1, 2026. ; https://doi.org/10.64898/2026.04.01.715781doi: bioRxiv preprint 13 The Neighbour Net analysis of the 11 mt-COX1 ASVs revealed two main nodes. ASV1 and 286 ASV2 were connected through the low-abundance ASVs 4, 6, 10, and 11, forming a genetic 287 bridge between host- and region-specific ASVs (Fig. 4a). The peripheral ASVs for ASV1 288 included ASV3, ASV7, ASV8, and ASV9. For ASV2, there was only ASV5. 289 For instance, ASV1 was most abundant in sheep fluke across Scotland and England, but was 290 also found in considerable amounts in cattle across the same regions. ASV2 was detected in 291 sheep and cattle fluke in England, as well as in sheep fluke in Scotland, with low abundance 292 in sheep fluke from Northern Ireland and in cattle fluke from Scotland. In contrast, other 293 ASVs showed strong geographic and host specificity. ASV3 and ASV6 were found in English 294 sheep, ASV4, ASV5, and ASV9 were found in Scottish sheep, ASV7 was specific to Scottish 295 cattle, ASV8 was restricted to English cattle, and ASV11 was found mainly in Northern 296 Ireland sheep (Fig. 4b). 297 Median Joining Network tree of mt-COX1 (Supplementary Fig. 3b) showed similar linkages 298 among ASVs as the Neighbour Net tree (Fig. 4a and b). 299 The PCA plot based on mt-COX1 sequence reads from 42 F. hepatica populations showed a 300 single cluster by host and region, with the first two principal components explaining PC1 301 19.07% (PC1) and 15.69%(PC2) of the total (34.76%) (Fig. 5a). All parasite populations in 302 sheep and cattle fluke across the UK show substantial overlap and indicating high gene flow. 303 The topology tree of 11 mt-COX1 ASVs showed two major ASVs (ASV1 and ASV2) at the ends 304 of the tree, and the remaining ASVs were mainly distributed along shorter branches in 305 between the main ASVs, indicating phylogenetic separation among ASVs (Fig. 5b). 306 Genetic diversity analysis of the mt-COX1 locus 307 The AMOVA results showed that the genetic variation occurred within populations (139.4%), 308 while variation among groups (0.75%) and among populations within groups (−40.2%) was 309 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 1, 2026. ; https://doi.org/10.64898/2026.04.01.715781doi: bioRxiv preprint 14 very low (Table 2, Supplementary Fig. 3b). Similar AMOVA results were found using Arlequin 310 (Supplementary File 2). Consistently, fixation indices were very low or negative (Phi ST = 311 −0.394, Phi SC = −0.405, Phi CT = 0.007), and none were statistically significant (p > 0.05). 312 Together, these results suggest that there is no significant genetic variation among groups 313 and populations. Genetic diversity was found only within populations, not across hosts or 314 geographic regions. Overall, nucleotide diversity (π = 0.00834) was low across all 315 populations. Tajima's D test was 0.375 and not significant (p = 0.343). 316

317 The present study developed a metabarcoding approach using multiplex mitochondrial 318 markers targeting the mt-ND1 and mt-COX1 loci to study the population genetics of F. 319 hepatica, using samples from natural infections collected across different areas of the UK. 320 The application of the metabarcoded multiplexed markers in deep amplicon sequencing can 321 provide an efficient tool for investigating genetic diversity patterns in multiple populations 322 of F. hepatica that can further inform transmission dynamics and infection rates. This study 323 demonstrated that F. hepatica populations with a small number of predominant ASVs were 324 circulating in both cattle and sheep with high gene flow across different regions of the UK. 325 Alongside this a few rare geographically restricted ASVs were also noted. The results 326 highlighted a pattern of widespread ASV flow across the UK, which may be driven by clonal 327 propagation of the parasite in intermediate snail hosts, livestock movement, grazing 328 practices, and parasite adaptation to the UK environment. 329 This study utilised mitochondrial markers to analyse the genetic diversity of F. hepatica 330 populations because mitochondrial DNA is useful for population genetics investigations, due 331 to its high copy number [33], maternal inheritance [32], and transmission without 332 recombination under neutral heteroplasmy conditions [35]. These properties make 333 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 1, 2026. ; https://doi.org/10.64898/2026.04.01.715781doi: bioRxiv preprint 15 mitochondrial DNA a powerful target for investigating evolutionary studies. Mitochondrial 334 markers have been used to study genetic diversity and infection dynamics in F. gigantica, 335 targeting the mt-ND1 locus where a single mt-ND1 variant of F. gigantica was predominant 336 in most of the hosts in Pakistan [12], and to investigate C. daubneyi, using the mt-COX1 337 locus. Notably, multiple variants of C. daubneyi infections, were detected across different 338 geographic regions of the UK [17]. A study from Malawi supported the suitability of both mt-339 ND1 and mt-COX1 loci for population genetics studies in F. gigantica and reported recent 340 population expansion [22]. Other studies also supported the use of mitochondrial markers 341 for genetic diversity and gene flow studies in F. hepatica [39–42]. 342 The amplification success across 90 samples, with 78 and 84 yielding sequence reads 343 for mt-ND1 and mt-COX1, respectively, highlighted the effectiveness of our method. 344 Further, the deep amplicon sequencing technique detected both dominant and rare ASVs, 345 enhancing our understanding of parasite population genetics. This ability to recover ASVs in 346 different parasite populations across geographical regions of the UK demonstrated that 347 multiplex PCR combined with next-generation sequencing can be a valuable tool for 348 studying fine-scale genetic structure. Previous studies used PCR-RFLP, conventional PCR, 349 and Sanger sequencing to amplify and analyse the mitochondrial genome regions to 350 investigate genetic variations [22,43,44]. However, these methods are low-throughput, 351 time-consuming, relatively expensive when handling medium to high numbers of samples 352 [45]. In contrast, high-throughput deep amplicon sequencing using the Illumina MiSeq 353 platform offers a convenient and cost-effective method [46] for handling a medium to high 354 number of samples, with thousands and millions of sequence reads generated per 355 population in a single run [12,17,47]. 356 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 1, 2026. ; https://doi.org/10.64898/2026.04.01.715781doi: bioRxiv preprint 16 This study identified some predominant ASVs across different regions of the UK, as well as 357 more locally restricted variants. For instance, the most common three variants of mt-ND1 358 (ASV1-3) were widely distributed across multiple regions and hosts and together accounted 359 for over 80% of infections. Similarly, two mt-COX1 variants were dominant. AMOVA results 360 showed that most genetic variation occurred within populations rather than between 361 counties or regions, which is consistent with high level of gene flow. Based on microsatellite 362 genotyping, previous research in the UK also found high genetic diversity and gene flow, as 363 well as an absence of defined population structures. This was believed to be due to the 364 clonal emergence of F. hepatica infections through the intermediate snail host [13], such 365 that a single miracidium infecting a G. truncatula mud snail can generate multiple 366 genetically identical cercariae [13]. However, metacercariae may mix on pasture, resulting in 367 a more varied genetic profile before ingestion by the definitive host [12]. 368 Findings from F. hepatica isolates from three geographical regions of China supported our 369 results, showing that genetic variation can occur within populations rather than between 370 populations [48]. In another example from the European context, in the Netherlands, 371 genetic diversity was also reported mostly within populations rather than between 372 populations [49]. In Algeria, low genetic diversity and a common origin for the parasite's 373 countrywide distribution were reported, with only two variants from the mt-COX1 gene 374 [50]. In Colombia, no genetic diversity was found among F. hepatica parasites. However, the 375 authors mentioned that this might be due to the low resolution of the molecular markers 376 used including nuclear markers (28S, β-tubulin 3, ITS1, ITS2), and mitochondrial marker (mt-377 COX1) [51]. 378 In contrast, high genetic diversity between populations was found in German dairy cattle 379 using mitochondrial (mt-ND1 and mt-COX1) and eight microsatellite markers [44], in cattle 380 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 1, 2026. ; https://doi.org/10.64898/2026.04.01.715781doi: bioRxiv preprint 17 and sheep in Spain and Peru with mt-ND1 marker [52], and in grazing cattle in Australia 381 using mt-ND1 and mt-COX1 markers [42]. Further, a high number of mt-COX1 mitochondrial 382 variants have been reported in cattle and horses in Chile [53]. A 2007 study reported that a 383 single animal in Ireland could harbour ten distinct mitochondrial variants of F. hepatica, and 384 the author related that this genetic diversity predates the last ice age [43]. A global-scale 385 analysis of NCBI data showed that both mt-ND1 and mt-COX1 locus-specific variants were 386 circulating in different parts of the world, with high Tajima's D values and a low likelihood of 387 future population growth [54]. However, from Armenia, Algeria, Brazil, Spain, and Ecuador, 388 negatively significant Tajima's D values were reported for mt-ND1 along with mt-COX1 389 showed deviation from neutrality, supporting recent population expansion [54]. This 390 contrast between the neutral, locally mixed parasite populations found in this work from UK 391 and the variable patterns observed globally highlights how local gene flow can influence in 392 future evolutionary processes in F. hepatica populations. 393 Unlike the predominant ASVs, the rare ASVs identified in this study exhibited a mostly regional 394 distribution, with only a few rare ASVs found in multiple populations across different areas, 395 for example, mt -ND1 (ASV8) and mt-COX1 (ASV9-ASV11). Th is showed the emergence of 396 localised variants , potentially linked to specific environmental factors, especially 397 environmental temperature and soil conditions, which can influence th e transmission 398 dynamics of F. hepatica [55]. For example, G. truncatula egg masses are rarely observed in 399 shaded areas [56], and this could affect opportunities for fluke transmission in different 400 environments. Fasciola eDNA was less frequently detected in dark brown or black soils, which 401 are rich in organic matter in the form of peat and have lower pH values [56]. This could impact 402 the distribution of intermediate snail hosts and the clonal expansion of rarer ASVs . Further, 403 Fasciola spp. egg hatching and development are optimal between 20 °C and 30 °C and 404 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 1, 2026. ; https://doi.org/10.64898/2026.04.01.715781doi: bioRxiv preprint 18 inhibited at temperatures below 10 °C. The time required for miracidia hatching decreased 405 with increasing temperature, and shedding of cercariae from snail hosts was most rapid at 406 around 27 °C [55]. Viability of metacercariae declined at higher temperatures but could be 407 prolonged under high humidity. Snails grow best at 25 °C, and their susceptibility to Fasciola 408 infection is also temperature dependent [55]. Climate variation and change could therefore 409 impact the geographic distribution of liver fluke, and drive adaptation to local conditions. For 410 instance, F. hepatica infections have historically been low in regions such as Southern Europe, 411 but climate shifts may increa se winter risk due to temperature and moisture fluctuations 412 falling into suitable ranges [11]. Increased fluke infections have been recorded in the EU and 413 the Northern Altiplano in South America [57]. Additionally, a study in New Zealand utilised 414 historic climate data (1972-2012) and predicted that areas with low initial risk, such as 415 Canterbury and Otago, could see a near 200% rise by 2090 [58]. Although we have currently 416 identified rare and locally restricted ASVs, these studies highlight how they could serve as 417 potential reservoirs of future genetic diversity that could become epidemiologically 418 significant under shifting environmental conditions. Further research is needed to determine 419 how genetic diversity measured using different markers to differences in biological responses 420 and environmental conditions by F. hepatica and other parasites. 421 In our study, Neighbour-Net, median joining networks and PCA showed that UK F. hepatica 422 populations are interconnected and dominated by a few ASVs. In the PCA plot, mt-ND1 423 revealed one cluster containing mainly sheep-hosted parasite populations and a second 424 cluster containing both sheep and cattle-hosted populations. In contrast, for mt-COX1, 425 parasite populations from both sheep and cattle across different regions of the UK formed a 426 single, overlapping cluster. Our study did not gather data on co-grazing of sheep and cattle. 427 However, sheep and cattle often co-graze in Northern Ireland and may be infected with the 428 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 1, 2026. ; https://doi.org/10.64898/2026.04.01.715781doi: bioRxiv preprint 19 clonal variants of parasites [43]. Prichard et al., 2005 also attributed fluke outbreaks in East 429 England to co-grazing with sheep imported from other areas and high rainfall during 430 summer [59]. This homogenisation of ASVs across populations in England, Scotland, and 431 Northern Ireland is likely due to livestock movement, a common practice in the UK 432 agricultural system and important for the spread of various diseases [60,61]. The role of 433 animal movement in parasite transmission and high gene flow has been well documented 434 [12,13,17,62]. Furthermore, greater genetic structure was observed in parasite populations 435 in sheep rather than cattle, which may be due to differences in grazing behaviour. Sheep 436 tend to feed closer to soil and waterlogged areas than cattle, potentially increasing 437 exposure [63]. In contrast, higher genetic diversity was reported in cattle fluke than in sheep 438 and goat fluke in Iran [64], but this may be because the prevalence of this disease is higher 439 in cattle than in sheep in Iran [65,66] 440 The present study has limitations. The use of only two mitochondrial markers (mt-ND1 and 441 mt-COX1) provided valuable resolution but captures only maternally inherited variation. The 442 use of nuclear markers or whole-genome data could further enhance understanding of 443 nuclear-level population genetics, and host adaptation [13,23,67]. Although 90 samples 444 from 17 counties were analysed, the sample coverage was uneven, with some regions 445 represented by only one or two populations, and no faecal samples were obtained from 446 Northern Ireland. Unequal sampling can overrepresent genetic structuring and overstate the 447 apparent dominance of a few ASVs in different regions. Finally, intermediate host snail and 448 livestock movement data were not included in this study, although both are known to 449 influence the spread of F. hepatica infection strongly. 450 Future research should aim to overcome the stated limitations by optimising nuclear 451 markers and whole-genome sequencing to complement mitochondrial markers and capture 452 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 1, 2026. ; https://doi.org/10.64898/2026.04.01.715781doi: bioRxiv preprint 20 genetic polymorphisms and adaptive behaviours (manuscript in preparation). A higher 453 number of samples across multiple years and seasons would help to track deeper dynamics 454 of ASV diversity. Linking parasite genetic data with phenotypic outcomes, such as flukicide 455 efficacy through faecal egg count reduction tests, infection intensity, and productivity losses 456 in cattle and sheep, will strengthen knowledge of fluke epidemiology and control options. 457 Assessment of genetic diversity of F. hepatica in intermediate snail hosts will provide 458 insights into bottlenecks and parasite persistence in wetlands, while combining parasite 459 genetics with livestock trade, and grazing movement data will enable causal inference about 460 how gene flow is maintained across the UK. These approaches can generate more 461 understanding of F. hepatica transmission and evolution, supporting targeted interventions 462 against fasciolosis. The methods optimised and described here provide an additional tool for 463 collecting genetic data and linking it with infection and disease outcomes, as well as 464 inferring patterns of parasite transmission and spread. 465

466 In conclusion, a multiplex mitochondrial metabarcoding approach has been developed here, 467 providing a platform for medium to large-scale population genetic studies of F. hepatica 468 infection. This study demonstrated that F. hepatica populations in the UK are largely 469 genetically interconnected and dominated by a small number of widespread variants in both 470 cattle and sheep. The findings confirm that cross-transmission of fluke between co-grazing 471 sheep and cattle is likely, although some genotypes seem to be more restricted to sheep 472 fluke. In addition, rare and region-specific variants were found at low frequencies, which 473 may contribute to the future emergence of new variants in the UK if not controlled. Our 474 findings support the idea that high gene flow can result from parasite adaptation in the UK 475 environment alongside high levels of livestock movement. These findings are important for 476 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 1, 2026. ; https://doi.org/10.64898/2026.04.01.715781doi: bioRxiv preprint 21 understanding transmission dynamics, detecting emerging variants, and informing effective 477 control strategies for F. hepatica infections to livestock farmers. 478

479 Field samples 480 A total of 90 field samples comprising 78 faecal egg samples and 12 adult worm samples 481 were selected, which identified as F. hepatica positive in our previous study [68]. Multiple 482 samples collected from the same sheep and cattle farms, veterinary practitioners, or 483 counties within close timeframes were merged into single parasite populations. 484 Additionally, adult fluke obtained from abattoirs and at post-mortem examination from the 485 same animal were treated as a single population. All populations were assigned by host 486 species, cattle (n=15) and sheep (n=27) for downstream analyses. 487 The samples were collected across 17 counties in the UK between December 2022 and May 488 2024 in collaboration with cattle and sheep farmers, as well as registered veterinary 489 practitioners, in accordance with ethical approval NASPA-2122-04 [68]. The populations 490 were further categorised into 11 regions across the UK including North West England 491 (Cheshire, Cumbria), East Midlands England (Derbyshire), West Midlands England 492 (Staffordshire), East of England (Essex), South East England (East Sussex, Kent, West Sussex), 493 South West England (Dorset, Devon, Gloucestershire, Wiltshire), Scottish Borders 494 (Peeblesshire), Southern Scotland (South Lanarkshire), Southeastern Scotland (West 495 Lothian), West of Scotland (Renfrewshire), Northern Ireland (County Tyrone). 496 Adult worm populations were obtained from the livers of infected animals at abattoirs 497 described [68]. All samples were transported to the School of Veterinary Medicine at the 498 University of Surrey and stored at –20°C for subsequent analysis. 499 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 1, 2026. ; https://doi.org/10.64898/2026.04.01.715781doi: bioRxiv preprint 22 DNA extraction 500 DNA was extracted from a pooled sample of head tissue obtained from all available adult 501 worms per host and from faecal egg samples , following the methods described [68]. Elution 502 was performed in 50 μL of molecular biology-grade water (Cytiva HyClone™), and the eluate 503 was stored at –80 °C for subsequent analyses. 504 Development of multiplex mitochondrial markers 505 PCR was conducted using Fasciola genus-specific mt-ND1 primers, resulting in a product size 506 of 311 bp [12]. In addition, mt-COX1 primers for F. hepatica were designed and tested in this 507 work (Supplementary Table 5), potentially aiming for a product length of 319-475 bp. These 508 primers were designed by aligning 435 mt-COX1 sequences from F. hepatica available on 509 Genbank and selecting a region showing variations and a conserved region after visualising 510 different sequences using Primer3 in Geneious Prime version 8.0.5. The primer pair was 20 511 bp lengths, Tm values range was 59.2-59.8 °C, GC content 50–60%, and a forward and 512 reverse Tm difference of 0.6 °C. Primers did not contain any hairpins, self-dimers, and cross-513 dimers. 514 PCR conditions were first optimised using positive control DNA in a gradient PCR with 515 annealing temperatures ranging from 52°C to 60°C for amplification of mt-COX1 516 (Supplementary Fig. 1c). PCR was performed in duplicate and twice with 2 μl of the DNA 517 template using DreamTaq Green PCR master mix (Thermo Scientific, USA) in a 25 μl reaction 518 mix and primer concentrations of 200 nM in the final volume. Final PCR conditions included 519 35 cycles of initial denaturation at 95°C for 5 minutes, followed by denaturation at 95°C for 520 1 minute, annealing of mt-ND1 primers at 50°C for 1 minute and mt-COX1 primers at 55°C 521 for 1 minute, and extension at 72°C for 1 minute, with a final extension at 72°C for 5 522 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 1, 2026. ; https://doi.org/10.64898/2026.04.01.715781doi: bioRxiv preprint 23 minutes. Water was used as a no-template control. The resulting PCR products were 523 processed for Sanger sequencing. 524 A multiplex PCR was developed using both mt-ND1 and the mt-COX1 primers, with 525 DreamTaq Green PCR master mix (Thermo Scientific, USA) with primer concentrations of 526 200 nM in 25 μl final reaction volume. The PCR was conducted for 35 cycles with conditions 527 as follows: an initial denaturation step at 95°C for 5 minutes, followed by denaturation at 528 95°C for 1 minute, annealing at 53°C for 1 minute, extension at 72°C for 1 minute and a final 529 extension at 72°C for 5 minutes. The multiplex PCR reaction was carried out in duplicate 530 using two μl of positive control F. hepatica DNA, and water as a negative control. 531 Deep amplicon sequencing of multiplexed mt-ND1 and mt-COX1 532 The metabarcoded mt-ND1 [12] and mt-COX1 (Supplementary Table 1) markers were used 533 to target the mitochondrial DNA of F. hepatica. A multiplexed first-round PCR was carried 534 out using the KAPA HiFi PCR Kit (KAPA Biosystems, South Africa) as described [12]. The 535 second-round primer sets, adaptors, barcoded PCR amplifications, magnetic bead 536 purification, and final library quantification were based on previously described methods 537 [12,17,69]. 538 Demultiplexing of mt-ND1 and mt-COX1 sequences and bioinformatics analysis 539 The Illumina MiSeq system demultiplexed the sequencing data based on sample-specific 540 barcoded indices (Supplementary Table 6), generating corresponding FASTQ files for each 541 sample (NCBI Bioproject: PRJNA1402908, accession No: SAMN54606237-SAMN54606325, 542 https://data.mendeley.com/datasets/822rxwph9t/1). The resulting FASTQ files were further 543 analysed using Mothur versions 1.41.0 and 1.48.1 [70] on the University of Surrey High-544 Performance Computing (HPC) cluster. Sequence analysis was performed following the 545 pipelines described in previous studies [12,17] with modifications as described below. 546 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 1, 2026. ; https://doi.org/10.64898/2026.04.01.715781doi: bioRxiv preprint 24 For reference database construction, mt-ND1 and mt-COX1 sequences were retrieved from 547 the NCBI database for F. hepatica and F. gigantica. A total of 363 mt-ND1 sequences and 548 462 mt-COX1 sequences were downloaded for F. hepatica, which were collapsed into 79 549 and 97 unique reference sequences, respectively. For F. gigantica, 351 mt-ND1 and 337 mt-550 COX1 sequences were downloaded, resulting in 117 and 108 unique collapsing sequences, 551 respectively. These reference sequence libraries for the mt-ND1 and mt-COX1 genes were 552 used for sequence demultiplexing and alignment. The Mothur pipeline joined paired-end 553 reads, filtered out ambiguous or low-quality sequences, and removed excessively long or 554 short sequences. Sequences were aligned against the reference library, unique sequences 555 were pre-clustered and abundant reads were then grouped. 556 ASVs were obtained from the filtered dataset using minimum read thresholds assigned to 557 eliminate noise and sequencing artefacts using the command "split.abund". A cutoff value 558 of 7,000 reads per ASV was applied to the mt-ND1 dataset, and 2,000 reads per ASV were 559 used for mt-COX1. These thresholds were determined empirically following visual inspection 560 of the output count table files. The aligned unique sequences were split into a high- and 561 low-abundance sequence read count table and FASTA files based on the defined cutoff 562 value. Sequences with low abundances and below threshold were separated and discarded, 563 while high-abundance sequences were selected for downstream analyses. Following ASV 564 extraction for all samples, downstream processing, including sample wise sequence cleaning 565 and sorting of FASTA files by specific ASV names, was performed in R using the phytools 566 [71], microseq [72], biostrings [73], dplyr [74], ggh4x [75] and tidyverse [76] packages. A 567 FASTA file and group file containing corresponding abundance data generated by Mothur 568 were used for downstream analysis. The dataset was standardised to prevent mismatches. 569 The data was used to generate population-specific sequence fasta files, where sequences 570 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 1, 2026. ; https://doi.org/10.64898/2026.04.01.715781doi: bioRxiv preprint 25 were replicated according to their abundance values. Biostrings-based function was used to 571 each FASTA file to improve data quality. This function removed ambiguous nucleotides (e.g., 572 'N') and non-ATCG characters and filtered out sequences shorter than 100 base pairs. The 573 cleaned sequences were saved into new FASTA files, and ASVs were named in each sample 574 based on the descending order of sequence read abundance. ASV sequences were 575 visualised using Geneious version 8.0.0 (https://www.geneious.com). Finally, sample files 576 are grouped into populations using an R script, and unique sequences were extracted for 577 each population for downstream analysis. All R scripts and the reference sequences library 578 used for this process are available at the Mendeley data repository 579 (https://data.mendeley.com/datasets/822rxwph9t/1). 580 Phylogenetic, network and Split Tree analysis 581 Phylogenetic trees were generated from unique reference sequences of mt -ND1 and mt-582 COX1 for F. hepatica and F. gigantica, downloaded from NCBI GenBank. The sequences were 583 aligned using MUSCLE in Geneious v8.0.5. Further, phylogenetic trees were constructed using 584 the Neighbour-Joining method [77]. The evolutionary distances were computed using the 585 Maximum Composite Likelihood method [78] in MEGA11 [79] with a bootstrap value of 2000 586 [80]. 587 Split trees were generated using SplitTrees4 CE 6.0.0 [81], employing the HKY85 Distance 588 [82] and Neighbor Net method [83,84]. The most appropriate nucleotide substitution model 589 for HKY85 Distance was identified using jModelTest 12.2.0 [85]. Split topology tree was 590 generated using UPGMA method [86], and with the 1000 Bootstraps [80]. Moreover, the 591 Median Joining Network tree, nucleotide diversity, Tajima D, and AMOVA analysis were 592 performed using popart-1.7 [87,88]. AMOVA analysis was further confirmed using Arlequin 593 [89]. 594 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 1, 2026. ; https://doi.org/10.64898/2026.04.01.715781doi: bioRxiv preprint 26 Data analysis 595 Pie charts, frequency distributions, and analyses of predominant and rare ASV patterns, as 596 well as the proportional distribution of ASVs within each county and across the population, 597 were generated to characterise F. hepatica populations in sheep and cattle in R using 598 packages readxl [90], ggplot2 [91], dplyr [74], and tidyr [92]. Location data points from 599 confirmed positive collection sites were plotted on the UK map. Geographic coordinates 600 (longitude and latitude) for each site were taken from Google Maps (Supplementary Table 601 7). Mapping was performed using spatial data sourced from the UK Data Service, including 602 Census Support Digitised Boundary Data (1840–present) and Postcode Directories (1980–603 present), which allowed for the accurate visualisation of ASV distribution patterns across the 604 UK. All data analysis and visualisations were performed in R version 4.3.3 (https://cran.r-605 project.org/). 606

607 1. Mas-Coma, S., Valero, M. A. & Bargues, M. D. Human and animal fascioliasis: origins and 608 worldwide evolving scenario. Clin. Microbiol. Rev. 35, e0008819 (2022). 609 2. Mas-Coma, S., Valero, M. A. & Bargues, M. D. Fascioliasis. Adv. Exp. Med. Biol. 1154, 71–103 610 (2019). 611 3. Novobilský, A. et al. Transmission patterns of Fasciola hepatica to ruminants in Sweden. Vet. 612 Parasitol. 203, 276–286 (2014). 613 4. Dreyfuss, G. & Rondelaud, D. Fasciola hepatica: a study of the shedding of cercariae from 614 Lymnaea truncatula raised under constant conditions of temperature and photoperiod. Parasite 615 Paris Fr. 1, 401–404 (1994). 616 5. Pfukenyi, D. M. & Mukaratirwa, S. A retrospective study of the prevalence and seasonal 617 variation of Fasciola gigantica in cattle slaughtered in the major abattoirs of Zimbabwe between 618 1990 and 1999. Onderstepoort J. Vet. Res. 71, 181–187 (2004). 619 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 1, 2026. ; https://doi.org/10.64898/2026.04.01.715781doi: bioRxiv preprint 27 6. Hotez, P. J. et al. The Global Burden of Disease Study 2010: Interpretation and Implications for 620 the Neglected Tropical Diseases. PLoS Negl. Trop. Dis. 8, e2865 (2014). 621 7. Sweeney, J. et al. Climate change refining the impacts for Ireland: Strive Report (2001-CD-C3-622 M1) ISBN: 978-1-84095-297-1. (2008). 623 8. Caminade, C., van Dijk, J., Baylis, M. & Williams, D. Modelling recent and future climatic 624 suitability for fasciolosis in Europe. Geospatial Health 9, 301–308 (2015). 625 9. Rapsch, C. et al. An interactive map to assess the potential spread of Lymnaea truncatula and 626 the free-living stages of Fasciola hepatica in Switzerland. Vet. Parasitol. 154, 242–249 (2008). 627 10. Howell, A. K. & Williams, D. J. L. The epidemiology and control of liver flukes in cattle and sheep. 628 Vet. Clin. North Am. Food Anim. Pract. 36, 109–123 (2020). 629 11. Kenyon, F., Sargison, N. D., Skuce, P. J. & Jackson, F. Sheep helminth parasitic disease in south 630 eastern Scotland arising as a possible consequence of climate change. Vet. Parasitol. 163, 293–631 297 (2009). 632 12. Rehman, Z. U. et al. Genetic diversity and multiplicity of infection in Fasciola gigantica isolates of 633 Pakistani livestock. Parasitol. Int. 76, 102071 (2020). 634 13. Beesley, N. J., Williams, D. J. L., Paterson, S. & Hodgkinson, J. Fasciola hepatica demonstrates 635 high levels of genetic diversity, a lack of population structure and high gene flow: possible 636 implications for drug resistance. Int. J. Parasitol. 47, 11–20 (2017). 637 14. Hodgkinson, J., Cwiklinski, K., Beesley, N. J., Paterson, S. & Williams, D. J. L. Identification of 638 putative markers of triclabendazole resistance by a genome-wide analysis of genetically 639 recombinant Fasciola hepatica. Parasitol. 140, 1523–1533 (2013). 640 15. Kelley, J. M. et al. Current Threat of Triclabendazole Resistance in Fasciola hepatica. Trends 641 Parasitol. 32, 458–469 (2016). 642 16. Vilas, R., Vázquez-Prieto, S. & Paniagua, E. Contrasting patterns of population genetic structure 643 of Fasciola hepatica from cattle and sheep: implications for the evolution of anthelmintic 644 resistance. Infect. Genet. Evol. J. Mol. Epidemiol. Evol. Genet. Infect. Dis. 12, 45–52 (2012). 645 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 1, 2026. ; https://doi.org/10.64898/2026.04.01.715781doi: bioRxiv preprint 28 17. Sargison, N. D., Shahzad, K., Mazeri, S. & Chaudhry, U. A high throughput deep amplicon 646 sequencing method to show the emergence and spread of Calicophoron daubneyi rumen fluke 647 infection in United Kingdom cattle herds. Vet. Parasitol. 268, 9–15 (2019). 648 18. Itagaki, T., Kikawa, M., Terasaki, K., Shibahara, T. & Fukuda, K. Molecular characterization of 649 parthenogenic Fasciola sp. in Korea on the basis of DNA sequences of ribosomal ITS1 and 650 mitochondrial NDI gene. J. Vet. Med. Sci. 67, 1115–1118 (2005). 651 19. Itagaki, T. et al. Genetic characterization of parthenogenic Fasciola sp. in Japan on the basis of 652 the sequences of ribosomal and mitochondrial DNA. Parasitol. 131, 679–685 (2005). 653 20. Chaudhry, U. et al. Molecular evidence shows that the liver fluke Fasciola gigantica is the 654 predominant Fasciola species in ruminants from Pakistan. J. Helminthol. 90, 206–213 (2016). 655 21. Rehman, Z. U. et al. High-throughput sequencing of Fasciola spp. shows co-infection and 656 intermediate forms in Balochistan, but only Fasciola gigantica in the Punjab province of 657 Pakistan. Infect. Genet. Evol. 94, 105012 (2021). 658 22. Mogha, L. et al. Genetic diversity and population structure of Fasciola gigantica isolated from 659 cattle in Malawi. Vet. Res. Commun. 49, 157 (2025). 660 23. Cwiklinski, K. et al. Characterisation of a novel panel of polymorphic microsatellite loci for the 661 liver fluke, Fasciola hepatica, using a next generation sequencing approach. Infect. Genet. Evol. 662 32, 298–304 (2015). 663 24. Ellegren, H. Microsatellites: simple sequences with complex evolution. Nat. Rev. Genet. 5, 435–664 445 (2004). 665 25. Stallings, R. L. et al. Evolution and distribution of (GT)n repetitive sequences in mammalian 666 genomes. Genomics 10, 807–815 (1991). 667 26. Johnson, P. C. D. et al. Abundant variation in microsatellites of the parasitic nematode 668 Trichostrongylus tenuis and linkage to a tandem repeat. Mol. Biochem. Parasitol. 148, 210–218 669 (2006). 670 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 1, 2026. ; https://doi.org/10.64898/2026.04.01.715781doi: bioRxiv preprint 29 27. Hodel, R. G. J. et al. The report of my death was an exaggeration: A review for researchers using 671 microsatellites in the 21st century1. Appl. Plant Sci. 4, apps.1600025 (2016). 672 28. Whitlock, M. C. G’ST and D do not replace FST. Mol. Ecol. 20, 1083–1091 (2011). 673 29. Jakobsson, M., Edge, M. D. & Rosenberg, N. A. The relationship between F(ST) and the frequency 674 of the most frequent allele. Genetics 193, 515–528 (2013). 675 30. Hoffman, J. I. & Amos, W. Microsatellite genotyping errors: detection approaches, common 676 sources and consequences for paternal exclusion. Mol. Ecol. 14, 599–612 (2005). 677 31. Ballard, J. W. O. & Whitlock, M. C. The incomplete natural history of mitochondria. Mol. Ecol. 13, 678 729–744 (2004). 679 32. Schwartz, M. & Vissing, J. Paternal inheritance of mitochondrial DNA. N. Engl. J. Med. 347, 576–680 580 (2002). 681 33. Bogenhagen, D. & Clayton, D. A. The number of mitochondrial deoxyribonucleic acid genomes in 682 mouse L and human HeLa cells: quantitative isolation of mitochondrial deoxyribonucleaic acid. J. 683 Biol. Chem. 249, 7991–7995 (1974). 684 34. Galtier, N., Nabholz, B., Glémin, S. & Hurst, G. D. D. Mitochondrial DNA as a marker of molecular 685 diversity: a reappraisal. Mol. Ecol. 18, 4541–4550 (2009). 686 35. Hagström, E., Freyer, C., Battersby, B. J., Stewart, J. B. & Larsson, N.-G. No recombination of 687 mtDNA after heteroplasmy for 50 generations in the mouse maternal germline. Nucleic Acids 688 Res. 42, 1111–1116 (2014). 689 36. Rehman, Z. U. et al. Genetic diversity and multiplicity of infection in Fasciola gigantica isolates of 690 Pakistani livestock. Parasitol. Int. 76, 102071 (2020). 691 37. Kiær, L. P. et al. Spontaneous gene flow and population structure in wild and cultivated chicory, 692 Cichorium intybus L. Genet. Resour. Crop Evol. 56, 405–419 (2009). 693 38. Moudgil, A. D. & Nehra, A. K. Mitochondrial genetic markers based phylogenetic analyses of 694 Hyalomma dromedarii Koch, 1844 (Acari: Ixodidae). J. Genet. Eng. Biotechnol. 23, 100460 695 (2025). 696 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 1, 2026. ; https://doi.org/10.64898/2026.04.01.715781doi: bioRxiv preprint 30 39. Amor, N., Farjallah, S., Merella, P., Alagaili, A. N. & Mohammed, O. B. Multilocus approach 697 reveals discordant molecular markers and corridors for gene flow between North African 698 populations of Fasciola hepatica. Vet. Parasitol. 278, 109035 (2020). 699 40. Shoriki, T. et al. Molecular phylogenetic identification of Fasciola flukes in Nepal. Parasitol. Int. 700 63, 758–762 (2014). 701 41. Salim, B. et al. Comprehensive mitochondrial genomics of Fasciola gigantica from Sudan: 702 insights into genetic diversity, evolutionary dynamics, and host adaptation. Front. Vet. Sci. 12, 703 (2025). 704 42. Elliott, T. et al. Evidence for high genetic diversity of NAD1 and COX1 mitochondrial haplotypes 705 among triclabendazole resistant and susceptible populations and field isolates of Fasciola 706 hepatica (liver fluke) in Australia. Vet. Parasitol. 200, 90–96 (2014). 707 43. Walker, S. M. et al. Evidence for multiple mitochondrial lineages of Fasciola hepatica (liver fluke) 708 within infrapopulations from cattle and sheep. Parasitol. Res. 101, 117–125 (2007). 709 44. Hecker, A. S., Raulf, M.-K., König, S., May, K. & Strube, C. Population genetic analysis of the liver 710 fluke Fasciola hepatica in German dairy cattle reveals high genetic diversity and associations 711 with fluke size. Parasit. Vectors 18, 51 (2025). 712 45. Calvani, N. E. D. & Šlapeta, J. Fasciola species introgression: just a fluke or something more? 713 Trends Parasitol. 37, 25–34 (2021). 714 46. Cho, Y. et al. Comparative analysis of Sanger and next generation sequencing methods for 16S 715 rDNA analysis of post-mortem specimens. Aust. J. Forensic Sci. 51, 426–445 (2019). 716 47. Zahid, O., Chaudhry, U. & Sargison, N. D. Development and field application of metabarcoding-717 adapted mt-ND4 markers shows substantial gene flow and varying local pressures on 718 Haemonchus contortus and Teladorsagia circumcincta populations in the UK. PLOS ONE 20, 719 e0327254 (2025). 720 48. Xifeng, W. et al. Molecular characteristics and genetic diversity of Fasciola hepatica from sheep 721 in Xinjiang, China. J. Vet. Res. 66, 199–207 (2022). 722 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 1, 2026. ; https://doi.org/10.64898/2026.04.01.715781doi: bioRxiv preprint 31 49. Walker, S. M. et al. Population dynamics of the liver fluke, Fasciola hepatica: the effect of time 723 and spatial separation on the genetic diversity of fluke populations in the Netherlands. Parasitol. 724 138, 215–223 (2011). 725 50. Chaouadi, M. et al. New insights into the genetic variability of Fasciola hepatica (trematoda) in 726 Algeria and relationships with other geographic regions revealed by mitochondrial DNA. 727 Helminthologia 59, 152–164 (2022). 728 51. Garcia-Corredor, D. et al. Molecular characterization of Fasciola hepatica in endemic regions of 729 Colombia. Front. Vet. Sci. 10, (2023). 730 52. Thang, T. N. et al. Genetic diversity of Fasciola hepatica in Spain and Peru. Parasitol. Int. 76, 731 102100 (2020). 732 53. Cabrera, G. et al. Molecular characterization of Fasciola hepatica obtained from cattle and horse 733 in Central Chile. Vet. Parasitol. Reg. Stud. Rep. 56, 101130 (2024). 734 54. Alvi, M. A. et al. Genetic variation and population structure of Fasciola hepatica: an in silico 735 analysis. Parasitol. Res. 122, 2155–2173 (2023). 736 55. Modabbernia, G., Meshgi, B. & Kinsley, A. C. Climatic variations and Fasciola: a review of impacts 737 across the parasite life cycle. Parasitol. Res. 123, 1–14 (2024). 738 56. Jones, R. A. et al. Identification of factors associated with Fasciola hepatica infection risk areas 739 on pastures via an environmental DNA survey of Galba truncatula distribution using droplet 740 digital and quantitative real-time PCR assays. Environ. DNA 6, e371 (2024). 741 57. Mas-Coma, S., Valero, M. A. & Bargues, M. D. Climate change effects on trematodiases, with 742 emphasis on zoonotic fascioliasis and schistosomiasis. Vet. Parasitol. 163, 264–280 (2009). 743 58. Haydock, L. A. J., Pomroy, W. E., Stevenson, M. A. & Lawrence, K. E. A growing degree-day model 744 for determination of Fasciola hepatica infection risk in New Zealand with future predictions 745 using climate change models. Vet. Parasitol. 228, 52–59 (2016). 746 59. Pritchard, G. C., Forbes, A. B., Williams, D. J. L., Salimi-Bejestani, M. R. & Daniel, R. G. Emergence 747 of fasciolosis in cattle in East Anglia. Vet. Rec. 157, 578–582 (2005). 748 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 1, 2026. ; https://doi.org/10.64898/2026.04.01.715781doi: bioRxiv preprint 32 60. Green, D. m, Kiss, I. z & Kao, R. r. Modelling the initial spread of foot-and-mouth disease through 749 animal movements. Proc. R. Soc. B Biol. Sci. 273, 2729–2735 (2006). 750 61. Volkova, V. V., Howey, R., Savill, N. J. & Woolhouse, M. E. J. Sheep movement networks and the 751 transmission of infectious diseases. PLoS ONE 5, e11185 (2010). 752 62. Semyenova, S. K. et al. Genetic differentiation in eastern European and western Asian 753 populations of the liver fluke, Fasciola hepatica, as revealed by mitochondrial NAD1 and COX1 754 genes. J. Parasitol. 92, 525–530 (2006). 755 63. Department for Environment, Food & Rural Affairs, UK. Graze with livestock to maintain and 756 improve habitats – Farming. https://defrafarming.blog.gov.uk/graze-with-livestock-to-maintain-757 and-improve-habitats/ (2025). 758 64. Javanmard, E. et al. Multigene typing and phylogenetic analysis of Fasciola from endemic foci in 759 Iran. Infect. Genet. Evol. 80, 104202 (2020). 760 65. Ashrafi, K. The status of human and animal fascioliasis in Iran: A narrative review article. Iran. J. 761 Parasitol. 10, 306–328 (2015). 762 66. Moghaddam, A. S. et al. Human and animal fascioliasis in Mazandaran province, northern Iran. 763 Parasitol. Res. 94, 61–69 (2004). 764 67. Cwiklinski, K. et al. The Fasciola hepatica genome: gene duplication and polymorphism reveals 765 adaptation to the host environment and the capacity for rapid evolution. Genome Biol. 16, 71 766 (2015). 767 68. Abbas, M. et al. Development of a qPCR assay for Fasciola spp. identification and deep amplicon 768 sequencing method for differentiation of fluke species in UK livestock. PLOS Negl Trop Dis 20, 769 e0014006 (2026). 770 69. Yasein, G. et al. A novel metabarcoded deep amplicon sequencing tool for disease surveillance 771 and determining the species composition of Trypanosoma in cattle and other farm animals. Acta 772 Trop. 230, 106416 (2022). 773 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 1, 2026. ; https://doi.org/10.64898/2026.04.01.715781doi: bioRxiv preprint 33 70. Schloss, P. et al. Introducing mothur: open-source, platform-independent, community-774 supported software for describing and comparing microbial communities. Appl. Environ. 775 Microbiol. 75, (2009). 776 71. Revell, L. J. phytools 2.0: an updated R ecosystem for phylogenetic comparative methods (and 777 other things). PeerJ 12, e16505 (2024). 778 72. Snipen, L. & Liland, K. H. microseq: basic biological sequence handling. 779 https://doi.org/https://cran.r-project.org/web/packages/microseq/index.html (2025). 780 73. Pagès, H., Aboyoun, P., Gentleman, R. & DebRoy, S. Biostrings: Efficient manipulation of 781 biological strings. R Package Version 2, 10–18129 (2019). 782 74. Wickham, H. et al. dplyr: A grammar of data manipulation. https://doi.org/https://cran.r-783 project.org/web/packages/dplyr/index.html (2023). 784 75. Brand, T. van den. ggh4x: Hacks for ‘ggplot2’. https://doi.org/https://cran.r-785 project.org/web/packages/ggh4x/index.html (2025). 786 76. Wickham, H. et al. Welcome to the tidyverse. J. Open Source Softw. 4, 1686 (2019). 787 77. Saitou, N. & Nei, M. The neighbor-joining method: a new method for reconstructing 788 phylogenetic trees. Mol. Biol. Evol. 4, 406–425 (1987). 789 78. Tamura, K., Nei, M. & Kumar, S. Prospects for inferring very large phylogenies by using the 790 neighbor-joining method. Proc. Natl. Acad. Sci. 101, 11030–11035 (2004). 791 79. Tamura, K., Stecher, G. & Kumar, S. MEGA11: Molecular evolutionary genetics analysis version 792 11. Mol. Biol. Evol. 38, 3022–3027 (2021). 793 80. Felsenstein, J. Confidence limits on phylogenies: an approach using the bootstrap. Evolution 39, 794 783–791 (1985). 795 81. Huson, D. H. & Bryant, D. Application of phylogenetic networks in evolutionary studies. Mol. 796 Biol. Evol. 23, 254–267 (2006). 797 82. Hasegawa, M., Kishino, H. & Yano, T. Dating of the human-ape splitting by a molecular clock of 798 mitochondrial DNA. J. Mol. Evol. 22, 160–174 (1985). 799 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 1, 2026. ; https://doi.org/10.64898/2026.04.01.715781doi: bioRxiv preprint 34 83. Bryant, D. & Moulton, V. Neighbor-net: An agglomerative method for the construction of 800 phylogenetic networks. Mol. Biol. Evol. 21, 255–265 (2004). 801 84. Bryant, D. & Huson, D. H. NeighborNet: improved algorithms and implementation. Front. 802 Bioinforma. 3, (2023). 803 85. Posada, D. jModelTest: phylogenetic model averaging. Mol. Biol. Evol. 25, 1253–1256 (2008). 804 86. Sokal, R. R. & Michener, C. D. A statistical method for evaluating systematic relationships. 805 (1958). 806 87. Bandelt, H. J., Forster, P. & Röhl, A. Median-joining networks for inferring intraspecific 807 phylogenies. Mol. Biol. Evol. 16, 37–48 (1999). 808 88. Leigh, J. W. & Bryant, D. popart: full-feature software for haplotype network construction. 809

Ecol. Evol. 6, 1110–1116 (2015). 810 89. Excoffier, L., Laval, G. & Schneider, S. Arlequin (version 3.0): An integrated software package for 811 population genetics data analysis. Evol. Bioinforma. Online 1, 47–50 (2007). 812 90. Wickham, H. et al. readxl: Read Excel Files. https://doi.org/https://cran.r-813 project.org/web/packages/readxl/index.html (2025). 814 91. CRAN: ggplot2 citation info. https://cran.r-project.org/web/packages/ggplot2/citation.html. 815 92. Wickham, H. et al. tidyr: Tidy Messy Data. https://doi.org/https://cran.r-816 project.org/web/packages/tidyr/index.html (2025). 817

818 Part of this work was carried out using computational HPC facilities and support provided by 819 the Research Computing Services team within IT Services at the University of Surrey, 820 specifically the Eureka2 HPC cluster 821 (https://docs.pages.surrey.ac.uk/research_computing/hpc/clusters/eureka2.html). 822 This research was funded by the UK Research and Innovation (UKRI), Biotechnology and 823 Biological Sciences Research Council (BBSRC) through the FoodBioSystems Doctoral Training 824 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 1, 2026. ; https://doi.org/10.64898/2026.04.01.715781doi: bioRxiv preprint 35 Programme (BB/T008776/1) and by the Sir Halley Stewart Trust (3153). For Open Access, the 825 authors have applied a Creative Commons Attribution (CC BY) public copyright license to any 826 Author Accepted Manuscript version arising from this submission. 827 We sincerely acknowledge all farmers and registered veterinary practitioners in the UK and 828 especially Dr. Iñaki Deza -Cruz (The Royal (Dick) School of Veterinary Studies and The Roslin 829 Institute, The University of Edinburgh, Easter Bush Veterinary Centre, Midlothian, EH25 9RG) 830 for reading the manuscript and sample collection. Dr. Sai Fingerhood (Department of 831 Veterinary Pathology, University of Nottingham, UK), and Dr. Mark W. Robinson (School of 832 Biological Sciences, Queen's University Belfast, UK) for their valuable assistance in sample 833 collection. 834 Contributions 835 Muhammad Abbas: conceptualisation, investigation, methodology, bioinformatics, 836 validation, visualisation, data curation and analysis, writing original draft, review and editing; 837 Kezia Kozel: methodology, writing review and editing; Nick Selemetas: writing review and 838 editing, supervision; Olukayode Daramola: writing review and editing , supervision; Eric R 839 Morgan: conceptualisation, funding acquisition , supervision, writing review and editing; 840 Umer Chaudhry: conceptualisation , writing review and editing, supervision; Martha Betson: 841 conceptualisation, writing review and editing, supervision, funding acquisition, project 842 administration. 843 Ethical statement 844 Non-invasive collection of faecal samples was approved by the NASPA (Non-Animal 845 Scientific Procedures Act) sub-committee of AWERB, University of Surrey, UK, under the 846

NASPA-2122-04 for the project "Developing Novel Rapid Diagnostics for 847 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 1, 2026. ; https://doi.org/10.64898/2026.04.01.715781doi: bioRxiv preprint 36 Neglected Parasitic Diseases." Adult F. hepatica were collected at licensed slaughterhouses 848 and through post-mortem examination. Completion of a University of Surrey SAGE-AR (ID 849 638929-638920-101535552) indicated that no formal ethical approval was required for 850 adult fluke sampling. 851 Supplementary information 852 Supplementary Fig. 1 to 4 853 Supplementary Files. 1 to 2 854 Supplementary Tables 1 to 7 855 Rights and permissions 856 All sequencing data reported in the paper are available under NCBI BioProject ID 857 PRJNA1402908 and accession numbers : SAMN54606237 -SAMN54606325, PX861700-858 PX861710 and PX902280-PX902290. 859 In addition, sequence data, R script, and codes are available at the Mendeley da tabase 860 https://data.mendeley.com/datasets/822rxwph9t/1 861 All other data are reported in the paper and associated supplementary material. 862 Funding 863 Muhammad Abbas received funding from the UK Research and Innovation (UKRI), 864 Biotechnology and Biological Sciences Research Council (BBSRC) through the FoodBioSystems 865 Doctoral Training Programme for project ID FBS2022 titled "New tools for sustainable control 866 of liver fluke in ruminants" Grant Ref: BB/T008776/1. Further, this research was funded by 867 the Sir Halley Stewart Trust under the project "Rapid Diagnostics for Neglected Parasites. 868 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 1, 2026. ; https://doi.org/10.64898/2026.04.01.715781doi: bioRxiv preprint 37 Competing Interest 869 The authors declare that no financial interests or personal relationships could have influenced 870 the work reported in this paper.871 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 1, 2026. ; https://doi.org/10.64898/2026.04.01.715781doi: bioRxiv preprint 38 Fig. 1. The relative frequencies of mt-ND1 and mt-COX1 ASVs in F. hepatica from sheep (S) and cattle (C) across 17 counties in the UK. (a) mt-ND1 ASVs and (b) mt-COX1 ASVs in 40 and 42 F. hepatica populations, respectively. Each pie chart represents a distinct population, originating either from adult worms (indicated with an asterisk *) or eggs purified from faeces. Individual ASVs are distinguished by different colours within the charts . The size and composition of each pie chart reflect the proportional distribution of ASVs within the respective population. Furthermore, each population is mapped to its geographical collection site, providing a visual representation of ASV diversity and distribution across the UK. .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 1, 2026. ; https://doi.org/10.64898/2026.04.01.715781doi: bioRxiv preprint 39 Fig. 2. Network tree and clustering of all ASVs based on region and host (a) mt-ND1 Neighbour Net tree in Split tree. Each pie chart shows regions represented by different colours, with representative ASVs displayed. The pie chart represents the ASV distribution, and its frequency in all populations found in the region. The branch lengths were calculated using the HKY85 Distance method, as determined to be best by jModeltest 2.1.10. (b) Each pie chart presents the countries of the UK, represented by different colours: England cattle (red), England sheep (light red), Northern Ireland sheep (yellow), and Scotland cattle (blue), and Scotland sheep (light blue), where representative ASVs were recorded. The pie chart shows the ASV distribution and its read frequency across all populations in the countries. .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 1, 2026. ; https://doi.org/10.64898/2026.04.01.715781doi: bioRxiv preprint 40 Fig. 3. Principal Component Analysis (PCA) plot and topology tree of all ASVs based on region and host. (a) PCA of F. hepatica populations based on 11 mt-ND1 ASV sequence abundance. Each point represents a population, with symbols representing the different regions and colour the host species. The axes represent the first two principal components (PC1 and PC2), which explain 18.31% and 15.18% of the variance, respectively. (b) Split topology tree of mt-ND1 with the UPGMA method. The pie chart circles in the tree represent the frequency of ASVs sequence reads in all populations. .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 1, 2026. ; https://doi.org/10.64898/2026.04.01.715781doi: bioRxiv preprint 41 Fig. 4. Network tree and clustering of all ASVs based on region and host (a) mt-COX1 Neighbor Net tree in Split tree. Each pie chart shows regions represented by different colours, with representative ASVs recorded. The pie chart circle represents the ASV distribution, and its sequence reads frequency in all populations found in the region. The branch lengths were calculated using the HKY85 Distance method, as determined to be best by jModeltest 2.1.10. (b) Each pie chart shows the countries of the UK, represented by different colours: England cattle (red), England sheep (light red), Northern Ireland sheep (yellow), and Scotland cattle (blue) and Scotland sheep (light blue), with representative ASVs recorded. The pie chart shows the ASV distribution and read frequency across all populations in the countries. .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 1, 2026. ; https://doi.org/10.64898/2026.04.01.715781doi: bioRxiv preprint 42 Fig. 5. Principal Component Analysis (PCA) plot and topology tree of all mt-COX1 ASVs (a) PCA of F. hepatica populations based on 11 mt-COX1 ASV sequence abundance. Each point represents a population, with symbols representing the different regions and colour the host species. The axes represent the first two principal components (PC1 and PC2), which explain 19.07% and 15.69% of the variance, respectively. (b) Split topology tree of mt-COX1 using the UPGMA method. The pie charts represent the frequency of ASVs' sequence reads in all populations. .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 1, 2026. ; https://doi.org/10.64898/2026.04.01.715781doi: bioRxiv preprint 43 Table 1. AMOVA for 11 ASVs of mt-ND1 Source of variation df Sum of squares Variance (σ²) % variation Among groups1 16 10.386 0.017 1.3136% Among populations within groups2 23 14.502 –0.508 –39.17% Within populations 54 96.633 1.790 137.853% Total 93 121.521 1.298 100% Note: The populations are grouped as follows: 1 all populations observed in a county were categorised as a group; 2 populations found within the group. The negative variance component observed in the AMOVA indicates the absence of genetic structure and should be considered as zero. Table 2. AMOVA for 11 AVS of mt-COX1 Source of Variation df Sum of Squares Variance (σ²) % Variation Among groups1 16 29.009 0.040 0.74717 % Among populations within groups2 25 50.794 -2.153 -40.16418 % Within populations 68 508.133 7.473 139.41701 % Total 109 587.936 5.360 100 % Note: The populations are grouped as follows: 1 all populations observed in a county were categorised as a group; 2 populations found within the group. The negative variance component observed in the AMOVA indicates the absence of genetic structure and should be considered as zero. .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 1, 2026. ; https://doi.org/10.64898/2026.04.01.715781doi: bioRxiv preprint .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 1, 2026. ; https://doi.org/10.64898/2026.04.01.715781doi: bioRxiv preprint .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 1, 2026. ; https://doi.org/10.64898/2026.04.01.715781doi: bioRxiv preprint .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 1, 2026. ; https://doi.org/10.64898/2026.04.01.715781doi: bioRxiv preprint .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 1, 2026. ; https://doi.org/10.64898/2026.04.01.715781doi: bioRxiv preprint .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted April 1, 2026. ; https://doi.org/10.64898/2026.04.01.715781doi: bioRxiv preprint