A major locus confers triclabendazole resistance in Fasciola hepatica and shows dominant inheritance

Fasciola hepatica infection is responsible for substantial economic losses in livestock worldwide and poses a threat to human health in endemic areas. The mainstay of control in livestock and the only drug licenced for use in humans is triclabendazole (TCBZ). TCBZ resistance has been reported on every continent and threatens effective control of fasciolosis in many parts of the world. To date, understanding the genetic mechanisms underlying TCBZ resistance has been limited to studies of candidate genes, based on assumptions of their role in drug action. Taking an alternative approach, we combined a genetic cross with whole-genome sequencing to localise a ~3.2Mbp locus within the 1.2Gbp F. hepatica genome that confers TCBZ resistance. We validated this locus independently using bulk segregant analysis of F. hepatica populations and showed that it is the target of drug selection in the field. We genotyped individual parasites and tracked segregation and reassortment of SNPs to show that TCBZ resistance exhibits Mendelian inheritance and is conferred by a dominant allele. We defined gene content within this locus to pinpoint genes involved in membrane transport, (e.g. ATP-binding cassette family B, ABCB1), transmembrane signalling and signal transduction (e.g. GTP-Ras-adenylyl cyclase and EGF-like protein), DNA/RNA binding and transcriptional regulation (e.g. SANT/Myb-like DNA-binding domain protein) and drug storage and sequestration (e.g. fatty acid binding protein, FABP) as prime candidates for conferring TCBZ resistance. This study constitutes the first experimental cross and genome-wide approach for any heritable trait in F. hepatica and is key to understanding the evolution of drug resistance in Fasciola spp. to inform deployment of efficacious anthelmintic treatments in the field.


Introduction
Amongst the helminth infections that pose a substantial risk to livestock and human health worldwide are the liver flukes Fasciola hepatica and F. gigantica. In livestock their impact can be extensive, reducing productivity through lower meat and milk yields, increasing liver condemnation, causing greater susceptibility to other infections, and as a cause of mortality [1][2][3][4][5][6]. In humans it is listed as a neglected tropical disease by the World Health Organisation and estimated that between 2.4 and 17 million people are infected with Fasciola spp. worldwide [7,8]. Historically, optimal control of fasciolosis has been through treatment with the highly effective anthelmintic, triclabendazole (TCBZ); the drug of choice in livestock (Fasinex, Novartis) and humans (Egaten, Novartis), respectively [9,10]. The rising threat of liver fluke infection driven by a changing climate, alterations in land use, enhanced movement of livestock and the ability to encroach into new territories is compounded by a growing problem of TCBZ resistance in livestock [11][12][13][14][15][16]. Similarly, there are increasing reports of the failure of TCBZ to effectively treat Fasciola spp. infections in humans [17][18][19].
Genetic linkage approaches offer a powerful means to map anthelmintic resistance loci, with distinct advantages over candidate gene studies, as no prior knowledge of drug mode of action is required [20]. In trematodes, linkage mapping has identified a sulfotransferase (SmSULT-OR) as the cause of oxamniquine resistance in Schistosome parasites, and in the process revealed its route of action, mode of inheritance and provided a path for future rational drug design [21]. This has allowed global mapping of oxamniquine resistance alleles in natural populations [22,23]. Similarly, genome-wide approaches screening populations of parasites phenotyped for their sensitivity to praziquantel have implicated a transient receptor potential channel (Sm.TRPM PZQ ) in praziquantel resistance in Schistosoma mansoni [24]. There have been similar successes in parasitic nematode species, with population genomic analyses revealing a single genomic quantitative trait locus (QTL) for ivermectin resistance [25] and monepantel resistance [26] in Haemonchus contortus, culminating in the identification of a putative ivermectin resistance gene, HCON_00155390:cky-1, a pharyngeal-expressed transcription factor [27].
Whilst the genetic basis of TCBZ resistance has been a focus of many studies the underlying mechanism remains elusive. A number of candidate genes have been proposed, including βtubulin, P-glycoprotein (Pgp)-linked drug efflux pumps, Flavin mono-oxygenase (FMO), Cytochrome P450 (CYP450), glutathione S-transferase (GST) and fatty acid binding proteins (FABP), as reviewed recently [14,28]. Currently we lack the understanding of whether there is a common mechanism or pathway involved in TCBZ resistance and how TCBZ resistance is inherited, or if the same mechanism is employed by both adult and immature parasites. This inhibits our ability to monitor development of resistance in the field and limits our capacity to effectively deploy anthelmintic drugs to control Fasciola spp. infections.
Herein we demonstrate the first genetic cross and subsequent genomic mapping of a phenotypic trait in Fasciola spp. [29]. We successfully generated an F2 cross between TCBZ resistant (TCBZ-R, FhLivR1) and TCBZ susceptible (TCBZ-S, FhLivS1) F. hepatica parental isolates. Following in vivo phenotyping of F2 parasites and subsequent bulk segregant analysis we identified a~3.2Mbp locus within the F. hepatica genome, comprised of 30 genes, that confers TCBZ resistance. Pooled genotyping of F. hepatica eggs pre-and post-TCBZ exposure in naturally infected sheep confirmed that this TCBZ resistance locus was also under selection in the field. Genotyping of individual parental, F1 and F2 recombinants, revealed that TCBZ resistance is primarily a single locus trait that shows dominant inheritance.

Genetic cross of Fasciola hepatica under experimental conditions
Our capacity to maintain the complete life cycle of F. hepatica in the laboratory and exploit clonal expansion within the snail means genetic crossing and linkage mapping studies are possible for this parasite. However, conducting a genetic cross with a parasite that has an indirect life cycle, is a hermaphrodite with the capacity to self-fertilise, whilst also being genetically diverse, is particularly challenging. There is a need to control for its complex reproductive biology and demography, which we did here using phenotypically defined clones and genotyping individual F1 from single miracidium infection of snails. Crossing of the FhLivS1 and FhLivR1 parentals yielded batches of metacercariae (n = 42), the majority of which (n = 36) were F1 crosses, based on the presence of at least two microsatellite markers from each parent. In most cases (n = 33), the FhLivR1 maternal parent was the source of eggs from which F1 crosses were derived. In total, F1 metacercariae from 28 snails were used to generate F1 adults in vivo and consequently a pool of F2 eggs (Fig 1). To maximise the number of F2 recombinants for in vivo phenotyping we a) performed multiple miracidial infection of snails, b) generated premixed pools of F2 metacercariae from multiple snails prior to infection, c) administered a large F2 metacercarial dose of 400 metacercariae per sheep and d) optimised infection recovery rates (total number of adult parasites recovered from untreated control animals as a proportion of total metacercarial dose administered), which were 21.1% and 22.75%, for Experiment 1 and 2, respectively. Importantly, to determine the impact of TCBZ on genome-wide allele frequency the two pools of F2 used to infect sheep within Experiment 1 and Experiment 2 had a common genetic composition. The number of adult flukes recovered from individual sheep pre-and post-treatment, was significantly different (Fig 1; Mann Whitney W = 25; P = 0.0119 (Experiment 1); 0.00794 (Experiment 2)). When considering all animals within a treatment group for each of the two experiments, drug selection resulted in lower numbers of parasites in TCBZ treated animals, a total of 164 and 119 flukes, compared to the 422 and 455 flukes in untreated hosts, for Experiment 1 and 2, respectively. This represented a reduction of 61% and 74% and an overall recovery rate of 8.2% and 5.95% in treated animals from Experiment 1 and 2 respectively, which constitutes a 2.57-and 3.8-fold reduction in survival of adult parasites in treated hosts compared to untreated controls.

Genome-wide analysis reveals the same scaffolds under selection in both experimental and naturally occurring recombinants
Genome-wide mapping of genetic determinants for phenotypic traits such as drug resistance relies on a well assembled reference genome. We enhanced our previous F. hepatica assembly, increasing scaffold N50 values from 204 Kbp to 1.9 Mb and reducing the number of scaffolds work to produce an F2 cross from FhLivS1 (a clonal population of susceptible parasites) and FhLivR1 (a clonal population of resistant parasites). The parental parasites FhLivS1 and FhLivR1 were produced separately and used to co-infect sheep (n = 2). Some of these parental parasites would cross-fertilise to produce an F1 cross of FhLivS1 and FhLivR1. Eggs were collected from the adult parasites within these sheep. A single miracidium (obtained from these eggs) was used to infect snails (n = 28) and produce clonal F1 populations. The metacercariae were genotyped to ensure they were from an F1 cross and then combined together and used to infect sheep (n = 4). Some of these F1 parasites would cross-fertilise to produce an F2 recombinant population. Eggs were collected from the adult parasites within these sheep. Snails (n = 41 and n = 44 for Experiment 1 and 2, respectively) were exposed to multiple miracidia obtained from these eggs and combined to produce a common pool of F2 metacercariae. For each experiment, two groups of animals were infected with metacercariae from this common pool. Once the infection had reached patency, one group of animals in each experiment was treated with triclabendazole (TCBZ) at a dose of 10mg/kg. At post mortem, those animals which received no treatment had a mixture of triclabendazole susceptible (TCBZ-S) and triclabendazole resistant (TCBZ-R) parasites, whilst those animals that were treated had only TCBZ-R parasites remaining. These parasites were then used for pooled genotyping. (B) A haplotype schematic to show the genetic principle behind the in vivo F2 cross. The F1 cross consists of one haplotype from the susceptible parent: FhLivS1 (FhLivS1.Hap1 or FhLivS1.Hap2) and one haplotype from the resistant parent: FhLivR1 (FhLivR1.Hap1 or FhLivR1.Hap2). In the subsequent F2 generation, recombination events take place and the resistant haplotype becomes introgressed amongst the susceptible haplotype producing an F2 recombinant population for study. (C) Plot to show the reduction in the number of F2 parasites recovered from treated animals compared to untreated animals in Experiments 1 and 2. Boxplot indicates the median number of parasites, upper and lower quartiles, and outliers; overlaid points indicate the number of parasites in each animal. In both experiments a significant difference (Mann-Whitney W = 25 p < 0.05) is seen in the number of F2 parasites from untreated and treated animals. https://doi.org/10.1371/journal.ppat.1011081.g001

PLOS PATHOGENS
A major locus confers triclabendazole resistance in Fasciola hepatica and shows dominant inheritance from 45,354 to 2816, with just 196 scaffolds covering 50% of the genome (Table 1; WormBase ParaSite BioProject PRJEB25283). The completeness of the annotation, as determined by BUSCO, is comparable to that of Schistosoma mansoni (WormBase ParaSite 10; BioProject PRJEA36577). Following discovery and filtering, we identified~9.1M SNPs that segregated between FhLivR1 and FhLivS1 parental clones.
Our approach to mapping loci conferring TCBZ resistance relied on bulk segregant analysis, quantitatively genotyping SNPs in pools of F2 progeny surviving TCBZ treatment and, by comparison with untreated controls, identifying regions of the genome enriched for alleles derived from the resistant parent. In contrast unlinked SNPs (neutral loci) show no difference in allele frequency. We examined the differences in allele frequencies between TCBZ treated (TCBZ+) and TCBZ untreated (TCBZ-) worm pools for each of 9.1M SNPs across the genome, using each sheep as a replicate. The median log-likelihood ratio (LRT) from the generalised linear models (GLM), following bulk segregant analysis from Experiments 1 and 2, is shown in Fig 2A. There was a high degree of concordance between the two experiments, evidenced by over-representation of moving windows of 1000 informative SNPs that independently fell within the 1% highest median LRT in both Experiments 1 and 2 (chi-square test, 1 d.f., 10.952, p < 0.001). We identified 6 scaffolds (13,157,166,324,1853 and 2049) of particular interest because they each had at least 10 moving windows in the top 1% of median LRT in both experiments, suggesting that these were due to a consistent signal of selection within the regions of the genome that they represent. Scaffold 157 showed the greatest evidence of selection (Fig 2A  and Table 2).
We then tested whether the same locus was subject to drug selection under natural field conditions. Bulk segregant analysis of naturally occurring F. hepatica recombinants under drug selection in the field (Field Isolate 1) demonstrated selection of genes on scaffolds 157 and 1853 (Table 3). The median LRT from the GLM following bulk segregant analysis of eggs pre-and post-TCBZ treatment identified scaffolds under selection; scaffold 157 with 206 out of 615 moving windows and scaffold 1853 with 41 of 104 moving windows in the top 1% (Fig 2B and  Table 3). Corroboration between the experimental cross and naturally occurring recombinants, indicating the same two scaffolds (1853 and 157) were under drug selection, supports inheritance of genes on these scaffolds as a means of conferring TCBZ resistance (Fig 2).

PLOS PATHOGENS
A major locus confers triclabendazole resistance in Fasciola hepatica and shows dominant inheritance

Triclabendazole resistance is conferred by a single genomic locus
Given that our F. hepatica genome is comprised of~2800 scaffolds we investigated whether the six mapped scaffolds are linked, by genotyping individual parasites and performing linkage analysis. We genotyped 485 x F2 TCBZ-parasites (S3 Table) with a subset of 48 SNPs derived from each of the six mapped scaffolds and 16 SNPs from neutral (not under selection) scaffolds of comparable size (S4 Table). Linkage is shown by a heat map using |D'| values (Fig 3); pairs of SNPs on scaffolds under selection had high |D'| values (median = 0.942, range = 0.264 to 1) which were usually significant, whilst pairs of SNPs that included neutral scaffolds generally   (Fig 4). SNP markers identified that a single genomic locus, including a 0.3Mbp region of scaffold 1853 and a 2.9Mbp region of scaffold 157, was consistently inherited in resistant parasites (S5 Table). Fasciola hepatica has a 1.25Gbp genome and this 3.2Mbp locus constitutes just 0.25% of the genome, encoding 30 genes.

Triclabendazole resistance shows dominant inheritance
Our experimental genetic cross between a TCBZ-R and a -S isolate confirmed that TCBZ resistance is a heritable trait. By SNP genotyping individual parents, F1 parasites, and 249 x F2 TCBZ+ and 485 x F2-parasites (S3 Table) we could track segregation and reassortment of SNPs through the generations to determine the mode of inheritance. We identified parental SNP genotypes; the TCBZ resistant parent was designated FhLivR1.Hap1/ FhLivR1.Hap2 and

PLOS PATHOGENS
A major locus confers triclabendazole resistance in Fasciola hepatica and shows dominant inheritance

Characterisation of candidate genes within the triclabendazole resistance locus
We identified 30 candidate genes for TCBZ resistance within the~3.2Mbp locus based on our annotated genome and by cross referencing with all available gene annotations ( Table 5). Any of these 30 genes may confer TCBZ resistance, but the data from Field Isolate 1 highlights a strong signal of selection at the start of scaffold 157, identifying a cluster of genes involved in membrane transport, signal transduction and cell signalling, and DNA/RNA binding and transcriptional regulation (Fig 5, genes 3-10; Table 5). Amongst this cluster of genes are those that have been the focus of previous studies on TCBZ action and/or resistance mechanisms, namely an ADP ribosylation factor (Gene 7: ARF, maker-scaffold10x_157_pilon-snap-gene-0.197; Fig 5 and Table 5), a Ras-related protein (Gene 10: Ras-RP, maker-  S4 Table). Analysis of informative resistant recombinant haplotypes (Rec4/5 and Rec7; S5 Table) found within surviving F2 parasites (i.e. those from treated animals) allowed us to further localise the area needed for a parasite to be resistant. In these recombinants, recombination between SNPs delineates a single genomic locus from 1853_3 to 157_6 (~3.2Mbp; 0.3Mbp region of scaffold 1853 and a 2.9Mbp region of scaffold 157) that was consistently inherited in surviving F2 parasites (S5 Table). https://doi.org/10.1371/journal.ppat.1011081.g004

PLOS PATHOGENS
A major locus confers triclabendazole resistance in Fasciola hepatica and shows dominant inheritance scaffold10x_157_pilon-snap-gene-0.182; Fig 5 and Table 5), and an ABCB1 gene (Gene 5: ABCB1, maker-scaffold10x_157_pilon-snap-gene-0.179; Fig 5 and Table 5). Lying within the mapped locus, albeit slightly outside the strongest signal of selection, is another gene that has been the focus of candidate gene studies for TCBZ resistance, a FABP V gene (Gene 17: FABPV, maker-scaffold10x_157_pilon-snap-gene-0.187; Fig 5 and Table 5).
Although several candidate genes fall within the mapped region, it is likely that only one gene within the locus is driving selection. To prioritise amongst gene candidates, we first determined differential expression across life cycle stage-specific transcriptomes (S8 Table). With the exception of the ABCB1 gene (maker-scaffold10x_157_pilon-snap-gene-0.179), two EGFlike proteins (maker-scaffold10x_157_pilon-snap-gene-0.185 and maker-scaffold10x_157_pilon-snap-gene-0.196) and a serine rich protein (maker-scaffold10x_157_pilon-snap-gene-0.188), all the candidate genes were transcribed by the three major F. hepatica life cycle stages, namely newly excysted juveniles (NEJ), immature fluke 21 days post infection and adult fluke (with TPM values ranging from 2-510). The most abundantly transcribed genes were the ADP ribosylation factor (maker-scaffold10x_157_pilon-snap-gene-0.197) and an uncharacterised protein (maker-scaffold10x_1853_pilon-snap-gene-0.13), with highest transcript levels present in adult parasites.
We prioritised candidate genes further by interrogating our genomic and genetic data. There was no evidence for difference in copy number variants (CNV) for the ABCB1 gene (maker-scaffold10x_157_pilon-snap-gene-0.179) and most of our prime candidates (genes 3-10; Table 5) were invariant or contained only synonymous mutations within coding regions. Three non-synonymous SNPs that segregated within the experimental crosses were noted, two within the ABCB1 gene (maker-scaffold10x_157_pilon-snap-gene-0.179 gene), T 830 A and S 852 G, and one within the ADP ribosylation factor (maker-scaffold10x_157_pilon-snap-gene-0.197), C 167 Y (S9 Table). On initial inspection the C 167 Y variant was conserved in other TCBZ-R isolates (FhLivR2, FhLivR3, FhLivR4pop).

PLOS PATHOGENS
A major locus confers triclabendazole resistance in Fasciola hepatica and shows dominant inheritance

PLOS PATHOGENS
A major locus confers triclabendazole resistance in Fasciola hepatica and shows dominant inheritance

A major locus, that shows dominant inheritance, confers triclabendazole resistance
We have demonstrated that TCBZ resistance exhibits Mendelian inheritance and is conferred by a dominant allele at a single locus. This is the first linkage mapping study for any phenotypic trait in F. hepatica and has important implications for our understanding of how drug resistance emerges and spreads in liver fluke populations. A particular strength of our work is the concordance of the classical genetic mapping approach with the field study. We chose to perform the experimental cross with TCBZ-R and -S parasites that were recently isolated from naturally infected sheep in the UK and rendered clonal by laboratory infection of snails [32]. This may explain why we found good agreement between the experimental approach and results from the outbred field populations under natural TCBZ selection. The provenance of both the FhLivR1 clone and Field Isolate 1 places them around 50 miles from one another in the Northwest of the UK, so perhaps such consistency might be expected. Analysis of further isolates will reveal if this genomic locus underpins TCBZ resistance in more geographically dispersed isolates within the UK and beyond. Our approach used pooled genotyping, which enhanced the statistical power and precision of the study [20]. Given the complexity of fluke biology it is difficult to know the final number of F2 recombinants used for in vivo phenotyping in the experimental cross but, based on our experimental design, we can estimate a minimum of 16 and 12 for resistant parasites (TCBZ+) and 42 and 47 for parasites from untreated animals (TCBZ-), in Experiment 1 and 2, respectively. This is broadly consistent with our observation of 113 unique F2 genotypes. The

PLOS PATHOGENS
A major locus confers triclabendazole resistance in Fasciola hepatica and shows dominant inheritance advantage of bulk segregant analysis in our field isolate was that we exploited natural recombination in wildtype populations. The fact that three replicates of relatively small numbers of eggs (500) pre-and post-TCBZ treatment was sufficient to detect signals of selection raises the exciting prospect of conducting similar studies for TCBZ selection in F. hepatica field populations in other geographical locations and provides a valuable approach for the study of flukicide resistance more broadly, e.g. for drugs such as closantel and albendazole. Our work was conducted with adult parasites and clearly shows TCBZ resistance is a heritable trait. One of the most important aspects of our classical linkage mapping is that it allowed us to determine that, in contrast to oxamniquine and praziquantel resistance in schistosomes, TCBZ resistance is a dominant trait [24,33,34]. This tells us that once resistance emerges or is introduced within liver fluke populations it has the potential to spread rapidly [35], and highlights the need for rapid detection and effective treatment to mitigate the impact of TCBZ-R liver fluke infections in livestock and humans.

PLOS PATHOGENS
A major locus confers triclabendazole resistance in Fasciola hepatica and shows dominant inheritance Linkage mapping has proved highly successful for identifying genetic determinants for phenotypic traits such as pathogenicity, host specificity and drug resistance in parasites of humans e.g. protozoa [36][37][38][39][40] and Schistosoma spp. [21]. More recently, population genomic analyses mapped a QTL for ivermectin resistance in the ruminant nematode, H. contortus [25]. Mapping studies in these parasites benefited from fully assembled genomes, something we are yet to achieve for F. hepatica, although recent publication of a chromosomal-level genome assembly for its sister species F. gigantica is encouraging [41]. Genomic resources for F. hepatica extend to two independent assemblies [42][43][44]. Our success at generating the first genetic cross and subsequent linkage mapping of drug resistance loci paves the way for studies on important phenotypic traits for F. hepatica in the future. The technical challenges presented by a fragmented genome have been highlighted elsewhere [25], and whilst we overcame many of these by our experimental cross and additional linkage experiments, a chromosome-level assembly will be a vital resource to progress future studies.

Genetic mapping pinpoints candidate genes conferring triclabendazole resistance
Anthelmintic resistance can occur due to increased efflux, enhanced metabolism and through efficient detoxification mechanisms. By integrating our mapping studies and our genomic and genetic data with genes previously implicated in TCBZ resistance or TCBZ mode of action, we can prioritise specific genes that may play a role. Our mapping studies have highlighted that any of 30 genes could be involved in TCBZ resistance, but strongest selection was placed on the region where ABCB1, RAS-RP, ARF and a few other genes cluster. ABCB1 (P-glycoprotein, Pgp), also known as MDR1, is implicated in drug resistance in multiple organisms. Overexpression of Pgp transporters, leading to increased drug efflux has been proposed as a potential route to drug resistance [45,46], and the observation that we are dealing with a dominant trait is consistent with a role for over expression of ABC transporters. In F. hepatica ABCB1 (Pgp)-linked drug efflux pumps have been the focus of altered drug uptake studies. Existing evidence of a role for Pgp in TCBZ resistance includes a) lower uptake of TCBZ and its metabolite TBCZ.SO in TCBZ-R compared to TCBZ-S flukes [47,48], b) reversal of the resistance phenotype in vitro by co-incubation with ivermectin (IVM), a known multidrug resistance (MDR) reversing agent and potential competitive substrate for Pgp [47], and c) potentiation of TCBZ action in vitro in TCBZ-R flukes in the presence of Pgp inhibitor R(+)-verapamil [49-51]. Therefore, our identification of ABCB1 is noteworthy, but based on our current annotation, there is no support for CNV that are thought to underlie overexpression of ABC transporters, and the lack of constitutive expression of this ABCB1 on scaffold 157 in adult parasites is inconsistent with a role in TCBZ resistance. Specific

PLOS PATHOGENS
A major locus confers triclabendazole resistance in Fasciola hepatica and shows dominant inheritance FABPs are small proteins that can bind anthelmintics [57] and they have been shown to be involved in drug storage and sequestration [58]. Upregulation of FABP mRNA was noted when drug resistant Anopheles gambiae were exposed to permethrin [59]. FABPs are known to be present in the tegument of Fasciola spp. [60-62]. In 2016, a systems approach revealed a FABP superfamily of seven clades, including the novel identification of the FABPV family, a representative of which is within our genomic locus [63]. The FABPV is closely related to FABP isoforms I-III [63]. Previous proteomic comparisons showed reduction in FABP synthesis (encoded by three FABP genes distinct from the FABPV gene located on scaffold 157) in a susceptible isolate exposed to TCBZ [64]. Moreover, a type I FABP Fh15 with the capacity for sequestration showed increased expression in resistant adult flukes exposed to TCBZ [64].
Ras-RP and ARF have not been implicated in TCBZ resistance based on previous candidate gene studies but are key regulators of important biological processes. The presence of a classical Ras gene, (Ras-RP) and another Ras superfamily member, ARF, within the major locus associated with TCBZ resistance is of interest for several reasons. A subfamily of Ras genes, Rabs, are small GTPases that have been linked to drug resistance in the protozoan parasite, Leishmania donovani [65]. In yeast, TCBZ has been shown to inhibit the production of cAMP by either direct inhibition of adenylate cyclase or by acting on the GTP-Ras-adenylyl cyclase pathway [66]. Fasciola hepatica adenylate cyclase is amongst the most active of any organism, its activity is thought to regulate carbohydrate metabolism and motility of the worms [67]. Adenylate cyclase in F. hepatica is activated by serotonin receptors that function through GTP-dependent transmembrane signalling pathways [68][69][70][71] and was identified as a potential therapeutic target in F. hepatica several decades ago [72]. Experiments with liver fluke tissue revealed that an endogenous ADP-ribosylation enzyme and its protein substrate were present and capable of regulating adenylate cyclase activity [73]. Our observation of a C 167 Y variant in ARF (maker-scaffold10x_157_pilon-snap-gene-0.197), that segregated in the genetic crosses and was conserved in other TCBZ-R adult fluke isolates (FhLivR2, FhLivR3, FhLivR4pop) is of particular interest and warrants further investigation.
Whilst these analyses may help us narrow down which gene might be responsible drug resistance mechanisms are not restricted to mutations in coding regions and changes in gene expression. It is possible that any one of the 30 genes within the locus is responsible for resistance and, given that RNAi has been optimised for F. hepatica [74], systematic knockdown of each candidate gene would be a sensible way forward. RNAi, combined with recent advances in the culture of juvenile parasites [75] and in vitro phenotyping for TCBZ resistance, offers a powerful platform with which to screen for the causal gene, and provides opportunity to investigate whether resistance mechanisms are stage specific. Interrogation of in vivo RNA-seq datasets from isolates of known phenotype would inform on whether any candidate genes show differential expression on TCBZ exposure. Similarly, as small non-coding microRNA (miRNA) are known to regulate gene expression it would be prudent to look for predicted miRNA binding sites in candidate genes, particularly given that they have been linked to drug resistance in nematodes [76] and a large dataset of miRNAs has been reported in F. hepatica [77].
Undoubtedly one or more of these approaches will allow us to pinpoint the causal gene, but the question remains as to whether this underlying mechanism explains all observations of phenotypic resistance. We note that the signal of selection in the experimental cross and field data are adjacent rather than coincident. Out of necessity the experimental and field data used different sets of SNPs so this could be a statistical artefact. Alternatively, it may be a biological effect, indicating that different mutations circulate in the field that target the same genetic locus. This raises the interesting possibility that resistance can evolve multiple times but is constrained in the number of genome targets that can confer resistance. With the genomic

PLOS PATHOGENS
A major locus confers triclabendazole resistance in Fasciola hepatica and shows dominant inheritance mapping approaches optimised here we now have the tools to address these more complex questions about TCBZ resistance in F. hepatica.

Gene content in the major locus does not support a direct role for many gene families formerly identified as candidates for triclabendazole
In the absence of a genome-wide approach, TCBZ resistance studies relied on assumptions about involvement of candidate genes or gene families [28]. Across the genome we identified multiple candidate genes: 14 tubulin genes, 25 ABC transporter genes, three CYP450 (-like) genes, seven FABP genes, 11 glutathione S-transferase genes and three thioredoxin peroxidase (-like) genes in the F. hepatica genome (S7 Table). Most of these genes were located in scaffolds that showed no evidence of being under selection in our experiment, with none of the moving windows appearing in the top 1%. Furthermore, only genes on scaffold 157 showed evidence of being under selection in both our experimental and field data (S7 Table). Amongst prime candidate genes that can be excluded based on their absence within the major locus of the populations studied here are β-tubulin, the microtubule fraction known to cause BZ resistance in nematodes [78]. Although they were initially implicated in TCBZ resistance based on changes typical of microtubule inhibition in TCBZ-S but not -R flukes, [reviewed by [79][80][81] no differences in β-tubulin isotypes sequences or expression levels were reported between TCBZ-S and TCBZ-R flukes [64,82,83]. Whilst it may still be the case that TCBZ acts via β-tubulin, the lack of a β-tubulin gene in our mapped region rules this gene out as a candidate for directly conferring TCBZ resistance. Similarly, it has been shown that drug metabolism is upregulated in TCBZ-R flukes [48,84] possibly involving FMO, CYP450 or GST, the mu type, specifically [85][86][87], and an amino acid substitution T 143 S of GST in the TCBZ-R flukes has been reported [88]. The absence of GSTs, FMO or CYP450 from the locus excludes the direct action of these molecules in TCBZ resistance at least in the populations studied here. It is important to note that the mapping approach was taken with populations within a restricted geographic region of the UK and it may be that resistance is driven by different processes in other locations. Our work provides the first means with which to address whether a common mechanism of resistance occurs in F. hepatica populations.

Conclusion
TCBZ is the drug of choice to treat fasciolosis in sheep and cattle, and is the only drug licenced to treat humans. Identifying genetic determinants for resistance, as we have here, is invaluable to our understanding of the mechanisms behind TCBZ resistance and how we might best mitigate its impact. In this study, we exploited the biological process of clonal expansion within the snail intermediate host and recent advances in large sequencing datasets for F. hepatica to further our understanding of the genetic mechanisms involved in TCBZ resistance. We have shown 1) that TCBZ resistance is primarily a single locus trait that shows dominant inheritance; 2) we have performed the first experimental genetic cross and linkage mapping study for any phenotypic trait in F. hepatica; 3) we successfully applied bulk segregant analysis of eggs pre-and post-treatment to detect signatures of selection within field isolates of F. hepatica and 4) we have conducted the first genome-wide analysis of TCBZ resistance. We have identified a small number of genes involved in membrane transport, (e.g. ABCB1), transmembrane signalling and signal transduction (e.g. Ras-RP, ARF and EGF-like proteins), DNA/RNA binding and transcriptional regulation (SANT/Myb-like DNA-binding domain protein) and drug storage and sequestration (e.g. FABP) as prime candidates for conferring TCBZ resistance. Detecting a signal of selection in naturally infected, live animals in the field provides a blueprint to determine if a common mechanism of TCBZ resistance is adopted by demographically

PLOS PATHOGENS
A major locus confers triclabendazole resistance in Fasciola hepatica and shows dominant inheritance distinct F. hepatica populations and paves the way for molecular tests to detect drug resistant parasites and more effectively target treatments in both livestock and humans.

Ethics statement
All applicable institutional, national, and international guidelines for the care and use of animals were followed. Experimental infection in sheep was conducted under Home Office Licence PPL 40/3621 and PE77BFD98 in accordance with Animal (Scientific Procedures) Act 1986 and ethical approval for the field study was provided by the University of Liverpool Veterinary Research Ethics Committee (VREC582).

Enhanced assembly and annotation of the Fasciola hepatica genome
The published assembly of F. hepatica [43] was improved using a combination of Hi-C and linked read data. High molecular weight DNA was prepared from FhLivS1 adult fluke [32] using either Genomic-tip (Qiagen, UK) or an adapted lithium chloride and Triton X-100 lysis and phenol-chloroform extraction [89,90]. Hi-C libraries were prepared by Dovetail Genomics (Santa Cruz, CA, USA) to generate 174 million Illumina paired-end reads and scaffolded using Hi-Rise [91]. Further scaffolding was performed using linked reads from a 10X Chromium platform using Illumina reads that mapped to within 20kbp of the end of a scaffold. Scaffolds were joined where they exhibited at least 10 linked reads connecting a pair of scaffolds and where the number of links between a pair of scaffolds was at least twice as many as the next best connection. Gaps within scaffolds were then filled where possible using 2x250 bp reads [43] assembled into contigs with Discovar [92], followed by further gap filling and polishing using Illumina 2x100bp and 2x250bp reads with Pilon [93].
Annotation was performed using MAKER2 [94]. RNA-seq data [43] were used to provide initial transcript predictions by running BRAKER [95] and StringTie [96,97] to generate low quality transcript predictions and (from BRAKER) to train AUGUSTUS [95]. RepeatMasker (http://www.repeatmasker.org) was used to identify repeat regions. SNAP [98] was trained in three iterative runs of MAKER2 [94]. The completeness of the set of predicted proteins was assessed using BUSCO v3 [99] against its set of Eukaryota reference proteins and compared with the predicted proteins from Schistosoma mansoni (WBPS10; PRJEA36577).

Pooled genotyping of phenotyped adult F2 populations derived from an experimental cross
The genetic crossing of a clonal TCBZ-R and -S isolate was carried out using the FhLivR1 and FhLivS1 isolates. Provenance, validation of phenotype and microsatellite genotyping of these two isolates was described previously [32]. The genetic cross experimental approach is shown in Fig 1. When generating F2 populations this approach required selection of F1 parasites from mating events between the two parental isolates, rather than mating between parasites of the same genotype, or self-fertilization. This was done by generating multiple single F1miracidium:snail infections of our laboratory-maintained Galba truncatula and screening F1 metacercariae by genotyping them for the presence of both parental multilocus genotypes [32,100]. In Experiment 1, two pools of adult flukes derived from a common population of F2 recombinants were generated by in vivo phenotyping in 10 sheep. This gave rise to one pool (F2 TCBZ-), comprised of a mixture of TCBZ-S and -R flukes (from untreated sheep, n = 5) and another pool (F2 TCBZ+), comprised of only TCBZ-R flukes (from treated sheep, n = 5). This process was repeated in Experiment 2, using a second, common, pool of F2 eggs (Fig 1). Each
High quality SNPs were identified from whole-genome resequencing of five isolates [43] using Bowtie2 [101] under sensitive settings and GATK [102]. SNPs were filtered to select high confidence SNPs (i.e. those that segregated within parental isolates FhLivR1 and FhLivS1 and F2 progeny [32,43,100], had a quality score >100 and a depth of between 6 and 50 for each isolate. Only biallelic SNPs were used. Following discovery and filtering we identified a panel of 9M SNPs (SNP panel 1). Illumina TruSeq libraries were generated from DNA from each F2 pool and sequenced with 2x125bp reads on an Illumina HiSeq2000. Illumina adapter was removed using Cutadapt v1.2.1 [103] and reads further trimmed with a minimum window quality score of 20 with Sickle v1.200 (https://github.com/najoshi/sickle). The counts for reference and alternate alleles in each F2 pool were generated using SAM Tools mpileup [104] and filtered to retain SNPs with coverage depth within the 10% and 90% quantiles (~7.7M SNPs). Generalised linear models (GLM) with a binomial error distribution were calculated for each SNP in R (https://www.R-project.org/) for each of the two experiments. Moving windows (containing 1000 SNPs and advanced by 100 SNPs) were calculated to give the median loglikelihood ratio (LRT) statistic associated with the difference in allele frequency between parasites in F2 TCBZ+ and F2 TCBZ-pools. Windows having a median LRT in the upper 1% quantile within each experiment were identified, and only those windows exhibiting a median LRT in the top 1% for both experiments 1 and 2 were taken forward for further analysis.

Pooled genotyping of eggs pre and post triclabendazole treatment from a field population of Fasciola hepatica
Our F. hepatica faecal egg count reduction test (FECRT) provided a common pool of eggs from three replicate groups of 10 sheep, with the same 10 sheep sampled pre-and 21days post-TCBZ treatment [105]. Field Isolate 1 (Cumbria, UK) had a total pre-treatment egg count of 15817 (4052, 2971 and 8794 across the three groups) and the total post-treatment egg count was 3187 (1037, 978 and 1172 across the three groups). This equates to an 80% reduction and indicates the presence of treatment failure. Five hundred eggs were collected from each of the six samples and washed five times in 1ml of ddH 2 O before being used for DNA extraction.
Given the genetic diversity inherent in fluke populations it was necessary to increase our panel of high-quality SNPs, by including SNPs previously identified in FhLivSP, FhLivR2,and FhLivR3 [32,43], and by resequencing the genome of six individual F. hepatica from isolate FhLivR4pop, a TCBZ-R F. hepatica population from South Wales, UK. SNPs were identified using BowTie2 [101] under sensitive settings and GATK [102] and filtered based on a quality score greater than 100 and a depth of between 6 and 50 for each isolate. This provided a~21M SNP panel (SNP panel 2). Sequencing of eggs was performed by the Centre for Genomic Research, University of Liverpool using the NovaSeq S2 Flowcell (Illumina). The GLM procedure described above was used to compare SNPs (~14M) from pre-treatment (eggs obtained from Day 0 faecal samples) and post-treatment (egg obtained from Day 21 faecal samples).

Linkage analysis of scaffolds under selection
We genotyped individual parasites: 249 x F2 TCBZ+ and 485 x F2 TCBZ-parasites, 45 x F1 parasites and ten parental (FhLivR1 and FhLivS1) parasites (S3 Table). To determine if the scaffolds under selection are linked, we analysed genotypes from 485 x F2 TCBZ-(untreated) parasites.

PLOS PATHOGENS
A major locus confers triclabendazole resistance in Fasciola hepatica and shows dominant inheritance Genotyping was performed on a subset of 48 SNPs, from scaffolds under selection, and 16 SNPs, from scaffolds of comparable size not under selection (neutral scaffolds), were selected from SNP panel 1 and further filtered for coverage depth within the 20% and 80% quantiles. There was a preference for SNPs in exons, they were selected along the entire length of the scaffold, and they had 50bp of conserved sequence either side, to allow primer design (S4 Table). Assay design and genotyping was conducted by LGC Genomics (Hertfordshire, UK) using KASP genotyping chemistry. It was not possible to design assays for nine SNPs and after genotyping three SNPs (13_5, 13_6 and 917_3) showed monomorphic results and were not included in subsequent analyses (S4 Table). PHASE 2.1.1 [30,31] was used to infer haplotypes from SNP data. After an initial analysis, scaffolds under selection were orientated to minimise recombinant events and PHASE was rerun with a 95% confidence cut-off, 1000 iterations, thinning interval of 10 and burn-in of 100. Haplotypes of the neutral scaffolds were inferred separately to those under selection and run with the same parameters. Arlequin 3.5.1.3 [106] was used to assess linkage disequilibrium and calculate |D'| values between all pairs of SNPs; each genotype was represented once to avoid duplication of genotypes from clones. The number of steps in the Markov chain was 100000 and the number of dememorization steps (burn-in) was 5000. False discovery rate correction [107] was used to correct p-values in R 3.0.1 (https://www.R-project.org/), a significance level of p < 0.05 was used. The R package ggplot2 was used to plot results.

Inheritance patterns (segregation) and finer scale mapping of triclabendazole resistance genes
To track segregation and association of SNPs from parental haplotypes to recombinant F2s, numbers of haplotypes across the region under selection and for each scaffold were identified in control (TCBZ-) and treated (TCBZ+) animals and assigned to a parental genotype. We used these to determine whether one or both parental haplotypes could confer resistance and whether resistance was a dominant or recessive trait. To further localise the region associated with resistance, recombinant haplotypes were identified and recombination between SNPs used to delineate a region always inherited by parasites that survived TCBZ treatment (i.e. resistant parasites).

Annotation of genes in region of genome under selection
The protein sequence of genes under selection (candidate genes) were run through UniProt Blast using the UniRef50 and UniProtKB_RefProtSwissProt databases [108], WormBase Para-Site Version WBPS14 (WS271) Blast against all species in the protein database [43,44,109,110] and OrthoDB version 9 against the Metazoan database [111] to determine an appropriate description and function for each candidate gene. InterPro [112] was used to identify domains as a confirmation of the protein function. WormBase ParaSite Version WBPS16 (WS279; [109,110]

PLOS PATHOGENS
A major locus confers triclabendazole resistance in Fasciola hepatica and shows dominant inheritance 10 ( Table 5). The focus was non-synonymous changes in segregating resistance alleles that were conserved in related parasites (FhLivR2, FhLivR3, FhLivR4pop).

Annotation of previously identified candidate genes
We interrogated the F. hepatica genome for candidate genes that had been the focus of previous TCBZ resistance studies [28]. Tubulin proteins were as previously characterised [116] and a list of ABC transporters was provided courtesy of A. Maule, personal communication (S7 Table). NCBI nucleotide and protein databases were searched for F. hepatica and either (i) cytochrome P450, (ii) fatty acid binding protein, (iii) glutathione S-transferase, or (iv) thioredoxin peroxidase. WormBase ParaSite Blast [109,110] was used to identify candidate genes within the F. hepatica genome (DNA and protein database of BioProject PRJEB25283). Genes were not included where the protein was only structurally related to the functional annotation, or only contained domains related to gene function.

In vivo experimental infections
Sheep infections were carried out essentially as described previously [32]. Briefly, >12 weekold Lleyn cross lambs were infected by oral administration of~200 (parental clones and F1) or 400 (F2) metacercariae per sheep. Infection status was monitored weekly by ELISA [117]

DNA isolation
For F2 pooled genotyping, genomic DNA was extracted from~20mg at the anterior end of each adult fluke using the DNeasy Blood and Tissue kit (Qiagen, UK) with elution in 100μl of buffer AE. This was followed by precipitation using 3M NaOAc and isopropanol at 4˚C. Each individual fluke DNA was checked for quality on a 2% agarose gel and quantified by Quant-IT PicoGreen (Life technologies/ThermoFisher Scientific). Equimolar concentrations of genomic DNA from each parasite was mixed and purified with GenomicTip (Qiagen, UK) to create an F2 pool of high molecular weight DNA per sheep, for sequencing. Egg DNA from field isolates was extracted using the DNeasy Blood and Tissue Kit (Qiagen, UK) with the following modifications: (i) a micropestle (Argos Technologies, USA) and Pellet Pestle Motor (Kontes) were used to homogenise the eggs; (ii) RNase A was used; (iii) elution was in 100μl of buffer AE. 5μl of egg DNA was subjected to whole genome amplification using a REPL-g Mini Kit (Qiagen, UK), followed by purification using QiaAmp Mini Column (Qiagen, UK) with the following modifications (i) the DNA was added to the spin column and then the wash steps were performed (ii) elution was with 65μl of buffer AE.

Maintenance of Fasciola hepatica in Galba truncatula
Galba truncatula snail stocks were maintained as described previously [32]. Briefly, they were maintained at 22˚C on pans of clay mud and fed on a diet of Oscillatoria spp. algae. F. hepatica eggs were

PLOS PATHOGENS
A major locus confers triclabendazole resistance in Fasciola hepatica and shows dominant inheritance embryonated at 27˚C for 14 days and stimulated to hatch by exposure to light. Each snail~4 mm in height was exposed to either one miracidium (F1) or 5-8 miracidia (F2) to generate pools of clonal or multi-genotype parasites, respectively. Following infection, snails were maintained on mud, fed every 2-3 days and stimulated to shed cercariae by sealing the snail into visking tubing containing water and exposing them to a drop in temperature. The cercariae then encyst on the visking tubing as metacercariae [32]. Metacercariae from multiple snails were pooled prior to infection, to provide a dose rate of~400 F2 parasites per sheep; 7-10 F2 metacerariae/snail (n = 41snails) for Experiment 1 and 4-10 metacerariae/snail (n = 47snails) for Experiment 2 (Fig 1).
Supporting information S1

PLOS PATHOGENS
A major locus confers triclabendazole resistance in Fasciola hepatica and shows dominant inheritance