Genetic diversity and population genetics of the warble flies Hypoderma bovis and H. sinense in Qinghai Province, China

Background Hypoderma bovis and H. sinense (Diptera: Oestridae) mainly parasitise cattle and yaks. The two parasites are pathogenic and cause economic losses that result from reduced amounts of livestock products, including milk, meat, and skin. Genetic diversity and population genetic structure of H. bovis and H. sinense have not been evaluated, but could be used to inform appropriate strategies to control these parasites. Methods We cloned and sequenced part of the mitochondrial cytochrome c oxidase subunit I (COI) gene from 60 H. bovis isolates and 52 H. sinense isolates from five locations in Qinghai Province, China, to identify polymorphisms, and infer their phylogenetic relationships, historical population expansions, and divergence time. Results We identified 17 COI haplotypes from the H. bovis samples, and 23 COI haplotypes from the H. sinense samples. The haplotype and nucleotide diversities were 0.738 and 0.00202 for H. bovis, and 0.867 and 0.00300 for H. sinense, respectively, which indicates rich genetic diversity in H. bovis and H. sinense populations. Bayesian phylogenetic analysis revealed that the two species are monophyletic, and geographical structuring of haplotypes was significantly different in H. sinense (P < 0.05), but not H. bovis. Neutrality tests and mismatch distribution statistical analysis revealed that populations of the two species have undergone demographic expansions. The divergence three Hypoderma spp. (H. bovis, H. lineatum, and H. sinense) was estimated to have occurred approximately 4.5 million years ago (Mya), which indicates that the rapid uplift of the Qinghai-Tibetan Plateau during the late Miocene-Pliocene was associated with divergence of Hypoderma species. Conclusions Results of the present study revealed that both H. bovis and H. sinense displayed high genetic diversity and widespread population genetic differentiation within and among populations; these data, along with the molecular phylogeny, demographic history, and divergence time estimation, provide new insight into evolutionary history of these species. These findings will help elucidate speciation in Hypoderma and provide theoretical basis for epidemiological surveillance and control of these species on the Qinghai-Tibetan Plateau.


Background
Hypoderma bovis and H. sinense (Diptera: Oestridae) are two species of flies in Oestridae and mainly parasitize cattle and yaks. The parasitizing flies are widely distributed in north and southwestern China [1,2]. The prevalence of Hypoderma spp. larval infection in yaks can reach up to 100 % in some areas of Qinghai Province [3]. Hypodermosis of cattle and yaks, caused by the larvae of Hypoderma spp. , is responsible for substantial economic losses in the livestock industry because it results in spontaneous abortion, reduced milk production, loss of weight, reduced fertility, and poor hide quality [4,5]. Therefore, there is a need to develop effective strategies to control this disease.
Therefore, in the present study, we analysed the genetic diversity of H. bovis and H. sinense based on mitochondrial COI sequences in samples collected from Qinghai Province, China. In addition, we investigated possible historical population expansions and divergence time of H. bovis and H. sinense. These findings are essential for understanding speciation of Hypoderma spp. and for epidemiological surveillance and control of these species on the Qinghai-Tibetan Plateau (QTP).

Sampling strategy
We collected 60 third-stage H. bovis larvae and 52 thirdstage H. sinense larvae from five localities in Qinghai Province from 2013 to 2014 (Fig. 1). The initial identification of H. bovis and H. sinense was mainly based on morphological characteristics [17], and confirmed by molecular methods using the mitochondrial COI gene [8]. All specimens were fixed by immersion in 70 % ethanol. The locations and sample numbers of H. bovis and H. sinense populations are shown in Table 1.

DNA extraction, amplification, cloning, and sequencing
The third-stage fly larvae were longitudinally cut to retrieve the internal organs. The genomic DNA was extracted from 10 mg of each internal organ using a commercial kit (TIANamp Genomic DNA Kit, TIANGEN Biotechnology, Beijing, China) in accordance with the manufacturer's recommendations.

Population haplotype diversity analysis
COI sequences were aligned using MEGA 5.2 [18]. Identical haplotypes were collapsed using DNASP 5.10 [19]. The number of haplotypes and standard diversity indices [haplotype and nucleotide diversities (h and π, respectively)] were calculated using DNASP 5.10 [19] for each population.

Phylogenetic analysis and haplotype network construction
Phylogenetic relationships of H. bovis and H. sinense COI haplotypes were inferred using Bayesian inference (BI). We selected the best-fit model (GTR + I + G) for BI analyses for each data partition using Modeltest 3.7 [20] in conjunction with PAUP 4.0b10 [21]. A Bayesian tree was constructed using MrBayes 3.1.2 [22], and Markov chain Monte Carlo was run for 10 million generations with sampling every 1000 generations. The first 25 % of generations were discarded as burn-in, and the remaining trees were used to estimate Bayesian posterior probabilities (PP).
COI sequences of H. bovis (AF497761) and H. sinense (AY350769) obtained from the GenBank database were used for phylogenetic analysis of the species in this study, and COI sequences from three other species of Hypoderma, H. lineatum (AF497762), H. tarandi (AF497764) and H. actaeon (AF497765), were used as in-group taxa for the phylogenetic analysis. Gasterophilus pecorum (AF497776) was selected as the out-group taxon to root the phylogenetic trees. Median-joining networks of all H. bovis and H. sinense haplotypes in this study were constructed using Network 4.6 [23] to visualize relationships among unique haplotypes.

Population genetic and demographic history analyses
Analysis of molecular variance (AMOVA) was used to evaluate H. bovis and H. sinense population genetic structure in Arlequin 3.11 with 1000 permutations [24]. This study implemented two levels of AMOVA for intraand inter-population analyses (Φ ST ). Phylogeographic structure of H. bovis and H. sinense populations was  (3) N number of individuals sequenced, NH number of different haplotype, h haplotype diversity, π nucleotide diversity; The number of individuals observed for each haplotype is given in parentheses estimated using PERMUT (http://www.pierroton.inra.fr/ genetics/labo/Software/Permut) with 1000 permutations. PERMUT tested phylogenetic structure by calculating G ST (only haplotype frequencies) and N ST (differences among haplotypes); phylogenetic structure is usually considered present when N ST is higher than G ST [25]. Fu's Fs and Tajima's D neutrality tests, and sum of squared deviation (SSD) and Harpending's raggedness (rg) test statistics of mismatch distributions were calculated to detect evidence of past population expansions in Arlequin 3.11 [26][27][28].
In addition, mismatch distributions of H. bovis and H. sinense populations were performed using DNASP 5.10 [19] to test whether demographic processes were consistent with the mismatch distribution test statistics. A population usually exhibits a uni-modal mismatch distribution when it has passed through a recent demographic expansion [29], whereas a multimodal mismatch distribution indicates that a population is comparatively stable [30]. Expansion time was estimated using the expectation τ = 2 ut [31], where τ represents the mode of the mismatch distribution, t represents time in generations since expansion, and u = 2 μk, where μ is the mutation rate (COI was estimated to be approximately 2 % per million years for Hypoderma spp. ) [7] and k is the length of the sequence [32]. The H. bovis and H. sinense generation time was estimated to be 1 y based on their life-cycle [33].

Divergence time estimates
Divergence times were estimated using the Bayesian phylogenetic method implemented in BEAST 1.8.0 [34]. The clock model was set to relaxed, uncorrelated lognormal [35], with a Yule speciation tree model. Because there is a lack of Hypoderma fossils, a secondary calibration approach was used. Based on previous research on divergence time in Hypoderma spp. [10], a mean of normal distribution with standard deviation were set as 8.2 million years ago (Mya) with 0.5 Mya for H. tarandi and H. actaeon, and 4.2 Mya with 0.5 Mya for H. bovis and H. lineatum, respectively. The Markov chain Monte Carlo chain length was set to 10 million generations and sampled every 1000 generations. Chain convergence was assessed to determine effective sample sizes greater than 200 for all parameters using Tracer 1.5 [36], and trees were summarized using TreeAnnotator 1.8.1 [34].

Sequence variation and haplotype diversity
There were no insertions or deletions of nucleotides in the 689-bp COI sequences amplified from any H. bovis and H. sinense. A total of 19 nucleotide polymorphisms (12 singleton variable sites and seven parsimony-informative sites) for H. bovis and 31 nucleotide polymorphisms (24 singleton variable sites and seven parsimony-informative sites) for H. sinense were detected. In total, 17 haplotypes were detected in H. bovis populations and 23 haplotypes were detected in H. sinense populations ( Table 1). Sequences of all haplotypes have been deposited in the GenBank under accession numbers KT600277-KT600316. The ratio of haplotypes relative to the total number of individuals sampled for each species was 0.28 for H. bovis and 0.44 for H. sinense. There was no significant difference in the number of haplotypes found in each sampling location between the two species (F 1, 8 = 0.459, P = 0.517; Table 1). Values of h and π were 0.738 and 0.00202 for H. bovis, and 0.867 and 0.00300 for H. sinense, respectively, which indicates rich genetic diversity in H. bovis and H. sinense populations; however, there were no significant differences between the two species (F 1, 8 = 0.897, P = 0.371; F 1, 8 = 0.809, P = 0.395; Table 1).
The network for clade HB showed that Haplotype HB3 was considered the central haplotype, to which a large number of private haplotypes (76.5 % of the total haplotypes of clade HB) were connected in a star-like manner (Fig. 2). The highest-frequency haplotype was HB3, followed by HB7, HB13, and HB8, which included 30, 5, 5, and 3 individuals, respectively. For clade HS, Haplotype HS6 was considered the central haplotype. The private haplotypes represented 87.0 % of all clade HS haplotypes. The highest-frequency haplotype was HS6, followed by HS1 and HS3, which included 16, 10, and 5 individuals, respectively, and they occupied a central position in the network (Fig. 2).

Genetic differentiation and population structure
AMOVA results showed that there was significant genetic differentiation in H. bovis and H. sinense populations ( Table 2). For H. bovis, AMOVA showed that 11.88 % of the variation was among populations and 88.12 % was within populations. High genetic structure was found (Φ ST = 0.119, P < 0.001), which indicates remarkable genetic differentiation in H. bovis. For H. sinense, AMOVA showed that 7.83 % of the variation was among populations and 92.17 % was within populations. We also detected high genetic structure (Φ ST = 0.078, P < 0.001), which likewise indicates significant genetic differentiation, in the H. sinense isolates. Large pairwise F ST values were found between H. bovis and H. sinense populations. For H. bovis, pairwise F ST ranged from 0.009 to 0.190, and most pairwise values were statistically significant (Table 3). For H. sinense, pairwise F ST ranged from − 0.037 to 0.210, and most pairwise values were also statistically significant (Table 3).
Demographic expansions were analysed for H. bovis and H. sinense populations using two neutrality tests and mismatch distributions. For H. bovis and H. sinense, Fu's Fs and Tajima's D values were significantly negative (Table 4), and mismatch distributions of both species each showed a unimodal curve (Fig. 3). From the τ value (2.652) calculated by Arlequin, expansion of H. bovis populations was estimated to have occurred about 0.049 Mya; expansion of H. sinense populations (τ = 3.098) occurred about 0.056 Mya. A permutation test showed that N ST (0.141) was significantly higher than G ST (0.048) for H. sinense (P < 0.05), whereas N ST (0.116) was less than G ST (0.122) for H. bovis, which indicates that significant phylogeographic structure is apparent in H. sinense, but not H. bovis.

Divergence times
The estimated evolutionary timescale of five Hypoderma species with the 95 % highest posterior densities (95 % HPD) intervals are presented in Fig. 4. Our analysis estimated that the most recent common ancestor of the five Hypoderma species existed approximately 8.1 Mya (95 %  Table 1

Discussion
The haplotype (h) and nucleotide (π) diversity are important indicators of genetic diversity in a population [19]. The results showed high h and π in H. bovis and H. sinense populations, which indicates rich genetic diversity of populations of the two species and might explain why the two species have a broad tolerance to environmental and habitat stresses; fast mutational processes inherent in individuals and populations may enable these two Hypoderma spp. to successfully adapt to complex environments. The Bayesian phylogenetic analysis strongly supported the coalescence of COI haplotypes within H. bovis and H. sinense (PP = 0.99; Fig. 2), which supports the notion that these are different species and is consistent with of the findings of a previous molecular study [8]. In this tree, H. sinense and H. lineatum (PP = 0.76; Fig. 2) were more closely related to each other than to H. bovis. This result was not consistent with that of a previous molecular study [37]. This is probably due to the use of different molecular markers and phylogenetic analysis methods. Therefore, more molecular markers and phylogenetic analysis methods should be applied in future studies to resolve this inconsistency. The medianjoining networks showed that private haplotypes of H.
bovis and H. sinense were derived from dominant haplotypes (Fig. 2), and the percentages of private haplotypes were high (H. bovis: 76.5 % of total haplotypes; H. sinense: 87.0 % of total haplotypes), which indicates that populations of the two Hypoderma species were closely related respectively, and speciation within Hypoderma might be relatively complex.
AMOVA indicated that genetic structure was substantially higher within than among populations (within populations: 88.12 % for H. bovis and 92.17 % for H. sinense; among populations: 11.88 % for H. bovis and 7.83 % for H. sinense). Therefore, the majority of H. bovis and H. sinense genetic differentiation was intra-population. This result may be caused by gene flow. QTP uplift, which resulted in topography changes [38], limited gene flow among the populations, and might have led to gene flow that primarily occurred within the populations; this was likely exacerbated because Hypoderma species cannot fly long distances, and adults live for a very short time (only 5-6 days) [39]. Genetic structure among populations was significant (P < 0.001) for H. bovis and H. sinense based on AMOVA (Table 2), which indicates that significant population differentiation occurred between populations of these two species. F ST is used to assess genetic differentiation among closely related populations [40]. Our assessment of population genetic structure using the F ST index revealed that the range of pairwise F ST values is narrower in H. bovis (0.009 to 0.190) than in H. sinense (−0.037 to 0.210), indicating greater genetic differentiation among H. sinense than in H. bovis populations. In addition, permutation test showed that H. sinense had greater geographic structure of h than H. bovis. Geographic structure of natural populations is determined by many factors, such as life history, population size, ecological traits, habitat, and historical events [41,42]. Our results may be caused, in part, by inconsistent habitat in this study for H. bovis and H. sinense. However, geographic genetic structure was previously shown to negatively correlate with dispersal abilities [43]. In the present research, H. sinense exhibited greater geographic genetic structure than H. bovis, indicating that dispersal ability is higher in H. bovis than H. sinense. This finding is consistent with the life   [1].
Overall, these results indicate that genetic structure may differ between Hypoderma species with differences in habitat preference and dispersal ability. In this study, the neutrality tests were significantly negative for H. bovis and H. sinense populations, which indicates that population expansion events may have occurred in the demographic history of H. bovis and H. sinense. Additionally, the mismatch population test statistics rg and SSD for both H. bovis and H. sinense populations were small and not statistically significant; this indicates that the sudden expansion model, which corresponded with a unimodal curve of the mismatch distribution analysis, could not be rejected (Fig. 3).   Table 1 Overall, these analyses indicated that a demographic expansion had occurred in the H. bovis and H. sinense populations. τ values reflected estimated expansion times of 0.049 Mya for H. bovis and 0.056 Mya for H. sinense, which correspond to the late Pleistocene [44,45]. These expansions probably occurred because H. bovis and H. sinense populations on the Tibetan Plateau experienced geological changes and climatic oscillations during the QTP uplift and Quaternary glaciation, which may have also led to population differentiation of H. bovis and H. sinense. After QTP uplift and Quaternary glaciation, many private haplotypes may have been derived from the dominant H. bovis and H. sinense haplotypes during expansion phases, which resulted in the starlike haplotype networks (Fig. 2).
The COI gene is a global molecular clock gene [46]. In this study, BEAST analyses of COI results estimated that the three Hypoderma species (H. bovis, H. lineatum, and H. sinense) diverged approximately 4.5 Mya, which indicates a late Miocene-Pliocene split among H. bovis, H. lineatum and H. sinense [47,48]. Climatic changes during the late Miocene-Pliocene might have played an important role in the divergence of Hypoderma spp. Furthermore, the divergence time between H. sinense and H. lineatum was estimated to be 3.7 Mya, which mainly corresponded to the rapid uplift of the QTP approximately 3.6 Mya [49] and therefore indicates that rapid uplift of the QTP could have greatly influenced the divergence of Hypoderma spp. More importantly, the rapid uplift of the QTP changed topography, strengthened the East Asia monsoon, and modified global climate [50][51][52], which may have led to the divergence and speciation of many organisms [53][54][55]. Therefore, our results indicate that the rapid uplift of the QTP led to habitat isolation, mutation accumulation, and fragmentation in Hypoderma populations, and eventually caused speciation of H. bovis, H. lineatum and H. sinense.
Hypoderma spp. play a critical role in production losses and susceptibility of cattle and yaks to disease [7,56,57]. Our study may be a first step toward a better understanding of Hypoderma evolutionary history and speciation, and provides important information for the future study of epidemiological surveillance and hypodermosis control on the QTP.

Conclusions
This is the first characterization of the genetic diversity and population structure using mitochondrial COI sequences of H. bovis and H. sinense populations in Qinghai Province, China. The results support the distinction of the two species of Hypoderma based on genetic diversity and divergence. Most genetic differentiation of H. bovis and H. sinense was found within populations, which may have been caused by QTP uplift and life history of the species. Further research including more molecular markers, increased sampling, and different phylogenetic analysis methods are necessary to elucidate genetic differentiation of Hypoderma spp. in more detail. In addition, the current findings provide fundamental evolutionary information regarding H. bovis and H. sinense. These findings also provide a molecular baseline for the control and elimination of Hypoderma spp. on the QTP.