Background of H. contortus populations and the original backcross
The MHco3(ISE), MHco4(WRS) and MHco10(CAVR) are H. contortus populations maintained and stored at the Moredun Institute. MHco3(ISE) was originally derived by multiple rounds of inbreeding of the SE population [46]. The precise provenance of SE is not clearly recorded but is thought to be originally obtained from East Africa. MHco10(CAVR) is derived from the Chiswick Avermectin Resistant (CAVR) strain, an ivermectin resistant field population originally isolated in Australia [54]. MHco4(WRS) is derived from the ivermectin resistant White River Strain (WRS) originally isolated from South Africa [14].
The experiment described here extends two previously-described backcross experiments. The construction, phenotypic validation, and initial microsatellite analysis of these backcrosses has been described previously [39], and is outlined in Fig. 1a and b. Briefly, two independent crosses were performed in parallel: (i) MHco3(ISE) x MHco10(CAVR), and (ii) MHco3(ISE) x MHco4(WRS). In the first generation of each cross, female worms of the ivermectin resistant parental populations [MHco10(CAVR) or MHco4(WRS)] were crossed with male MHco3(ISE) worms, to generate lines designated MHco3/4 and MHco3/10, respectively (Fig. 1a). After the first cross, the F1 generation females were mated to male MHco3(ISE), resulting in backcross generations designated MHco3/10.BCn and MHco3/4.BCn, where n is the number of backcross generations (BC). This involved the recovery of immature worms from the abomasum at day 14 post infection by necropsy of donor sheep; worms were washed in physiological saline (0.85% NaCl) after which their sex was determined [55], followed by surgical transfer of 45–100 male MHco3(ISE) and 50–100 female F1 L4/immature adult worms into worm-free recipient lambs all within 2 h of the original collection. This process of backcrossing was repeated for a total of 4 generations, i.e., MHco3/10.BC4 and MHco3/4.BC4 (Fig. 1b). In vivo ivermectin selection was applied after mating and before collection of eggs to enrich for ivermectin resistant adults to be used in the subsequent round of the backcross. A controlled efficacy test was performed on the three parental-, MHco3/4.BC4, and MHco3/10.BC4 populations to determine the initial levels of ivermectin efficacy, and the degree of resistance acquired in the introgressed lines. After four rounds of backcrossing, a selection experiment was performed. MHco3/10.BC4 and MHco3/4.BC4 populations were used to infect recipient sheep, after which eggs were collected and in vitro culture to L3. These infective larvae were used to infect two sheep per population, one that received 0.1 mg/kg ivermectin and one that remained drug naive. At 7 days post treatment, treated and naive L4/immature adults were recovered by necropsy and stored for molecular analysis.
The backcrosses described above were extended by performing a further four rounds of in vivo ivermectin selection, during which mating within the BC4 population continued (Fig. 1c). For each cross, progeny of the BC4 generation were cultured to L3 and used to infect a donor sheep, which was subsequently treated with ivermectin (0.1 mg/kg in round 1, followed by 0.2 mg/kg in following rounds). Eggs from ivermectin-treated survivor adults were collected post-treatment, cultured to L3, and used to infect a new donor sheep. L4/immature adults were collected by necropsy after rounds three (i.e., BC4.IVM.P3) and four (i.e., BC4.IVM.P4) of selection and stored for molecular analysis.
The sheep used in this study were born and raised indoors at the Moredun Research Institute under worm free conditions. For those recipient sheep requiring surgery: anaesthesia was induced either by intravenous thiopentone injection, or using a halothane and oxygen mask. The sheep were then intubated and anaesthesia maintained using halothane and oxygen. Sheep were routinely injected with a non steriodal anti-inflammatory agent (1 mg/kg meloxicam; Metacam 20 mg/ml solution for injection; Boehringer Ingelheim) and antibiotic (7 mg/kg amoxicillin/1.75 mg/kg clavulanic acid, Synulox ready-to-use injection; Pfizer) and closely monitored on completion of the surgery. No adverse effects were noted. Sheep were humanely killed under schedule 1 to the Animals (Scientific Procedures) Act.
Library preparation and sequencing
Whole genome sequencing was performed on the three parental populations (MHco3[ISE], MHco10[CAVR] & MHco4[WRS]), and from each cross, pre and post ivermectin treatment following four rounds of backcrossing (MHco3/10.BC4.noIVM, MHco3/4.BC4.noIVM, MHco3/10.BC4.IVM, MHco3/4.BC4.IVM), and after 3 and 4 rounds of additional selection after the backcross (MHco3/10.BC4.IVM.P3, MHco3/4.BC4.IVM.P3, MHco3/10.BC4.IVM.P4, MHco3/4.BC4.IVM.P4). Pools of male and female worms were included for the parental (n = 50–60 worms) and BC4 samples (n = 25–30 worms), whereas only female worms were used in the passage 3 & 4 samples (n = 60 worms). Genomic DNA from each pooled sample was prepared using a phenol chloroform extraction protocol as previously described. Sequencing libraries were prepared using a PCR-free protocol [56], and sequenced as described in Additional file 12: Table S1. We generated approximately 6.14 × 1011 bp of sequence data, which equates to approximately 199.65× unmapped genome coverage per sample. Raw sequence data quality was analysed using FASTQC (
http://www.bioinformatics.babraham.ac.uk/projects/fastqc/
) and was assessed using MultiQC [57] prior to further processing.
Mapping and variant analysis
Raw sequence data was mapped per lane to the version 3 reference genome (available here: ftp://ngs.sanger.ac.uk/production/pathogens/Haemonchus_contortus) using BWA-MEM [58](default parameters with -k 15). Samples for which multiple lanes of mapped data were available were merged, duplicate reads were marked and removed using Picard v2.5.0 (https://github.com/broadinstitute/picard), and mapped perfect read pairs were extracted per sample (samtools-1.3 view -f 2).
Genome-wide variants were determined for each sample using samtools-1.3 mpileup (−F 0.25 -d 500). Nucleotide diversity and Tajima’s D was determined using npstat [59], which required pileup files (which we derived from the mpileup generated above) that were split per sample per chromosome as input. We determined short-range linkage disequilibrium between pairs of variants in single and paired reads using LDx [60]. We compared these data with estimates of LD decay over genetic distance as proposed by Hill and Weir [60, 61]. Popoolation2 [62] was used for the analysis of pairwise genetic diversity throughout the genome. Briefly, the mpileup file was converted into a synchronised file (popoolation2 mpileup2sync.jar --min-qual 20), after which indels + 5 bp (popoolation2 identify-indel-regions.pl −−min-count 2 −−indel-window 5) and repetitive and difficult to map regions (separately identified by repeat masking the genome [http://www.repeatmasker.org/]) were excluded (popoolation2 filter-sync-by-gtf.pl). The synchronised file was used as input to determine pairwise FST which was calculated in 10 kbp windows throughout the genome (popoolation2 fst-sliding.pl --pool-size 50 --window-size 10,000 --step-size 5000 --min-count 4 --min-coverage 50 --max-coverage 2%). Similarly, per base comparisons were determined using a Fisher’s exact test (popoolation2 fisher-test.pl --min-count 4 --min-coverage 50 --max-coverage 2% --suppress-noninformative).
Copy number variation between parental strains was determined using cnv-seq [63]. Best read mapping hits were determined per bam (samtools-1.3 view -F 4 bam | perl -lane ‘print “$F[2]\t$F[3]”‘> bam.hits), before read counts for pairwise comparisons determined in 10 kbp sliding windows using cnv-seq.pl. The R package cnv (v 0.2–8) was used to determine pairwise chromosome-wide normalised log2 ratios. Characterisation of structural variation in the parental populations was performed using the speedseq sv pipeline [64] to identify putative duplications, deletions and inversions in each population. This approach exploits all reads mapped using BWA-MEM, including split and discordant read pairs, which are subsequently scored using LUMPY and SVTyper. Conservative filtering was applied to retain only homozygous variants (1/1) with a minimum quality score of 100.
Population genetic modelling
Bespoke single-locus and multi-locus models were used for population genetic inference from the data. Common to both methods, variants in the genome potentially associated with resistance were identified. At each locus in the genome, a nucleotide {A,C,G,T} was defined as existing in the susceptible parental population if it was observed at a frequency of at least 1% in that population. Loci in the resistant parental population were then identified for which exactly one nucleotide that did not exist in the susceptible population was observed at a frequency of at least 1%; without prejudice as to its phenotypic effect, this was denoted the ‘resistant’ allele at that locus, being associated with the resistant parental population.
Single-locus model
A single-locus population genetic model was then used to identify variants with frequencies that were inconsistent with selective neutrality. The neutral expectation was calculated using a Wright-Fisher evolutionary model to simulate the progress of a variant allele, modelling the probability distribution of the frequencies of homozygous resistant (q1), heterozygous (qh), and homozygous susceptible (q0) individuals throughout the course of the experiment under the assumption of selective neutrality. The Wright-Fisher model uses a simple multinomial process for propagation; if the next generation of the population contains N individuals, the probability of obtaining n1, nh, and n0 individuals with each diploid variant is given by:
$$ P\left({n}_1,{n}_h,{n}_0\right)=\frac{N!}{n_1!{n}_h!{n}_0!}{q}_1^{n_1}{q}_h^{n_h}{q}_0^{n_0} $$
(1)
At the stages in the experiment where the population was comprised of eggs, we assumed N to be large, such that the processes of crossing and backcrossing could be described deterministically. Under these circumstances the genotype frequencies following a cross are given by:
$$ {\displaystyle \begin{array}{c}{q}_1^{\prime }={q}_1^2+{q}_1{q}_h+0.25{q}_h^2\\ {}{q}_h^{\prime }={q}_1{q}_h+2{q}_1{q}_0+0.5{q}_h^2+{q}_h{q}_0\\ {}{q}_0^{\prime }={q}_0^2+{q}_h{q}_0+0.25{q}_h^2\end{array}} $$
(2)
Where a prime denotes the frequency in the next generation. Similarly, the genotype frequencies following a backcross with the homozygous susceptible population are given by:
$$ {\displaystyle \begin{array}{c}{q}_1^{\prime }=0\\ {}{q}_h^{\prime }={q}_1+0.5{q}_h\\ {}{q}_0^{\prime }={q}_0+0.5{q}_h\end{array}} $$
(3)
Evaluating this process gives a probability distribution for the frequency of the resistant allele at any given point in the experiment, dependent upon the number of resistant alleles in the initial resistant population. This number of resistant alleles, which we denote nr, is unknown; sequence data from the resistant parental population were used to generate a prior distribution for this value.
We first suppose that the frequency of the resistant allele in the resistant parental population is equal to some value, pr. Given that the initial experimental population contains 50 resistant worms, with 100 alleles, collected from this larger population, the distribution of nr.
is then given by:
$$ P\left({n}_r=j\right)=\left(\begin{array}{c}100\\ {}j\end{array}\right){p}_r^j{\left(1-{p}_r\right)}^{100-j} $$
(4)
We now suppose that sequencing the resistant parental population gave ar resistant alleles at the locus in question out of a total read depth of Ar. The initial frequency pr is unknown; however, using the data an estimate can be made for this statistic, expressed as a distribution of the allele frequency. Specifically, under the assumption of a uniform prior, the underlying probability pr can be said to be distributed as a beta distribution with parameters ar and Ar - ar + 1. We therefore obtain the following distribution for the statistic nr:
$$ P\left({n}_r=j\right)={\int}_0^1\frac{A_r!}{\left({a}_r-1\right)!\left({A}_r-{a}_r\right)!}{p}_r^{a_r-1}{\left(1-{p}_r\right)}^{A_r-{a}_r}\left(\begin{array}{c}100\\ {}j\end{array}\right){p}_r^j{\left(1-{p}_r\right)}^{100-j}{dp}_r $$
(5)
for any value of j between 0 and 100. Values of this distribution were calculated using numerical integration. This process generated the neutral expectation of the resistant allele frequency conditional on the observation of this frequency in the resistant population. If at the sampling point t, a total of Nt worms were collected, the distribution of the number of worms nt with the resistant allele at that point is given by:
$$ P\left({n}_t=k\right)=\sum \limits_jP\left({n}_t|{n}_r=j\right)P\left({n}_r=j\right) $$
(6)
for all 0 ≤ k ≤ Nt.
The extent to which observed allele frequencies were consistent with the neutral model was calculated using a likelihood model:
$$ {L}_t=\sum \limits_{k=0}^{N_t}P\left({n}_r=k\right)\left(\begin{array}{c}{A}_t\\ {}{a}_t\end{array}\right){\left(k/{N}_t\right)}^{a_t}{\left(1-\left(k/{N}_t\right)\right)}^{A_t-{a}_t} $$
(7)
where at time t, at resistant alleles from a read depth of At were observed. A low likelihood Lt indicates deviation from the neutral expectation; Bonferroni correction was used to identify significance at the 95% level. Data from each sample were analysed; four-fold non-neutral sites, being loci for which a significant likelihood was calculated from all four samples collected throughout each cross experiment, were identified.
Multi-locus model
A multi-locus model was developed to describe the manner in which allele frequencies would be expected to change over the course of the experiment; this model exploits information about the location of variant loci provided in the H. contortus reference genome, and the genome-wide map of recombination rate in the worms [21]; this map was inferred by characterising recombination breakpoints in F1 L3 progeny by whole genome sequencing, from which recombination rates were determined by comparing genetic distances between SNPs against their physical location in the genome. Where fixed genetic differences occur between the parental populations, the dynamics of adaptation in the resulting cross are relatively straightforward [48, 50]; this has been exploited to identify quantitative trait loci in yeast and malaria cross populations [49, 65]. Here a heuristic approach was used to identify fixed genetic differences between the parental populations before modelling evolution at these loci to identify the location of alleles conveying drug resistance.
Three filters were used to identify putatively fixed differences between parental populations. Firstly, we identified genomic sites for which the frequency of the resistant allele was 95% or greater in the resistant parent, requiring a read depth of at least 50x coverage for such sites to achieve a good level of statistical certainty. Secondly, we noted that, following the backcross performed at the start of the experiment, no locus can be homozygous for the resistant allele and as such, the resistant allele frequency can be no greater than 50% in the population; by way of reducing noise, we removed any locus with resistant alleles at a frequency of 60% or greater in the MHco3/10.BC4.noIVM sequence data. Thirdly, we noted that at homozygous separating sites in a population following a cross, the allele frequency will change smoothly over time; where there are N genomes in the population, the allele frequency, considered as a function of chromosome position, will change by 1/N wherever a crossover recombination event occurs within a genome. A diffusion model of allele frequency change, described in a previous publication [49], was used to identify allele frequencies across the genome that were consistent with this pattern. This analysis fits a posterior distribution to the allele frequencies across the genome, also inferring the extent of noise in the sequence data via a beta-binomial model intrinsic to the fitting process. Conservatively, sites no more than eight standard deviations from the mean of the posterior frequencies were included in the analysis. These approaches resulted in data from a total of 1368 loci in the MHco10(CAVR) dataset, and from 219 loci in the MHco4(WRS) dataset, were retained; the full set of allele frequencies are shown in Additional file 11: Figure S5.
Considering these data, we used an individual-based Wright-Fisher simulation to model the outcome of the experiment, taking into account the backcrossing, selection, and bottlenecking according to the experimental design. This approach accounted for the stochastic nature of allele frequency changes due to the effect of repeated population bottlenecking. We suppose that selection acts in favour of the resistance allele at a given locus when worms are in a sheep being treated with ivermectin, such that the fitness of individuals that are homozygous for the resistant allele (w1), heterozygous for the resistant allele (wh), or homozygous for susceptible alleles (w0) are given by:
$$ {\displaystyle \begin{array}{c}{w}_1=1+s\\ {}{w}_h=1+ hs\\ {}{w}_0=1\end{array}} $$
(8)
where s is the selection coefficient (fitness advantage of the homozygous resistant compared to homozygous susceptible) and h is the dominance coefficient (proportion of fitness change conveyed by a single copy of the allele). We assumed that all worms have equal fitness in the absence of the drug.
Our model is defined in terms of the manner in which selection acts upon the population, by the location in the genome of the selected allele, and the extent to which that allele conveys a fitness advantage to worms in the presence of drug treatment, this fitness advantage being characterised in terms of s and h. The likelihood of any given model was calculated using the beta-binomial likelihood function obtained from the diffusion process used above. Accounting for the stochasticity of the system, a set of at least 250 simulations were generated for each set of parameters, calculating the mean likelihood fit to the data across the simulations. A bootstrapping process was applied to quantify the uncertainty in this likelihood. Parameters giving the maximum likelihood fit to the data were identified.
Confidence intervals for the maximum likelihood locations of the selected allele were calculated. At the maximum likelihood position, we took the 250 simulations giving the maximum mean likelihood, calculating the mean of their likelihoods and the standard deviation of these values. We defined a threshold likelihood as the mean minus one standard deviation. For other allele positions, we also identified the parameters giving the maximum mean likelihood, constrained by allele position; we calculated a test likelihood equal to the mean plus one standard deviation of the replicate likelihoods generated by these parameters. Confidence intervals were then generated as the sites nearest the maximum likelihood position for which the test likelihood was lower than the threshold likelihood. We believe this gives a conservative estimate of the uncertainty in the location of the selected allele, given the stochasticity inherent in the experiment.
In an extension to this model, a two-driver scenario was considered, in which the resistant variant at each of two loci was under selection; to allow the exploration of model space in this case within reasonable computational time the assumption was made of additive selection at each locus. To distinguish this from the single-driver model, a requirement was imposed that the selected alleles be separated by at least 2 Mbp in the genome; this distance reflects the extent to which the process of recombination during the experiment allows the selective effects of distinct alleles to be observed.