Population structure of Haemonchus contortus from seven geographical regions in China, determined on the basis of microsatellite markers

Background Studying genetic variation within and among Haemonchus contortus populations can inform some aspects of this parasite’s population genetics and epidemiology. However, almost nothing is known about such variation in China. Methods Adult males of H. contortus (n = 184) representing seven distinct populations in China were collected, and genetic variation within and among these populations was explored using eight distinct microsatellite markers. Results Genetic parameters, such as heterozygosity and inbreeding coefficient (F IS) indicated that all eight microsatellites were highly polymorphic. Various analyses (AMOVA, F ST, phylogenetic, structure, mantel test and population dynamics) revealed high within-population variation, low population genetic differentiation and high gene flow for H. contortus in China. Conclusions This study provides a first snapshot of the genetic substructuring of H. contortus populations in China using polymorphic markers, and might provide a starting point for assessing genetic changes over space and time during or following the implementation of particular treatment or control strategies, or changes as a consequence of environmental, management and climatic factors. Electronic supplementary material The online version of this article (doi:10.1186/s13071-016-1864-z) contains supplementary material, which is available to authorized users.

Although microsatellite analyses have been applied specifically to H. contortus populations in Australia, some European countries (i.e. France, the Netherlands, Sweden and the UK) and, more recently, India [17][18][19][20], this is not the case for China. This is surprising, given that this nematode is likely the most economically important parasite of livestock in China [12]. Indeed, there is no detailed information on the population genetic structure(s) of H. contortus populations in this country using polymorphic markers. Therefore, in the present study, we investigated the structures of seven H. contortus populations from disparate geographical regions across China, including Inner Mongolia, using eight microsatellite markers.

H. contortus populations
A total of 184 individual adult males of H. contortus were collected from the abomasa of sheep or goats from seven geographical locations (23-29 specimens per location/ population) in Heilongjiang and Liaoning in northeastern China; Shaanxi and Inner Mongolia in northwestern China; Hubei in central China; and Yunnan and Guangxi in southwestern China (Table 1; Additional file 1: Figure  S1). These locations are in temperate to subtropical climate zones and are separated by distances of 800 to 5000 km (Additional file 1: Figure S1). Each individual adult male worm was identified morphologically to genus level [23], and specific identification was based on the sequence of the second internal transcribed spacer of nuclear ribosomal DNA (ITS-2) ( [24]; GenBank accession no. X78803) (see below). Samples were stored in 70% ethanol and frozen at -20°C until use.

Isolation of genomic DNA
Genomic DNA was extracted from single worms using sodium dodecyl-sulfate/proteinase K treatment, followed by spin-column purification (Wizard DNA Clean-Up, Promega, Madison, USA) [25] and then stored at -20°C until use.

PCR amplification of ITS-2 and sequencing
To confirm the identity of the nematodes as H. contortus, the ITS-2 of each individual worm was amplified using the primers NC1 (5′-ACG TCT GGT TCA GGG TTG TT-3′) and NC2 (5′-TTA GTT TCT TTT CCT CCG CT-3′) [24]. Briefly, PCR was performed in 25 μl containing 10 mM Tris-HCl (pH 8.3), 50 mM KCl, 4 mM MgCl 2 ,250 μM of each dNTP, 100 pmol of each primer, 1 U Taq DNA polymerase [TaKaRa, Dalian, China]) and 30-50 ng of DNA template (except for notemplate controls). The cycling protocol was: initial denaturation at 94°C for 5 min, followed by 30 cycles of denaturation at 94°C for 30 s, annealing at 55°C for 30 s and extension at 72°C for 1 min, with a final extension of 72°C for 5 min. Amplicons were examined on agarose gels (1.5%) to verify that they represented single bands and then column-purified (PCR-Preps, Promega) and sequenced directly (using primer NC2) in an automated DNA sequencer (ABI 3730).

PCR amplification of microsatellite loci
Eight microsatellite markers (designated Hcms 15,19,22co3,28,33,36,37 and 40) were selected based on previous publications [17,18]. PCR was performed in 20 μl containing 10 mM Tris-HCl (pH 8.3), 50 mM KCl, 4 mM MgCl 2 , 250 μM of each dNTP, 100 pmol of each primer, 1 U of Taq DNA polymerase (TaKaRa) and 30-50 ng of DNA template (except for no-template controls). The cycling protocol was: initial denaturation at 94°C for 5 min, followed by 30 cycles of denaturation at 94°C for 30 s, annealing at 50-60°C (temperatures for individual primer sets, see refs. [17,18]) for 30 s and extension at 72°C for 1 min, with a final extension of 72°C for 5 min. Each forward primer was labelled at their 5-ends with the fluorescent dye FAM (Hcms 22co3, 28, 33 and 37) or HEX (Hcms 15, 19, 36 and 40). The sizes of individual amplicons were established by capillary electrophoresis in an ABI Prism 3730 XL analyzer. Individual chromatograms were assessed using GeneMarker 1.51 software (Applied Biosystems, USA) to assign allele fragment lengths to individual samples. The data set was assessed for errors using the Micro-Checker program [26], and any errors were corrected.

Data analyses
Genetic characteristics were estimated by calculating observed and expected heterozygosities (H o and H e ), the number of alleles (A) and population inbreeding coefficients (F IS ) for individual loci using the program GenA-LEx [27]. Exact tests for deviation from Hardy-Weinberg equilibrium and pairwise linkage disequilibrium were conducted using Genepop [28]. To correct for multiple tests, Bonferroni corrections were employed to adjust significance values [29].
Pairwise F ST values were calculated using the program Arlequin 3.5 [30] to evaluate the genetic differentiation among populations. Analysis of molecular variance (AMOVA) was used to estimate genetic diversity within and among populations using the same program (Arlequin 3.5). Genetic structure was inferred using Bayesian Markov Chain Monte Carlo (MCMC) approach, and the Bayesian-based analysis was run in the program STRUCTURE [31]. Principal coordinate analysis (PCoA) was performed using GenALEx software, preserving individual worm genotypes to plot individuals [27]. To determine the correlation of geographical distance and genetic distance, the Mantel test was conducted using the TFPGA package [32]. The Bottleneck program [33] was used to assess any possible, recent reduction in effective population size.

Genetic diversity and Hardy-Weinberg equilibrium
All 184 adult specimens of Haemonchus were identified individually as H. contortus based on their ITS-2 sequence. Individual adult male worms from each of the seven populations were then genetically characterized using the eight microsatellite markers (Hcms 15, 19, 22co3, 28, 33, 36, 37 and 40). All of these markers displayed a high degree of polymorphism, with the number of 'alleles' per locus ranging from 6 to 15 (Table 1). The observed heterozygosity (H o ) ranged from 0.103 to 0.923, and the expected heterozygosity (H e ) from 0.237 to 0.848 (Table 1). For each population, the average number of alleles at all loci ranged from 4.750 to 5.375, and the average, expected heterozygosity for all loci was > 0.5 for all populations (Table 1). Among all seven populations of H. contortus, those from Shaanxi and Heilongjiang were most polymorphic, with an average number of 5.375 alleles, whereas the population from Liaoning had the smallest mean number of alleles (4.750).
Population genetic analysis on departure from Hardy-Weinberg equilibrium (HWE) was undertaken to assess the suitability of the eight selected markers to assess genetic variability in the seven H. contortus populations. For some marker/population combinations, there was no significant departure from HWE after sequential Bonferroni correction, and F IS values were low (see Table 1). However, for each of five markers (Hcms15, Hcms40, Hcms28, Hcms22co3 and Hcms37), there was a significant departure from HWE in one to five populations with heterozygotic deficiency or excess (see Table 1). In the populations with heterozygosity deficiency, the observed heterozygosity (H o ) was less than the expected heterozygosity (H e ), and F IS values were high (see Table 1), suggesting the presence of null alleles at four loci (Hcms15, Hcms40, Hcms28 and Hcms22co3). There was also evidence that null alleles were observed for the same four loci (assessed using Micro-Checker software; data not shown). Although there were deviations from HWE for five markers in some H. contortus populations (Hcms15/Heilongjiang; Hcms40/Yunnan/Liaoning/Heilongjiang; Hcms28/Guangxi; Hcms22co3/Hubei/Yunnan/Guangxi/Liaoning/Inner Mongolia; Hcms37/Hubei/Guangxi/Shaanxi; Table 1), there was no evidence to support linkage disequilibrium for any combination of loci for individual populations, suggesting that these loci were not genetically linked and could be used independently.

Genetic structure
Pairwise F ST values were calculated to estimate the degree of genetic differentiation among the seven H. contortus populations using the eight microsatellite markers ( Table 2). Although most pairwise F ST values were low (< 0.05), ranging from 0.0024 to 0.0640, the population from Inner Mongolia showed the highest level of genetic divergence from other populations, with F ST values ranging from 0.0352 to 0.0640 (Table 2). Low F ST values indicated limited differentiation and high gene flow among most of the populations.
Bayesian clustering analysis showed no obvious subdivision at the population level (results not shown). We explored K values between 1 and 10. For K = 3, corresponding lnP(D) was greatest, and the peak of ΔK was highest (see Fig. 1). To determine whether the H. contortus populations comprised a single panmictic population with a high degree of gene flow, an analysis of molecular variance (AMOVA) among populations was conducted (Table 3). An analysis using the STRUCTURE program supported the division of H. contortus populations into three groups: group 1 (Guangxi, Hubei and Yunnan), group 2 (Heilongjiang, Liaoning and Shaanxi) and group 3 (Inner Mongolia). Most (96.1%) of the genetic variation was distributed within populations: 3.3% of variance components among groups (F CT ) and 0.6% of variance components among populations within groups (F SC ). These findings suggest a very low genetic substructuring of H. contortus populations in China.

Geographical sub-structuring
Evidence of sub-structuring was assessed by Principal coordinate analysis (PCoA) analysis (Fig. 2). The two axes (i.e. 21.2 + 24.5%) accounted for 45.7% of the Nei's genetic distance [34] was calculated to evaluate whether there was an effect of geographical separation on population differentiation. The analysis of the relationship between genetic distance and geographical distance for the seven H. contortus populations indicated that the genetic differentiation of these populations did not follow a pattern of isolation by distance (Mantel test: r = 0.1654, P = 0.284) (Fig. 3), indicating that genetic differentiation was not correlated with geographical distance for all populations, as a consequence of high gene flow among populations.

Population dynamics
The analysis of the seven H. contortus populations showed that the populations from Inner Mongolia and Yunnan exhibited heterozygosity-excess, according to the IAM and TPM/SMM models of Sign test and the Wilcoxon signed-rank test ( Table 4). The Mode-shift analysis demonstrated that the allele frequencies had a normal Lshaped distribution. Although the heterozygosity-excess was detected in two populations (i.e. Inner Mongolia and Yunnan), the results indicated that the populations studied did not appear to deviate from a mutation-drift equilibrium (Table 4).

Discussion
In the present study, we explored the population genetic structure of seven H. contortus isolates from China using a panel of eight polymorphic microsatellite markers. The allelic richness and overall heterozygosity of these microsatellites were shown to be high, supporting previous findings using mtDNA (nad4) nucleotide sequence data [12] and revealing high levels of genetic diversity within and among H. contortus populations in China.
Deviation from HWE was observed for five of eight markers for some populations (Hcms15/Heilongjiang; Hcms40/Yunnan/Liaoning/Heilongjiang; Hcms28/ Guangxi; Hcms22co3/Hubei/Yunnan/Guangxi/Liaoning/ Heilongjiang; Hcms37/Hubei/Guangxi/Shaanxi) from China, following Bonferroni correction, which is likely due to the presence of null alleles. As the repeated "genotyping" from the same samples gave consistent results, and individual samples that failed to amplify for one locus could be genotyped robustly for other loci, the null alleles might not relate to allelic drop-out or poor template quality; otherwise, the apparent 'null alleles' might relate to sequence variation at primer annealing sites and/or due to selection, mutation or genetic drift [7]. The high frequency of null alleles in Chinese H. contortus populations is in accordance with the previous studies of H. contortus [18] and other trichostrongylid nematodes including Teladorsagia circumcincta and T. tenuis [6,7,35].  No linkage disequilibrium was detected across the entire population. It is proposed that linkage disequilibrium results from a low level of genetic exchange, population sub-structuring, recombination and/or inbreeding [36]. The absence of linkage disequilibrium from Chinese populations of H. contortus was supported by AMOVA analysis results, which showed that the majority (96.1%) of genetic variation was within H. contortus populations, suggesting a high level of gene flow among populations, and PCoA results indicated that no discernable geographical sub-structuring was evident over large distances in China (Fig. 2). Therefore, both AMOVA and PCoA results demonstrated high gene flow and low geographical sub-structuring among the seven Chinese populations of H. contortus. The high rate of gene flow appears to relate to the frequent trade and transport of the sheep and goats with H. contortus infection across vast distances in China (cf. [12,37]).
Although the numbers of samples studied here were relatively small, genetic differentiation between H. contortus populations from sheep and goats were not apparent, showing that there is no reproductive isolation [13]. Compared with other H. contortus populations in China, however, the population from Inner Mongolia showed a higher F ST value (Table 2), which might relate to factors including unique local agricultural management and cultural factors in the Inner Mongolia autonomous region, the second largest plateau in China. The nomadic life  In subtropical areas (Guangxi, Hubei and Yunnan) of China, the main climatic characteristics are warm temperatures and high rainfall in summer. Otherwise, winters are cold (-20 to -30°C), and summers are hot (20 to 30°C) and dry, with limited rainfall throughout much of the year in temperate zones (Heilongjiang, Liaoning and Inner Mongolia). Previous studies [38,39] have shown that temperature and moisture have a dominant effect on the development and survival of the free-living stages of H. contortus; hence, the availability of infective larvae and rates of infection are affected. For this parasite, abundant rainfall and warm temperatures are favourable for development, but hot, dry and extreme cold conditions are usually lethal to infective larvae on pastures [40]. Consequently, the larval stages of H. contortus are often poorly adapted to cooler and drier climates, such that the perpetuation of the life-cycle is largely linked to the survival of parasitic stages in host animals to carry infection through from year to year [41]. The Autonomous Region of Inner Mongolia has a monsoon climate at medium latitudes, but temperature and precipitation appear to have changed in recent years, leading to a drier climate [42]. The ecological system in Inner Mongolia might limit the opportunity for the transmission as well as recombination, dispersal and genetic exchange of H. contortus with other populations in China. Thus, the low level of genetic exchange and higher genetic diversity in H. contortus might be impacted by environmental and/or microclimatic factors in Inner Mongolia.
The results of the low level of genetic differentiation among the seven H. contortus populations in China by microsatellite analysis were essentially in agreement with those of a previous study based on mitochondrial nad4 gene sequence data [12]. These findings are also supported by the proposal that limited genetic structuring of H. contortus populations might be ascribed to a large effective population size and low rates of genetic drift [2,4]. However, it is noteworthy that there appears to be lower level of gene flow and population sub-structuring among H. contortus populations on a global scale, supported by the studies employing both microsatellite markers and the nad4 gene sequence data [9,18]. Apparently, the restricted movement of livestock associated with the trade barriers and agricultural policies can lead to reduced gene flow and the formation of 'isolated' H. contortus populations.

Conclusions
This is the first population genetic study of H. contortus in China using microsatellite markers. The findings provide a snapshot of the genetic make-up of H. contortus populations in parts of China using these polymorphic markers, revealing high within-population variation, low population genetic differentiation and high gene flow for H. contortus. This study should provide a reference for assessing genetic change in H. contortus over space and time during or following the implementation of particular treatment or control strategies, or as a consequence of environmental, management and climatic factors.

Additional file
Additional file 1: Figure S1.