Coevolutionary analyses of the relationships between piroplasmids and their hard tick hosts

Host–parasite coevolution is a key driver of biological diversity. To examine the evolutionary relationships between piroplasmids and their hard tick hosts, we calculated the molecular clock and conducted phylogenetic analyses of both groups. Based on our results, we conclude that the divergence time of piroplasmids (∼56 Mya) is later than divergence time of their hard tick hosts (∼86 Mya). From analyses of the evolution of both piroplasmid and vector lineages and their association, we know that hard ticks transmit piroplasmids with high genus specificity and low species specificity.


Introduction
Coevolution between parasites and their hosts is a widely recognized phenomenon in many host-parasite systems. It is also an essential component of evolution and the diversity of life. Piroplasmids (Babesia spp. and Theileria spp.) are protozoan hemoparasites of great economic, veterinary, and medical importance worldwide. They infect a diverse array of mammalian and a few avian hosts (Lack et al. 2012;Schnittger et al. 2012). Their definitive hosts, species of hard ticks, are obligate ectoparasites, and seem to be second in importance only to mosquitoes as vectors of human and animal diseases ( Fig. 1) (Fivaz et al. 1992;Jongejan and Uilenberg 2004). Hard ticks are exclusively hematophagous and are distributed worldwide. Each piroplasmid species is transmitted by different hard tick species, even those piroplasmids with close phylogenetic relationships (Yin et al. 2002).
The diversity of the relationships between host and parasite population structures (such as host switching, cospeciation, duplication, and extinction) usually makes it difficult to draw general conclusions and underlines the need for more studies on host-parasite systems (Garamszegi 2009;Ricklefs et al. 2004). Therefore, systematic studies of the diversification of piroplasmids in relation to their vector can provide insight into the coevolutionary events that shaped present-day interactions among these taxa. (Agudo et al. 2010).
Previous studies estimated that the date of origin of the Apicomplexa was between 1100 and 500 Mya based on molecular clock estimations (Berney and Pawlowski 2006;Escalante and Ayala 1995), predating by far the presumed origin of the ticks as well as the radiation of mammals and birds (120-90 Mya) (Escalante and Ayala 1995). Given that apicomplexan phylogeny correlates with the taxonomy of their definitive host, hard ticks probably represent the primordial host of piroplasmids (Barta 1989). Thus, the question arises as to when an ancient piroplasmid ancestor might have selected the hard tick as its vector. MtDNA offers several advantages in comparison to nuclear DNA: rapid evolution, limited exposure to recombination, lack of introns, and high copy number. These characteristics of mtDNA are important for routine amplification by polymerase chain reaction (PCR) and use as a molecular marker for phylogeny (Gou et al. 2012). In this study, we used the cytochrome oxidase I (COI) gene to calculate the molecular clock of piroplasmids and hard ticks, and also reconstructed molecular phylogenies for the parasites and their hosts to assess the extent of coevolution in this association.

Taxa sampling
In total, 11 piroplasmid species were collected from five provinces and three autonomous regions (Gansu, Henan, Yunnan, Hebei, Sichuan, Xinjiang, Ningxia, and Inner Mongolia) of China. Cattle and sheep were inoculated with 10-mL cryopreserved infectious blood containing each isolate. When the level of parasitemia reached more than 5%, blood from the jugular vein was collected into tubes using heparin as an anticoagulant. Parasites were isolated as described in previous studies, and stored at À70°C for DNA extraction.
Correspondingly, 10 species of hard tick were also collected from field sites and different hosts in various regions (Gansu, Xinjiang, Xizang, Guangxi, Qinghai, Hubei, and Henan) of China. After morphological identification, these ticks were stored in 100% ethanol and conserved at 4°C. To eliminate interference from the specimens, only female and unfed adult hard ticks were used in the study.
DNA of hard ticks was extracted using a tissue kit (Qiagen). 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 elution buffer provided in the kit and stored at À20°C. The universal primers used for PCR were as follows: LCO1490: 5′-GGT-CAACAAACATA AAGATATTGG-3′ and HCO2198: 5′-TAAACTTCAGGG TGACCAAAAAATCA-3′ (Vrijenhoek 1994).
PCR reactions were run using the following thermal cycling program: 94°C for 5 min, 35 cycles at 94°C for 1 min, 45-54°C for 1 min, and 72°C for 1 min, followed by a final extension at 72°C for 8 min. The 50-lL PCR reaction mixture included 37.5 lL of ultrapure water, 5 lL of 59 PCR buffer, 1 lL of each primer (20 lmol/ L), 4 lL of dNTPs (10 mmol/L), 0.5 lL of Taq polymerase (5 U), and 1 lL of the DNA template. PCR products were purified using a PCR purification kit (Takara Bio Inc, Japan). Sequencing reactions were resolved on an ABI 3730 automated DNA sequencer.

Sequence alignment
Edited sequence contigs for both COI genes (from piroplasmids and from hard ticks) were aligned using the software program MEGA 4.0 (Tamura et al. 2007) with default settings. To avoid bias in refining alignments, ambiguous alignment positions, including the beginning and end regions, were excluded. Furthermore, nucleotide sequences were translated into amino acid sequences using the invertebrate mitochondrial code in MEGA 4.0 to prove the correctness of the DNA sequences.

Saturation test
Given that saturation in substitutions can lead to incorrect phylogenetic inferences, the COI genes from piroplasmids and hard ticks were evaluated for transition and transversion substitutions by DAMBE 5.3.3 (Xia and Xie 2001).

Calibration point
Owing to a lack of fossil piroplasmids and potentially many unsampled lineages, we generated a combined alignment of piroplasmid and Plasmodium COI sequences and calibrated the phylogeny with an age of 12.82 million years (95% highest posterior density [HPD] 9.93-19.49 Mya) for the three most recent common ancestors (MRCA) of Plasmodium spp., as reported previously (Ricklefs and Outlaw 2010).

Divergence time and phylogenetic analyses
ModelTest 3.7 (Posada 2005) was used to determine the appropriate DNA substitution model and gamma rate heterogeneity with the akaike information criterion (AIC). For both matrices, ModelTest was performed for each gene region. As a result, the GTR + I + G was determined as the best-fitting statistical model. BEAST 1.7.1 (Drummond and Rambaut 2007) was used to conduct relaxed-clock Bayesian analyses to estimate topology and divergence times simultaneously. Evolutionary rates along branches followed an uncorrelated lognormal distribution, and a Yule speciation process was imposed for all analyses. Two independent MCMC (Markov chain Monte Carlo) runs were performed, of 10 million and 5 million generations, with sampling every 1000 generations. These two separate runs were then combined (following the removal of 10% burn-in) using LogCombiner 1.7.1. Adequate sampling and convergence of the chain to stationary distribution were confirmed by inspection of MCMC samples using Tracer 1.5 (Drummond and Rambaut 2007). The effective sample size (ESS) values of all parameters were greater than 200, which were considered a sufficient level of sampling. The sampled posterior trees were summarized using TreeAnnotator 1.7.1 to generate a maximum clade credibility tree (maximum posterior probabilities [PP]) and calculate the mean ages, 95% HPD intervals, PP, and substitution rates for each node. The BEAST topology was visualized with FigTree 1.3.1 (Drummond and Rambaut 2007).

Data characteristics
We collected 14 species of piroplasmid (eight species of Babesia and six species of Theileria). Nucleotide sequence editing and alignment resulted in COI data partitioning of 950 characters. In the aligned sequences, 466 sites were variable and 388 sites were parsimony informative. First-, second-, and third-codon positions showed AT biases of 58.6, 60.6, and 84.2%, respectively. Codon usage frequency was 33.3%.
In total, 31 species of hard tick were used in this study (six species of Ixodes, five species of Haemaphysalis, three species of Dermacentor, seven species of Hyalomma, three species of Amblyomma, and seven species of Rhipicephalus). All sequences were aligned with a consensus length of 586 bp, 303 sites were variable and 260 sites were parsimony informative. First-, second-, and third-codon positions showed AT biases of 57.5, 57.5, and 85.9%, respectively. Codon usage frequency was 32.3%.

Saturation effects and substitution patterns
Comparison of the transition and transversion frequencies of piroplasmids and hard ticks revealed that the frequency of transition was lower than transversion (Fig. 2). For the test of substitution saturation, which examined whether the observed Iss (index for saturation substitution) was significantly less than Iss.c (critical Iss), we found that the observed Iss (0.58) was less than the Iss.c (0.76) for piroplasmids, and Iss (0.27) was less than the Iss.c (0.74) for hard tick. These results indicated that the two COI genes were appropriate for phylogenetic reconstruction.

Estimating divergence time
The analyses estimated the divergence time of piroplasmids to be during the Paleocene,~56.48 Mya (95% HPD 28.17-86.87); node A represented Theileria with a node age of 23.38 Mya (95% HPD 11.11-36.71) and node B represented Babesia with a node age of 25.74 Mya (95% HPD 12.75-40.73) (Fig. 3, Table 1). The Babesi clade formed slightly earlier than that of Theileria, and the two events both began during the late Eocene.
The analyses estimated that hard ticks began to diverge during the Cretaceous,~86 Mya (95% HPD 61.11-104.16 Mya). The nodes C, D, E, F, G, and H in Figure 4 represent the times at which the following genera formed: Ixodes (

Cophylogenetic analyses
When the two phylogenetic trees of the ticks and parasites are compared (Fig. 5), we can see that the genus divergence time of the parasite (Theileria and Babesia) is later than their invertebrate hosts (Ixodes, Haemaphysalis, Dermacentor, Hyalomma, Amblyomma, and Rhipicephalus). Moreover, this trend is also lasting at the species divergence time for majority of parasites and their definitive hosts. Furthermore, the Bayesian topologies of the combined parasite sequence data set and host data set were used in the cophylogenetic analysis and considered as species trees. This revealed that parasite and host phylogenies were not congruent (Fig. 6).

Discussion
In the study of evolution, estimates of divergence time can provide crucial information to speciation events (Hayakawa et al. 2008;Ricklefs and Outlaw 2010). In previous studies, the divergence time of the Apicomplexa was estimated to be between 1100 and 500 Mya based on molecular clock estimates (Berney and Pawlowski 2006;Escalante and Ayala 1995). Published divergence time of piroplasmid varies from 820 Mya to 17 Mya and is surrounded by a degree of uncertainty. We calculated the divergence time to be 56.48 Mya (95% HPD 28.17-86.87), which is earlier than the result from a study based on 18S rRNA (17.11 Mya) (Lack et al. 2012). This discrepancy might be because of the difference in evolutionary rate of the two genes. The gene encoding the 18S rRNA is more conservative than is COI, and has a slower evolutionary rate (Gou et al. 2012;Nassonova et al. 2010). Published divergence time of ticks ranges across two evolutionary epochs, including the Triassic (248-206 Mya) (Balashov 1989(Balashov , 1994 and Cretaceous (146-65 Mya) (Black and Piesman 1994). However, there are few studies about the divergence time of hard ticks. We estimated the divergence time of hard ticks to be 86 Mya (95% HPD 61.11-104.16) based on three dominica amber samples, which is later than previous estimates (90-120 Mya) (Klompen et al. 2005;Mans et al. 2002). We can see that the both genus and species divergence time of the parasite are later than their hosts hard ticks after compared two phylogenetic trees.
The other objective of this study was to assess the degree of congruence between the parasite phylogeny and   2010; Mazyad et al. 2010;Schwint et al. 2008), these hard tick species are always from the same genus. This means that hard ticks, as a vector, transmitted piroplasmids with high genus specificity and low species specificity. Although environmental effects are often excluded from experimental coevolutionary studies and are included in theoretical models as 'noise', some studies under laboratory conditions have shown that the environment in which hosts and parasites interact can substantially affect the strength and specificity of selection (Wolinska and King 2009). However, geographical isolation as an element of the environment should also be discussed in this study of coevolution between piroplasmids and hard ticks. For instance, both B. major and its vector H. punctata are only found in the Xinjiang Uygur Autonomous Region of China. This region is located in the Junggar basin and is isolated from other areas by desert and high mountains. A similar case is also observed with T. sinensis and its vector H. qinghaiensis. These two organisms are mainly distributed in Gansu and Qinghai provinces in northwest of China (Yin et al. 1996), areas characterized by high altitude, drought, and loess. These two examples of geographical isolation have resulted in piroplasmids and hard ticks forming a single biological system over a long period of evolutionary time.

Supporting Information
Additional Supporting Information may be found in the online version of this article: