INTRODUCTION
Hard ticks are obligate ectoparasites, and seem to be second in importance only to mosquitoes as vectors of human and animal diseases [1]. Tick-borne diseases cause a huge loss to the livestock industry and increase the risk of disease such as Lyme disease, babesiosis, human granulocytic ehrlichiosis, forest encephalitis, spotted fever, anaplasmosis, and Crimean–Congo hemorrhagic fever [2–4]. All species are exclusively hematophagous in all feeding stages. Hard ticks are distributed worldwide with their hosts range from wild to domestic vertebrates except fishes.
Traditionally, classifications and phylogenetic inferences for Ixodidae were based on morphological, biological and ecological characteristics, often suggesting host specificity as the main factor [5,6]. However, methods for species determination are limited when taxa are morphologically very similar, specimens are damaged, and in frequent cases where immature stages are not described or are engorged [7].
Molecular systematics offered new possibilities to improve species recognition in hard ticks. ITS, 18S rDNA, 28S rDNA and other mitochondrial rDNA genes have been used to study these organisms and have played an important role in analyzing the classification and phylogenetics of hard ticks [8–10]. However, compared to the number of species of hard ticks, the extent of these studies are very limited [11].
Until recently, there has been little effort to standardize the methods for molecular identification of hard ticks, and no one gene has been formally selected as an admitted DNA marker to deal with problems of classification and phylogenetics in hard ticks. So, we chose the mitochondrial cytochrome c oxidase subunit 1 (CO1) gene fragment as a candidate molecular marker, and collected 194 samples (from 67 species of 7 genera) of hard ticks. Intra- and interspecies genetic divergences were assessed using the Kimura 2-parameter (K2P) distance model. Phylogenetic tree were performed to analyse their relationship at evolutionary level.
MATERIALS AND METHODS
Sample collection
Ticks used in this study were collected from field sites and different hosts in various regions of China (Table 1). After morphological identification, ticks were stored in 100% ethanol and conserved at 4°C. Only male and unfed adult specimens were used.
DNA extraction, PCR amplification, and sequencing of CO1
DNA was extracted from the ticks using a tissue kit (Qiagen AG, Basel, Switzerland) according to the manufacturer’s instructions. Each sample was cut with sterile scissors within a 1.5 ml microtube. After digestion with proteinase K (20 mg/ml), the samples were applied to the columns for DNA absorption and washing. Finally, the DNA was eluted in 100 ml of eluting buffer provided in the kit and stored at −20°C. The primers used for PCR were LCO1490 (5′-GGTCAACAAATCATAAAGATA-TTGG-3′) and HCO2198 (5′-TAAACTTCAGGGTGACCAAAAAATCA-3′) [12]. The 5′ region of CO1 was amplified using the following thermal cycling program: 94°C for 5 min, 35 cycles at 94°C for 1 min, 53°C for 1 min, and 72°C for 1 min, followed by a final extension at 72°C for 8 min. PCR products were purified using a PCR purification kit (Takara, Shiga, Japan). Sequencing reactions were resolved on automated DNA sequencer.
Data from GenBank
Some CO1 sequences from the hard ticks were downloaded from GenBank. Sequences <500 bp in length, with ambiguous bases (more than 15 ‘Ns’), or those belonging to unnamed species (sequences with ‘sp.’ in the species name) were removed from the analysis. In addition, we checked all the sequences using BLAST analysis (E-value<0.001) to make sure that there were no host sequences in our data. The selected sequences were used to construct analysis datasets.
Sequence analysis
The CO1 sequences were translated into amino acids with the program MEGA 4.0 in order to exclude sequencing errors and to avoid the inclusion of pseudogene sequences in the datasets. The sequences were aligned using ClustalW [13], and genetic distances were computed using MEGA 4.0 according to the K2P distance model. The maximal/mean/minimum intra- and interspecies distances were used to represent species level divergence. Meanwhile, the maximal/mean/minimum intra-and intergenus distances were calculated to evaluate the genus level variation. Then a neighbor joining (NJ) tree was constructed and the genetic K2P distance was calculated within species and genera using MEGA 4.0. Evaluation of statistical confidence was based on 1,000 non-parametric bootstrap replicates. One soft tick isolate was used as the outgroup of phylogenetic tree.
RESULTS
Data acquisition
We collected 194 samples (36 from this study, 158 from GenBank) from 67 species and 7 genera of hard ticks (Table 1 and Supplementary Table S1). The mitochondrial CO1 region of all samples collected in China was successfully amplified using PCR. The length of the PCR product was 707 bp. As some sequences of the CO1 gene obtained from GenBank were shorter than 707 bp, all sequences were aligned with a consensus length of 586 bp, and no insertions, deletions, or stop codons were observed in any sequence. The sequences acquired in this study have been deposited in the GenBank database with accession numbers JQ737066-JQ737128.
Genetic divergence and gap
Using the K2P model, sample divergences at various taxonomic levels are shown in Tables 2 and 3. As expected, the genetic divergence increased with higher taxonomic ranking: 0.001±0.001 to 0.016±0.003 at intraspecies level, 0.002±0.001 to 0.248±0.023 at interspecies level, 0.005±0.002 to 0.175±0.011 at intragenus level, and 0.186±0.012 to 0.243±0.016 at intergenus level. The Bothriocroton showed the lowest mean intraspecies divergence (0.005±0.002), while Rhipicephalus showed the highest mean intraspecies divergence (0.062±0.039) (Fig. 1). The largest ratio between the average intra- and interspecies divergence was in the Ixodes with a 7.5-fold difference, and the lowest ratio was in the Dermacentor with a 2.4-fold difference. As shown in Fig. 1, there was not a distinct gap between the distribution of the intra- and interspecies divergence. The overlapping regions were mainly distributed in the R. turanicus, Hya. dromedarii, D. marginatus, D. silvarum, and A. testudinarium.
Phylogenetic tree
The NJ tree of the overall analysis is shown in Fig. 2. The phylogenetic relationship at the genus level was well resolved with the exception of Amblyomma. From the tree, Hyalomma, Dermacentor, Amblyomma, and Rhipicephalus formed 1 clade. Bothriocroton and Haemaphysalis formed another clade. Ixodes as distinct difference at morphous to other hard ticks, formed a third clade. However, at a species level, 9 species (13.4%) did not form a monophyletic group. They were Hya. dromedarii, Hya. marginatum, Hya. asiaticum asiaticum, D. marginatus, D. silvarum, A. testudinarium, R. microplus, and H. flava.
DISCUSSION
In this study, the mean sequence divergence in hard ticks (0.197±0.063) is higher than that observed in other organisms [14–16]. Such high values of genetic distance reflect possible biological diversity within the Ixodidae. Such as the distance between Amblyomma testudinarium (HM193893) and A. testudniarium (HM193895) was 0.112±0.010, and they were in different clades of the phylogenetic tree. However, Rhipicephalus microplus and Dermacentor marginatus also gave similar results. The reason may be geographic variation or comprise cryptic species [17]. Additionally, the distance between the species Dermacentor everestianus (JQ737079) and D. niveus (JQ737080) was only 0.004±0.002, and also formed into 1 clade. Therefore, these analyses might indicate hybridization or a misidentification among these species.
The CO1 gene appears to be an informative molecular marker on several taxonomic scales, but particularly at the species level [18]. Our analysis shows a general increase in the molecular divergence of CO1 with taxonomic rank. The diversity within species is especially high, with a maximum of 0.116±0.015. It makes CO1 suitable for investigating intraspecies variation. DNA barcoding assumes that the genetic distances between species are greater than within species. In that way, clusters of similar sequences represent species, clearly separated from other clusters (species) [19]. However, there also 30 samples where the maximum interspecies distance was larger than the minimum interspecfic distance. This means that the gap might be absent in these samples because of insufficient variation between them [20,21]. From the NJ phylogenetic tree, nine of the 67 species (13.4%) examined in this study (Hya. dromedarii, Hya. marginatum, Hya. asiaticum asiaticum, Hya. truncatum, D. marginatus, D. silvarum, A. testudinarium, R. microplus, and H. flava.) did not form a monophyletic group. Hya. asiaticum asiaticum and Hya. dromedarii shared similar morphologic characters from capitulum, scutum, Haller’s organ, peritrematal plate, the first caruncle, coax and spur of feet of adults and larval stages. Ecologically, these 2 species also share the same desert intertidal area. They are 2 different species proved by previous studies [22–24]. However, they formed one clade in this study. This phenomenon was also found for other hard ticks. For example, Hya. dromedarii, Hya. marginatum and Hya. truncatum formed a complex clade. These results agreed with some studies using mt 12S rDNA, 16S rDNA or ITS, in which Hyalomma spp. shown high divergence distance and low bootstrap value [25,26]. As many results indicated that there is a high diversity in hard ticks [27,28].
This study provides that using the CO1 gene is a potential tool for species identification in Ixodidae. However, it is inadequate to use a single mitochondrial gene (CO1) for DNA taxonomy. Therefore, an integrative approach is needed to combine nuclear and mitochondrial genes, morphological characters, and ecological information into further studies of hard ticks.