Edinburgh Research Explorer Population genetics of benzimidazole-resistant Haemonchus contortus and Haemonchus placei from buffalo and cattle: implications for the emergence and spread of resistance mutations

The population genetics of nematode parasites are poorly understood with practical reference to the selection and spread of anthelmintic resistance mutations. Haemonchus species are important to study the nematode population genetics due to their clinical importance in ruminant livestock, and the availability of genomic resources. In the present study, it has been examined that Haemonchus contortus and Haemonchus placei populations from three buffalo and nine cattle hosts. Seventy-three individual adult worms of H. contortus and 148 of H. placei were analysed using a panel of seven microsatellite markers. The number of alleles per locus in H. contortus and H. placei indicated that all populations were polymorphic for the microsatellites used in the present study. Genetic diversity parameters included high levels of allelic richness and heterozygosity, indicating effective population sizes, high mutation rates and high transmission frequencies in the area. Genetic structure parameters revealed low genetic differentiation between and high levels of genetic variation within H. contortus and H. placei populations. Population dynamic analyses showed an absence of heterozygosity excess in both species, suggesting that there was no deviation from genetic drift equilibrium. Our results provide a proof of concept for better understanding of the consequences of specific control strategies, climatic change or management strategies on the population genetics of anthelmintic resistance alleles in Haemonchus spp. infecting co-managed buffalo and cattle.


Introduction
High levels of genetic diversity have been described in parasitic nematodes of the superfamily Trichostrongyloidea, attributed to large effective population sizes and high mutation rates (Blouin et al. 1995;Gilleard and Beech 2007;Prichard 2001). Studies in small ruminants have reaffirmed the contribution of high parasitic population size to genetic diversity in Haemonchus contortus (Chaudhry and Gilleard 2015;Prichard 2001;Redman et al. 2015). A number of studies have investigated how population genetic structure is partitioned between parasitic nematodes (Archie and Ezenwa 2011;Redman et al. 2008;Redman et al. 2015;Silvestre et al. 2009;Troell et al. 2003). A previous study of the global population genetic structure of H. contortus in small ruminants revealed a high level of genetic differentiation between continents, with populations from each essentially forming monophyletic groups (Troell et al. 2003). Similarly, high levels of genetic divergence were found between laboratory strains of H. contortus derived from different part of the world (Redman et al. 2008). These results are not surprising because of the limited potential for gene flow among geographically isolated populations. Within the countries, an early study of H. contortus field isolates in the USA suggested that there was essentially no genetic differentiation within North America due to very large effective population sizes combined with high gene flow by animal movement (Blouin et al. 1995). However, subsequent studies using more discriminatory markers suggest that relatively low but significant levels of regional genetic differentiation occur between H. contortus populations in Australia, UK, France and Sweden (Hunt et al. 2008;Redman et al. 2015;Silvestre et al. 2009;Troell et al. 2003). The results from the UK in particular show that genetic differentiation exists, even where high levels of animal movement are known to occur .
Few studies have investigated the population dynamics of parasitic nematodes due to their complex epidemiology involving free living or parasitic stages and effects of climate, animal management and host responses (Yin et al. 2016). For example, in case of H. contortus in small ruminants, the vast majority of environmental stages die during cold winters in Sweden and the UK, or dry and hot periods in Australia. Consequently, parasite survival at these particular times of year is confined to adult worms and inhibited early 4th stage larvae in their hosts. Subsequently, there is significant potential for population bottlenecks, which may allow sufficient genetic drift to account for genetic structure in these regions.
To date, most population genetic analyses have been performed in H. contortus from sheep and goats (Brasil et al. 2012;Chaudhry et al. 2015a;Chaudhry et al. 2016;Hunt et al. 2008;Hussain et al. 2014;Redman et al. 2015;Silvestre et al. 2009;Yin et al. 2016). Despite the economic importance of H. contortus and H. placei in buffalo and cattle in many subtropical regions, the only population genetic studies in these hosts are based on rDNA ITS-2 and ND4 mitochondrial markers (Ali et al. 2014;Brasil et al. 2012;Hussain et al. 2014). Although these markers are informative, they are not always sufficiently variable to show the population genetics of the parasites. Here for the first time, we use a panel of seven microsatellite markers to show the population genetic structure of H. contortus and H. placei from buffalo and cattle in six different locations of the Punjab province of Pakistan.

Field parasite samples
In the present study, several different regions were chosen from the Punjab province of Pakistan, where it was anticipated the prevalence of Haemonchus to be high. Adult Haemonchus spp. worms were obtained from the abomasa of three buffalo and nine cattle, immediately following slaughter at six different abattoirs (Lahore, Faisalabad, Sargodha, Sahiwal, Okara and Gujranwala). The numbers of H. contortus and H. placei collected from the three buffalo (Pop3B, Pop11B and Pop12B) and nine cattle (Pop1C,Pop2C,Pop4C,Pop5C,Pop6C,Pop7C,Pop8C,Pop9C,Pop10C) populations are described in Supplementary Table S1.

Genomic DNA extraction
Adult worms were fixed in 80% ethanol immediately following removal from the host abomasa. The heads of individual worms were dissected and lysed in single 0.5-uL tube containing 40 μl lysis buffer and proteinase K (10 mg/ml, New England BioLabs). One microliter of 1:5 dilution of neat single worm lysate was used as PCR template and identical dilutions of lysate buffer, made in parallel, were used as negative controls. Pooled lysates of each population were prepared using 1μl aliquots of each individual neat adult worm lysate. 1 μL of a 1:20 dilution of pooled lysates was used as PCR template (Chaudhry et al. 2016).
Pyrosequence specie-specific genotyping for the position 24 SNP of ITS-2 rDNA Genotyping of the SNP at position 24 (P24) of the rDNA ITS-2 region was used to confirm the identity of Haemonchus spp. in buffalo and cattle parasite populations as previously described by Chaudhry et al. (2015b). Briefly, the rDNA ITS-2 region was amplified from individual Haemonchus adult worm lysates using a Buniversal^forward primer complementary to 5.8S rDNA coding sequence and biotin labelled reverse primer complimentary to the 28S rDNA coding sequence. Final PCR conditions were 1× thermopol reaction buffer, 2 mM MgSO 4, 100 μM dNTPs, 0.1 μM forward and reverse primers and 1.25 U Taq DNA polymerase (New England Biolabs). Thermo-cycling parameters were 95°C for 5 min followed by 35 cycles of 95°C for 1 min, 57°C for 1 min and 72°C for 1 min with a single final extension cycle of 72°C for 5 min. Following PCR amplification of rDNA ITS-2, the single nucleotide polymorphism (SNP) at P24 was determined by pyrosequence genotyping using the PryoMark ID system (Biotage, Sweden). The sequencing primer used was Hsq24 (5′-CATATACTACAATGTGGCTA-3′) and the nucleotide dispensation order was CGAGTCACA. Peak heights were measured using the SNP mode in the PSQ 96 single nucleotide position software (Biotage, Sweden). Worms were designated as H. contortus, H. placei or putative hybrids based on being homozygous A, homozygous G or heterozygous A/G, respectively at P24 (Chaudhry et al. 2015a).
Microsatellite genotyping of H. contortus and H. placei populations Seven previously published microsatellite markers (Hcms3561, Hcms53265, Hcms36, Hpms43, Hpms52, Hpms53, Hpms102) were selected based on known polymorphism in H. contortus and H. placei worms (Chaudhry et al. 2015a;Chaudhry et al. 2016;Santos et al. 2017). A summary of primers sequences and allele ranges are given in Supplementary Table S1. PCR amplification was performed using a 25μl master mix containing final concentrations of 1X thermopol reaction buffer, 2 mM MgSO 4, 100uM of each dNTP, 0.1μM forward and reverse primers and 1.25 U Taq DNA polymerase (New England Biolabs). Thermo-cycling parameters were 94°C for 2 min followed by 40 cycles of 94°C for 15 s, 54°C for 30 s and 72°C for 1 min with a single final extension cycle of 72°C for 15 min. The forward primer of each microsatellite primer pair was 5′ end labelled with fluorescent dye (IDT,Canada) and the GeneScan ROX 400 internal size standard was used on the ABI Prism 3100 genetic analyser (Applied Biosystems, USA). Individual chromatograms were analysed using Gene Mapper software version 4.0 for the accurate size of the amplicons and determine genotypes (Applied Biosystems, USA) previously described by Chaudhry et al. (2016).
The multi-locus microsatellite genotype data analysis was previously described by Chaudhry and Gilleard (2015). A log likelihood ratio test statistic (G test) was used to estimate the linkage equilibrium, using Arlequin 3.11 (Excoffier et al. 2005). Expected and observed heterozygosities (H e and H o ) and allele variation (A C ) for each locus were calculated using Arlequin 3.11. An exact test was used to statistically evaluate deviations from Hardy-Weinberg equilibrium for all populations (Guo and Thompson 1992). Significance levels were adjusted using the sequential method of Bonferroni correction for multiple comparisons in the same dataset (Rice 1989). Analysis of molecular variance (AMOVA) was estimated through partition of genetic variation between and within populations (Excoffier et al. 1992). Fixation index (pairwise F ST ) values were calculated from the multi-locus microsatellite genotype data, by random permutation in Arlequin 3.11. Principal coordinate analysis (PCoA) was performed using GenALEx software preserving individual worm genotypes to plot coordinates of individuals (Peakall and Smouse 2012). The bottleneck software was used to assess any possible recent reduction in effective population size (Luikart and Cornuet 1998).

H. contortus and H. placei co-infections are common in buffalo and cattle
Up to 20 individual worms from each host were pyrosequence genotyped for the rDNA ITS-2 P24 SNP (228 worms genotyped in total) (Supplementary Table S1). All worms were identified as H. contortus (homozygous A at rDNA ITS-2 P24) in population number Pop1C of cattle and Pop3B of buffalo and H. placei (homozygous G at rDNA ITS-2 P24) in population number Pop4C, Pop5C, Pop6C and Pop7C of cattle and Pop11B of buffalo. The remainder five populations (Pop2C, Pop8C, Pop9C, Pop10C and Pop12B) of the buffalo and cattle examined contained a mixture of H. contortus (homozygous A at rDNA ITS-2 P24) and H. placei (homozygous G at rDNA ITS-2 P24) indicating co-infection with the two species (Supplementary Table S1). Three of the individual hosts (two cattle and one buffalo) also each contained a single worm with a heterozygous A/G genotype at the rDNA ITS-2 P24 position, suggesting that they were H. contortus/H. placei hybrids.  (Table 1). Departure from Hardy-Weinberg equilibrium was used to assess genetic variability. There was some significant departure from Hardy-Weinberg equilibrium, even after Bonferroni correction, in addition to relative P values for 3 out of 35 loci combinations for H. contortus and 7 out of the 63 loci combination for H. placei populations (Table 1). In addition, H o was less than the H e and inbreeding coefficient (F IS ) values were high ( Table 1). The presence of null alleles (N 0 ) for microsatellite loci has been previously reported and is likely associated for these departures from Hardy-Weinberg equilibrium high level in both species (Chaudhry et al. 2015a;Chaudhry et al. 2016). There was no evidence to support linkage disequilibrium for any combination of loci across all populations, indicating that alleles at these loci were randomly associating and not genetically linked. Overall results indicated high level of genetic diversity was attributed to a large population     were not geographically sub-structured (Fig. 1).

Population dynamics of H. controtus and H. placei
The population bottlenecking analysis of both H. contortus and H. placei showed that there was no heterozygosity excess according to the Sign Test and Wilcoxon Test (Table 3). The mode shift analysis established that all populations had a normal L-shaped distribution. The overall results suggest that population studied did not appear to deviate from genetic drift equilibrium.

Discussion
Anthelmintic resistance in Haemonchus spp. of large ruminants now represents a serious challenge to the livestock industry worldwide (Kaplan and Vidyashankar 2012). Hence, it is important to understand the population genetics of H. contortus and H. placei of large ruminants and their implications for the emergence and the spread of resistance mutations. The data reported in the present study reveal a high level of genetic diversity in H. contortus and H. placei from buffalo and cattle hosts in Pakistan. Similar results were seen in previous Indian and Pakistani studies of sheep and goats (Chaudhry et al. 2015a;Chaudhry et al. 2016). Our results differ from those of Blouin et al. (1995), Brasil et al. (2012) and Jacquiet et al. (1995), which showed higher genetic diversity in H. contortus than in H. placei in North America, Brazil, and Mauritania. This may be due to different evolutionary rates or population demography (Brasil et al. 2012). The implication of the high level of genetic diversity on the emergence of benzimidazole resistance has been described by number of population genetic studies. The multiple time emergence of the F200Y (TAC) and F167Y (TAC) SNPs at the isotype 1 ß tubulin locus (Kwa et al. 1994) through either recurrent, or pre-existing mutations; and the single time emergence of the E198A (GCA) SNP at the same locus through a single mutation have been reported, associated with high levels of genetic diversity in H. contortus populations of small ruminants (Chaudhry et al. 2015a).
In the present study, genetic differentiation occurred in H. contortus and H. placei from Pakistani buffalo and cattle at low but significant levels. Similar findings have been described in Australia and the UK in H. contortus of small ruminants (Hunt et al. 2008;Redman et al. 2015). The consequences of low levels of population genetic structure in H. contortus and H. placei within the region, at least in part reflects high levels of animal movement (Blouin et al. 1995;Chaudhry et al. 2015a;Chaudhry et al. 2016;Hunt et al. 2008;Redman et al. 2015). Conversely, if the gene flow is high, the anthelmintic resistance mutations potentially spread within regions. Low levels of genetic variation within H. contortus and H. placei populations also imply that there is no reproductive isolation. Interestingly, the Pop9C population of H. contortus from Okara showed a high F st value, which might be related to limited trade, brought about by known lower levels of communication and economic isolation, hence lower levels of animal movement compared to the other regions. Moreover, confounding environmental management, animal breed, host density and climatic conditions may also influence the epidemiological pattern of this population of H. contortus.
The free-living stages of parasitic nematodes are influenced by environmental factors, in particular temperature and humidity (O'Connor et al. 2006), which vary depending on geographical location and season, while the parasitic stages are influenced by host responses. The success of parasitic nematodes is due to their adaptation to these factors . Haemonchus spp. are adapted to warm temperatures and plentiful precipitation, while dry, hot and cold weather is usually fatal for the free-living larval stages (Besier et al. 2016). Punjab is a subtropical region of Pakistan, where climatic conditions are warm with moderate precipitation, throughout the year, while winter temperatures are mild (between 10 and 20°C) providing consistent opportunities for genetic exchange, transmission and dispersal in both H. contortus and H. placei. Our findings suggest that these conditions may have resulted in little population bottlenecking or genetic drift, potentially even for those parasite populations under drug selection and high gene flow.

Conclusions
This is the first report of the population genetic study of H. contortus and H. placei from buffalo and cattle hosts using a panel of microsatellite markers. The results of the low level of genetic differentiation among H. contortus and H. placei populations of Punjab, Pakistan, may be explained by large effective population sizes and low rates of genetic drift. Furthermore, limited animal movement associated with trade may have led to the reduction of gene flow and isolation of H. contortus and H. placei populations in this region. Pakistan for his great support to collect samples from abattoirs. We would also like to thank Moredun Research Institute Scotland for their kind support to use pyro-sequencer. The article is the sequel of Dr. Umer Chaudhry thesis, which is available in full at the University of Calgary, Canada website.

Funding information
The study was financially supported by the Higher Education Commission of Pakistan.

Compliance with ethical standards
Conflict of interest The authors declare that they have no competing interests.
Open Access This article is distributed under the terms of the Creative Comm ons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.