- Open Access
Phylogeography and population differentiation in Hepatozoon canis (Apicomplexa: Hepatozoidae) reveal expansion and gene flow in world populations
Parasites & Vectors volume 14, Article number: 467 (2021)
Hepatozoon canis is a protozoan transmitted to dogs and other wild carnivores by the ingestion of ticks containing mature oocysts and is considered the principal cause of canine hepatozoonosis in the world. Here, we examined ribosomal RNA 18S gene sequence variation to determine the genetic differences and phylogeographic diversity of H. canis from various geographical areas around the world.
We used 550 publicly available sequences of H. canis from 46 countries to assess haplotype relationships, geographical structure, genetic diversity indices, and relationships among populations. We performed neutrality tests and pairwise comparisons of fixation index (FST) values between groups and pairwise comparisons of FST values between populations. To determine whether populations are structured, analyses of molecular variance (AMOVAs) and spatial analysis of molecular variance (SAMOVA) were performed.
The dataset of H. canis yielded 76 haplotypes. Differentiation among populations indicated that there is no phylogeographical structure (GST = 0.302 ± 0.0475). Moreover, when samples were grouped by continents a significant FST was obtained, meaning that populations were genetically differentiated. The AMOVA showed that 57.4% of the genetic variation was explained by differences within populations when all locations were treated as a single group and revealed that there is no population structure when populations are grouped into two, three, and four groups (FCT, p > 0.05), suggesting that dispersal between populations is high. SAMOVA revealed significant FCT values for groups K = 5. The Tajima’s D and Fu’s Fs show that populations have undergone recent expansion, and the mismatch distribution analysis showed population expansion (multimodal distribution).
The current molecular data confirmed that H. canis does not show phylogeographic or population structure. The haplotypes exhibit low genetic differentiation, suggesting a recent expansion due to gene flow among populations. These results provide pivotal information required for future detailed population genetic analysis or to establish control strategies of this parasite.
The vector-borne pathogens of the genus Hepatozoon infect a wide range of species of mammals, reptiles, amphibians, and birds [1,2,3,4,5]. These parasites are transmitted by different groups of arthropods including ticks, lice, mosquitoes, mites, sand flies, tsetse flies, kissing bugs, and leeches [6,7,8]. In the past years, the prevalence of Hepatozoon species affecting domestic and wildlife animals has increased worldwide; therefore, the study of Hepatozoon has gained relevance in the veterinary field . Hepatozoon americanum and H. canis are protozoa transmitted to dogs and other wild carnivores by ingestion of ticks containing mature oocysts [9, 10]. Although both H. americanum and H. canis are reported in canids, the latter is widely spread and is considered the principal cause of canine hepatozoonosis in the world . The main vector for H. canis is the brown dog tick Rhipicephalus sanguineus ; however, other tick species are suspected to be possible vectors, for example Haemaphysalis flava and Haemaphysalis longicornis in Japan , Amblyomma ovale in Brazil , Amblyomma mixtum and Rhipicephalus turanicus in Mexico [1, 14], and Ixodes ricinus ticks in Italy . Following ingestion of infected ticks, H. canis sporozoites spread via the bloodstream and lymph to several organs including the spleen, bone marrow, lung, liver, and kidney, infecting leukocytes and parenchymal tissue cells.
Wild canids may be an important reservoir for H. canis, thus representing a possible health hazard to domestic dog populations . Wild canids do not develop clinical signs; in contrast, infection in dogs affects several organs resulting in anemia and lethargy [16, 17]. Currently, H. canis infects a wide range of carnivorous hosts around the world including dogs, jackals, foxes, opossums, and domestics cats [1, 15, 16, 18,19,20].
Molecular sequences obtained from collected samples can provide information about past evolutionary events. Relevant probabilistic models of molecular evolution enable reconstructing ancestral sequences to a sample of taxa, and phylogenetics provides an adequate framework for the reconstruction of ancestral sequences . Phylogeographic studies of several taxa have revealed a pattern of lineage divergence associated with contraction to and expansion populations associated with the glacial and interglacial cycles  as well as with the formation of geographic barriers . Mutation rates and speed of selection differ among taxa, and in many cases, their dispersion to new areas is mediated by anthropogenic factors such as the mobility of the hosts towards different geographic regions, as has been shown by several studies (i.e. Vibrio vulnificus , canine distemper virus , Corynosoma australe , Angiostrongylus cantonensis, and A. malaysiensis ) including those performed with some apicomplexans (Plasmodium knowlesi  and Toxoplasma gondi ).
Despite the relevance of many parasitic diseases, the studies in genetic and phylogeographic diversity are few compared with other taxa like mammals, amphibians, birds, and plants. Molecular detection of H. canis is widely reported in different studies in several regions around the world; however, the genetic diversity and phylogeography of these hemoparasites is understudied [1, 15, 16]. In this article, we examined ribosomal RNA 18S gene sequence variation to determine the genetic differences and phylogeographic diversity of H. canis from various geographical areas around the world. For H. canis, we can infer that dispersion around the world could be related to the movement of dogs associated with human migration; hence, we hypothesized there is no population or phylogeographic structure. Thus, we aimed to test whether the migratory movements in humans is associated with the population structure of eukaryotic parasites, using Hepatozoon as a model of study.
To define the final dataset of sequences included in this study, several filtering steps were applied as follows: (1) a total of 1170 sequences were downloaded from the ENA/GenBank database using as searching criteria “Hepatozoon canis” to discard any sequences not identified to species level (Hepatozoon sp.); (2) 55 sequences without complete metadata lacking host and country of origin were discarded; (3) sequences with duplicated accession number (redundant) or missed annotations from parasites different from H. canis were not selected; (4) considering an observed wide range in sequence length between 140 bp and 3.1 kbp, sequences < 500 bp were omitted from further analysis. The remaining sequences had a median of 602 bp (533–677). During this step, sequences from genetic regions other than the ribosomal 18S gene were discarded. (5) Finally, the remaining 618 sequences were aligned, and sequences with coverage < 50% (in the final alignment with length of 589 bp) were discarded. The custom R scripts used for the sequence selection are available at: https://github.com/VictorHJD/Hepatozoon_Phylogeography. The script for selection was run in R version 4.0.3 .
The final cured dataset contained 550 sequences (samples) from the rRNA 18S gene spanning the V4 hypervariable region of H. canis (Table 1). Sequences belonged to populations from 46 countries (populations) on four continents (Fig. 1, Table 1). These included 109 sequences from Africa, 61 sequences from America, 116 sequences from Asia, and 264 from Europe (Table 1, Additional file 1: Dataset S1). The collected sequences corresponded to isolates from several hosts: dogs (Canis lupus familiaris) (n = 339), wildcats (Felis silvestris) (n = 2), red foxes (Vulpes vulpes) (n = 94), pampas fox, (Lycalopex gymnocercus) (n = 1), black-backed jackals (Canis mesomelas) (n = 7), golden jackals (Canis aureus) (n = 15), opossums (Didelphis albiventris) (n = 2), and capybara (Hydrochoerus hydrochaeris) (n = 1) and several tick vector species: Rhipicephalus sanguineus (n = 46), Rhipicephalus (Boophilus) microplus (n = 1), Rhipicephalus tiranicus (n = 7), Rhipicephalus spp. (n = 3), Haemaphysalis adleri (n = 2), Haemaphysalis bispinosa (n = 6), Haemaphysalis concinna (n = 2), Haemaphysalis longicornis (n = 1), Ixodes canisuga (n = 2), Ixodes hexagonus (n = 1), Ixodes ricinus (n = 10), Dermacentor marginatus (n = 2), Dermacentor reticularis (n = 1), Dermacentor spp. (n = 1), Amblyomma mixtum (n = 1), and three unidentified ticks. All sequences were automatically aligned using Clustal X  and manually edited with PhyDE-1v0.9971 .
To infer genealogical relationships among H. canis haplotypes, statistical parsimony networks (haplotypes network) were constructed using the program TCS 1.2.1 . The aligned sequences were mapped to estimate gene genealogies with gaps treated as missing data and a 95% connection probability limit. Loops were resolved following the criteria given by Pfenninger and Posada . The full dataset was processed blindly and treated independently from their host of origin.
The most likely number of genetically differentiated clusters was estimated using BAPS (Bayesian Analysis of Population Structure, version 6.0) . BAPS assesses the most likely number of genetically different clusters using the module for linked molecular data. The codon linkage model, appropriate for our sequence data, was applied. For the chosen marker, the probability of a different number of genetic clusters (K = 2 to K = 10) under the independent loci model was evaluated. Two independent runs with ten replicates for each K, accepting the partition with the K value that had the highest likelihood and posterior probability (PP) provided the final result.
Geographical structure, genetic diversity indices, and relationships among populations
Population diversity indices for unordered (hS, hT) and ordered haplotypes (vS, vT) and differentiation parameters (GST, NST) were estimated using software PERMUT 1.0 . The within-population diversity (hS), total diversity (hT), geographical average haplotype diversity (vS), geographical total haplotype diversity (vT), level of population differentiation at the species level (GST), and an estimate of population subdivisions for phylogenetically ordered alleles (NST) were calculated. GST and NST are often used to assess the geographical structure affecting population differentiation. Significant differences between NST and GST parameters were tested with 10,000 permutations. If NST is significantly larger than GST, this means that two alleles sampled from the same region are phylogenetically more closely related than two alleles sampled from different regions, evidencing the presence of a significant phylogeographical signal in the data .
Neutrality tests and pairwise comparisons of FST values between continents were calculated using ARLEQUIN 3.11 . Pairwise comparisons of FST values between populations were calculated with 1000 permutations. Populations with three samples or fewer were excluded as required for this analysis.
To determine whether populations were structured, analyses of molecular variance (AMOVAs)  were performed based on pairwise differences. The AMOVA was carried out considering all the populations across the sampling area as a single unit. Furthermore, hierarchical AMOVA was performed with populations treated as (i) grouped into two, New world and Old world (America, Europe + Africa + Asia), or (ii) three (America, Africa, Asia + Europe), according to the hypothesis that H. canis arises in Eurasia and spreads to Africa and America, or (iii) four groups (America, Africa, Asia, Europe) similar to the episode where the four geographic regions (continents) have been isolated long enough to differentiate populations. The AMOVAs were performed using ARLEQUIN 3.11 with 16,000 permutations to determine the significance of each AMOVA. Additionally, spatial analysis of molecular variance (SAMOVA) was performed using SAMOVA 1.0  to identify clusters of populations that are geographically homogeneous and genetically differentiated.
Molecular diversity indices, including the number of haplotypes (NH), gene diversity (h), nucleotide diversity (\(\pi\)), and pairwise comparisons of FST values between populations, were calculated using ARLEQUIN 3.11 . We used mismatch distributions to assess population expansion and time variation in effective population size (Ne). The historical population expansions were analyzed with Tajima’s D (DT)  and Fu’s (FS) . Negative values in DT and FS that result from an excess of rare alleles present at low frequencies indicate that populations have undergone recent expansion.
Neutrality tests were performed in ARLEQUIN 3.11, with 10,000 permutations . Mismatch distributions  were calculated using the sudden expansion model of Schneider & Excoffier , with 1000 bootstrap replicates. Mismatch distributions  and Harpending’s raggedness index (Hri; ) were calculated in DnaSP v.6 . The validity of the sudden expansion assumption was determined using the sum of squares differences (SSDs) and Hri , both of which are higher in stable, non-expanding populations . Low and non-significant values of Hri and SDD indicated a good fit between the observed and the expected values of the sudden expansion model.
Ramos-Onsins and Rozas R2 statistics were calculated based on the difference between the number of singleton mutations and the average number of nucleotide differences, where small positive values of R2 are expected under a scenario of population expansion . Significance was evaluated by comparing observed values with null distributions generated by 10,000 replicates, using the empirical population sample size and the observed number of segregating sites in the “pegas” package of R version 4.0.3 .
Finally, the Bayesian skyline plot  method implemented in BEAST 2.6.4  was performed to further explore the demographic history and to estimate the timing of the population expansion in H. canis populations (four continents and all continents). The 550 sequences were used to run the analysis using the coalescent Bayesian skyline tree prior. This coalescent-based approach estimates the posterior distribution of effective population size (Ne) at intervals along a phylogeny, thereby allowing one to infer population fluctuations over time. The haplotype sequence alignments were imported in the program BEAUti2 (included in BEAST 2.6.4) in nexus format. The best nucleotide substitution model was estimated using Jmodeltest 2.0 . The HKY + G substitution model was the best model identified for African, American and European H. canis, while that for Asia was TPM3uf. Finally, when populations were treated as a single group, the best substitution model was TPM3uf + G. Since the literature lacks estimations of substitution rates for H. canis, a substitution rate of 1.5% per 100 million years for the 18S rRNA used to study the lineage evolution of Toxoplasma gondii  was implemented in the analysis. Three independent runs of 10 million generations each were run, with trees and parameters sampled every 1000 iterations, with a burn-in of 10%. Results of each run were visualized using TRACER to ensure that stationarity and convergence had been reached .
A highly abundant haplotype of H. canis is widely spread worldwide
Statistical parsimony retrieved a well-resolved haplotype network that identified 76 haplotypes in the world representing four haplogroups (Fig. 1, Additional file 2: Figure S1). The aligned 18S rRNA gene dataset for 550 samples of H. canis (mean length: 586 bp) yielded 76 haplotypes with different abundance along the geography (Table 1, Fig. 1a). The haplotype H1 (20.3%, 112/550 sequences) was the most common and widespread of the individuals and in 43.4% (20/46 total populations, considering each set of sequences from the same country as a different population), retrieved in populations located in the four continents. It was followed by haplotypes H7 (18%, 99/550 sequences) in populations located in the Old World and H9 (13.3%, 73/550 sequences) retrieved in populations located in the four continents. Haplotypes H13, H24, and H16 were exclusively found in Europe, and haplotype H41 was only retrieved in Africa and Asia (Fig. 1). The haplotypes did not cluster based on their host of origin, suggesting the lack of variability of the sequences between invertebrate and vertebrate hosts (Additional file 2: Table S1, Figure S2).
Out of total sequences in H1, 13.4% (15/112) correspond to tick while 84.8% (95/112) belong to dogs and two sequences (1.7%) to red foxes. In H7, 9.09% (9/99) belong to ticks, and the 90.91% are distributed in sequences from other vertebrate hosts (74/99; 74.74% belonged to dogs, 10/99; 10.1% to red foxes, and 6/99; 6.06% to golden jackals). For H9, 12.32% (9/73) belong to ticks, 53.42% (39/73) to dogs, 30.13% (22/73) to red foxes, 2.73% (2/73) to opossum, and 1.38% (1/73) to pampas fox. Moreover, the haplotype network shows various singletons for each host (Additional file 2: Figure S2).
BAPS analyses with sequences indicated the existence of four genetic clusters (K = 4) (log marginal likelihood = − 4667.6928, PP = 1.0) (Fig. 1c). The four genetic clusters from the BAPS analysis distributed in different proportions among continents. The yellow cluster includes individuals from the Asian populations, exclusively. The light blue cluster gathers individuals from the African, Asia and European populations exclusively. The pink and magenta clusters include individuals from populations distributed in the four continents, confirming the lack of genetic structure (Fig. 1c). BAPS analysis was in agreement with the shallow genetic divergence and haplotype distribution shown in Fig. 1b.
Lack of phylogeographic structure of H. canis suggests gene flow between populations
Differentiation among populations based on 18S rRNA gene variation (mean ± SE, GST = 0.302 ± 0.0475) indicated that H. canis is genetically subdivided. Genetic diversity across all populations (hT = 0.913 ± 0.0136; vT = 0.914 ± 0.0616) was higher than the average within-population value (hS = 0.637 ± 0.0463; vS = 0.550 ± 0.0682), indicating that populations are similar among them. PERMUT analysis showed that NST was not significantly different (PERMUT: NST = 0.399 ± 0.0578 vs. GST = 0.302 ± 0.0475, perm: 10,000, p > 0.05) than GST, indicating no phylogeographical structuring.
Hepatozoon canis is genetically differentiated between populations, as pairwise comparisons of FST values were significant for the data set when samples were grouped by continents (Table 2). Europe was the continent that most differed from the other continents. The lowest differentiation was observed between Asia and Africa, while the highest difference was observed between Asia and Europe (Table 2).
The AMOVA results showed that 57.4% of the genetic variation was explained by differences within populations when all locations were treated as a single group (Table 3). The hierarchical AMOVA revealed that there is no population structure, with the lowest FCT value obtained when populations are grouped into two, three, and four groups with differences that were not statistically significant (Table 3), suggesting that dispersal of H. canis between populations is high, thus increasing gene flow. SAMOVA results revealed significant FCT values for groups between K = 5. The highest FCT value was for K = 5 (Table 3). When K = 6, FCT was smaller than the FCT of K = 5, and an additional increase in the number of K led to a dissolution of group structure, and single-population groups were formed.
Hepatozoon canis populations have a recent demographic history
Considering the hypothesis of enough divergence time to form independent populations by continent, the Tajima’s D and Fu’s Fs show a small, statistically significant, negative value for America (Tajima's D test: D = − 1.5908, SD = 0.9046, p = 0.0288; Fu’s Fs test: Fs = − 4.0733, SD = 0.18721, p = 0.0288), Asia (Tajima’s D test: D = − 2.7120, SD = 1.0353, p = 0.0000; Fu’s Fs test: Fs = − 5.2295, SD = 0.18721, p = 0.08600), and Europe (Tajima’s D test: D = − 1.6212, SD = 1.0353, p = 0.0240; Fu’s Fs test: Fs = − 13.3581, SD = 5.9440, p = 0.0020), and a small, statistically non-significant, negative value for Africa (Tajima’s D test: D = − 0.1857, SD = 1.0353, p = 0.5100; Fu’s Fs test: Fs = − 1.1153, SD = 5.9440, p = 0.4230). When populations are treated as a single group, the Tajima’s D and Fu’s Fs had a small and statistically significant negative value (Tajima’s D test: D = − 2.5002, SD = 0.0000, p = 0.0000; Fu’s Fs test: Fs = − 24.7563, SD = 0.0000, p = 0.0000), indicating that populations have undergone recent expansion, often preceded by a bottleneck (Table 4). In the mismatch distribution, SSD and Hri results were low and non-significant (Mismatch analysis test: SSD = 0.0146, SD = 0.0000, p = 0.0500; Hri = 0.0348, SD = 0.0000, p = 0.0900), considering all populations as one group (global) and grouping America (mismatch analysis test: SSD = 0.0853, SD = 0.0427, p = 0.1900; Hri = 0.2013, SD = 0.0963, p = 0.1500), Asia (mismatch analysis test: SSD = 0.0145, SD = 0.0427, p = 0.1100; Hri = 0.0589, SD = 0.0963, p = 0.0600), and Europe populations Asia (mismatch analysis test: SSD = 0.0154, SD = 0.0427, p = 0.4100; Hri = 0.3628, SD = 0.0963, p = 0.2900), providing evidence for sudden demographic expansion in the past for all groups of populations (Table 4). The mismatch distribution analysis showed a multimodal distribution when observed frequencies were compared against expected frequencies, implying population expansion (Fig. 2). Finally, R2 had small, statistically significant, positive values for Africa, America, Asia, Europa, and the global populations, indicating recent demographic expansion (Table 4, Fig. 2).
The Bayesian skyline plots suggest that population size effectively increased around 75,000–5000 years ago in all populations (Fig. 3). The African population was relatively stable over time, showing a slight increase around 5000 years ago. In contrast, the American population and Asia population increased around 75,000–5000 years ago. The Bayesian skyline plots from the European population and World population indicated a marked Ne increase around 25,000–5000 years ago, consistent with a more recent and faster population expansion (Fig. 3d, e).
Our results showed no structured genetic patterns in H. canis populations, suggesting high levels of gene flow. Likewise, our results provide evidence of past and recent population expansions that could be due to human migrations carrying their cats and dogs as pets. These results are discussed in detail below.
Since its first description in 1905 , H. canis has been reported in several countries of America, Europe, Asia, Africa, and Australia, [1, 15, 52,53,54]. Most studies lack genetic information , reporting the presence of H. canis by histopathological or microscopic diagnosis , which complicates the better understanding of the genetic structure and gene flow of their populations. Our survey of 18S rRNA gene sequence data in populations of H. canis identified four haplogroups in the world (Fig. 1). However, haplogroup structure was not geographically congruent nor did it fit previous hypotheses proposed in the past [25, 26]. The more frequently found haplotypes (H1, H4, H7, H9, H16, H23, H28, and H41) are distributed in multiple populations in the four continents, and they come from different hosts. This differs from other studies performed with parasites like Corynosoma australe  (cytochrome c oxidase subunit I (cox 1) sequences), Angiostrongylus cantonensis and A. malaysiensis  (partial 66-kDa protein gene sequence), or Plasmodium knowlesi  (cox-1 and 18S RNA sequences), which showed haplogroups defined in the network and with geographic congruence. However, these studies used DNA fragments shorter than the data used in our present study. On the other hand, our results agreed with a previous study that reported that dogs carried the same or very similar H. canis haplotypes as red foxes, suggesting that both hosts participate in the same epidemiological cycle . In contrast, in this study, the similarities in the haplotypes between host species suggest cross-species transmission, which could be due to (i) the mobility of diverse tick species and their host across the world and (ii) shared habitats of the different host species, both promoting greater gene flow between populations (Fig. 1).
Haplotype diversity (h), a measure of species evenness, was high for all populations (Table 4), which may indicate a sustained transmission of H. canis in the world over long time periods and high gene flow. Similar high haplotype diversity values have been reported for other parasites such as P. knowlesi , widely distributed in Southeast Asia. Nucleotide diversity was low in our study, irrespective of the samples geographic or host origins, suggesting that only minor differences occurred between the haplotypes observed, as has been suggested for other protozoa [57, 58]. The haplotype diversity and distribution could be related to multiples factors associated with H. canis life cycle, transmission cycle, and dispersion capacity. Our results showed that genetic diversity across all populations was higher than within-populations and genetic differentiation analyses indicated no phylogeographical structuring. These results could be due to the existence of strong gene flow between populations as shown in other species (Cryptosporidium, Giardia, Fasciola) [59, 60]. Inferred patterns of genetic differentiation in H. canis do not appear to be correlated with the geographical regions (Fig. 1) and could be associated to diverse human migration events (travel and trade) as reported with other dog pathogens [24, 60].
FST pairwise comparisons, AMOVAs, and BAPS (Tables 2, 3, Fig. 1C) suggest low levels of genetic structuring and the existence of four genetic groups. Altogether these results suggest that dispersal of H. canis between populations is high. This pattern is commonly reported for other taxa (e.g. Psittacanthus calyculatus, canine distemper virus, Corynosoma austral, Fasciola hepatica) [24,25,26, 59]. Our results showed low structure and high genetic flow among domestic and wildlife hosts, and among populations worldwide. However, more studies are required to predict the most probable reservoirs for this parasite.
Human activities have affected the geographical ranges of many host species, shrinking by loss of habitat or increasing through introduction of individuals into new regions and promoted spillover events. In parasitic nematodes, the impact of human activity is well documented, offering an opportunity to study how changes in host population size and connectivity shape the population genetics structure . In this study, we infer that the geographical structure, genetic diversity, and relationships among populations in H. canis are affected by anthropogenic activities.
Our results showed that H. canis populations had a recent expansion (Table 4, Fig. 3). A signature of population expansion was also evident from the multi modal shape of the pairwise mismatch distribution of the 18S rRNA genes (Fig. 2) and Bayesian skyline plots (Fig. 3). Despite our results confirming the hypothesis of population expansion, these findings cannot define the geographical origin of this expansion. A further analysis dating the evolutionary divergence within Apicomplexa phylum could provide a refined overview regarding the origin of H. canis. Historical demographic expansions were determined by analyzing the frequency distributions of pairwise differences between sequences [36, 40]. Neutrality tests with Tajima’s D and Fu’s Fs statistics estimate the deviation from neutrality, based on the expectation of a constant population size at mutation-drift equilibrium. Here, a negative Tajima’s D signifies an excess of low frequency polymorphisms relative to expectation, indicating population size expansion or positive selection . Hepatozoon canis populations displayed a genetic pattern typical of a population that has undergone a recent expansion, and the range of expansion could be a recent phenomenon (5000 years ago), at least that shows the Bayesian skyline plots analysis for African population and when all population were analyzed like a group (Fig. 3). These populations may not have achieved the migration-drift equilibrium, as shown by the lack of phylogeographical structure [39, 40, 62]. The discrepancy in detecting population expansion, based on the Tajima’s D, Fu’s FS, and R2 statistics, reflects different responses to past changes (e.g. population reduction, population subdivision, a recent bottleneck) or high gene flow between populations . However, these changes could be recent in evolutionary times and could be related to anthropogenic factors, as reported in other studies carried out with canid pathogens . Moreover, the ability of H. canis to spillover among host species is feasible given the wide range of carnivores that can be infected by these hemoparasites worldwide [1, 16, 18,19,20].
In the mismatch distribution result, H. canis presented a multimodal pattern suggesting recent expansion. Furthermore, non-significant values for SSD mean that the observed data do not deviate from that expected under the model of expansion, and the non-significant raggedness index also indicates population expansion. These results suggest that population expansion occurred recently, except to Africa . Demographic expansion supported by the statistics of demographic analysis and the mismatch distribution as well as the pattern shown by the haplotype network reflects the existence of many different and rare haplotypes (high value of h, Table 4, Fig. 1), possibly because not enough time elapsed to accumulate genetic differences. Metazoonosis, like hepatozoonosis, involves complex interactions among the hosts, parasite, and environmental factors . Owing to the historical relationship between dogs and humans (since domestication in Siberia ∼ 23,000 years ago) , high mobility from the H. canis host favors dispersion and gene flow among populations. The latter is supported by several studies done with dogs and wild canids worldwide, which observed the interactions networks among parasites, vectors, and hosts [1, 20, 52, 54].
The detection of H. canis in felids and capybaras highlights the lack of information on the transmission dynamics of H. canis into new hosts [6, 8, 66]. Moreover, it shows how disregarded the role of other little-studied tick species, such as A. mixum, as potential vectors has been . Although these isolates in non-canid hosts represent incidental findings, this supports the hypothesis that H. canis moves between different host geographic areas and vectors favoring gene flow.
This study provides novel and significant insights into the phylogeography and population differentiation of H. canis in the world. The current molecular data confirmed that H. canis does not show phylogeographic and population structure, which could be due to its wide range of definitive host (domestics and wild canids), lack of clear understanding of the role of other non-canid species in the transmission dynamics, and the impact of human activities. Our results showed identical H. canis haplotypes co-occurring in several geographical regions and host species, indicating wide distribution of the parasite. The haplotypes exhibit low genetic differentiation, suggesting a recent expansion due to gene flow among populations. The reasons are mainly still unknown. However, we consider that the role of human migrations might cause the dispersion of vector and host. Finally, using more variable target regions such as mitochondrial and apicoplast genomes might improve resolution. We suggest these more variable regions should be applied in future studies to understand the natural history of H. canis and approach it from an eco-epidemiological vision.
Availability of data and materials
Data supporting the conclusions of this article are included within the article and its additional files.
Jarquín-Díaz VH, Barbachano-Guerrero A, Maldonado-Rodríguez R, Vásquez-Aguilar AA, Aguilar-Faisal JL. First molecular evidence of Hepatozoon canis in domestic dogs and ticks in fragmented rainforest areas in Mexico. Vet Parasitol Reg Stud Rep. 2016;6:4–8. https://doi.org/10.1016/j.vprsr.2016.11.001.
Perles L, Roque AL, D’Andrea PS, Lemos ER, Santos AF, Morales AC, Machado RZ, André MR. Genetic diversity of Hepatozoon spp. in rodents from Brazil. Sci Rep. 2019;9:1–9. https://doi.org/10.1038/s41598-019-46662-2.
Bouer A, André MR, Gonçalves LR, Luzzi MD, Oliveira JP, Rodrigues AC, Varani AD, Miranda VF, Perles L, Werther K, Machado RZ. Hepatozoon caimani in Caiman crocodilus yacare (Crocodylia, Alligatoridae) from North Pantanal, Brazil. Rev Bras Parasitol Vet. 2017;26(3):352–8. https://doi.org/10.1590/s1984-29612017041.
Harris DJ, Damas-Moreira I, Maia JP, Perera A. First report of Hepatozoon (Apicomplexa: Adeleorina) in caecilians, with description of a new species. J Parasitol. 2014;100(1):117–20. https://doi.org/10.1645/13-203.1.
Merino S, Martínez J, Masello JF, Bedolla Y, Quillfeldt P. First molecular characterization of a Hepatozoon species (Apicomplexa: Hepatozoidae) infecting birds and description of a new species infecting storm petrels (Aves: Hydrobatidae). J Parasitol. 2014;100(3):338–43. https://doi.org/10.1645/13-325.1.
Baneth G. Perspectives on canine and feline hepatozoonosis. Vet Parasitol. 2011;181(1):3–11. https://doi.org/10.1016/j.vetpar.2011.04.015.
Smith TG. The genus Hepatozoon (Apicomplexa: Adeleina). J Parasitol. 1996;82(4):565–85. https://doi.org/10.2307/3283781.
Schäfer I, Kohn B, Volkmann M, Müller E. Retrospective evaluation of vector-borne pathogens in cats living in Germany (2012–2020). Parasit Vectors. 2021;14(1):1–9. https://doi.org/10.1186/s13071-021-04628-2.
Baneth GA, Samish M, Alekseev E, Aroch I, Shkap V. Transmission of Hepatozoon canis to dogs by naturally-fed or percutaneously-injected Rhipicephalus sanguineus ticks. J Parasitol. 2001;87(3):606–11. https://doi.org/10.1645/0022-3395(2001)087[0606:TOHCTD]2.0.CO;2.
Otranto D, Dantas-Torres F, Weigl S, Latrofa MS, Stanneck D, Decaprariis D, Capelli G, Baneth G. Diagnosis of Hepatozoon canis in young dogs by cytology and PCR. Parasit Vectors. 2011;4(1):1–6. https://doi.org/10.1186/1756-3305-4-55.
Allen KE, Li Y, Kaltenboeck B, Johnson EM, Reichard MV, Panciera RJ, Little SE. Diversity of Hepatozoon species in naturally infected dogs in the southern United States. Vet Parasitol. 2008;154(3–4):220–5. https://doi.org/10.1016/j.vetpar.2008.03.027.
Murata T, Inoue M, Taura Y, Nakayama S, Abe H, Fujisaki K. Detection of Hepatozoon canis oocyst from ticks collected from the infected dogs. J Vet Med Sci/Jpn Soc Vet Sci. 1995;57(1):111–2. https://doi.org/10.1292/jvms.57.111.
Forlano M, Scofield A, Elisei C, Fernandes KR, Ewing SA, Massard CL. Diagnosis of Hepatozoon spp. in Amblyomma ovale and its experimental transmission in domestic dogs in Brazil. Vet Parasitol. 2005;134(1–2):1–7. https://doi.org/10.1016/j.vetpar.2005.05.066.
Giannelli A, Lia RP, Annoscia G, Buonavoglia C, Lorusso E, Dantas-Torres F, Baneth G, Otranto DA. Rhipicephalus turanicus, a new vector of Hepatozoon canis. Parasitology. 2017;144(6):730–7. https://doi.org/10.1017/S003118201600250X.
Gabrielli S, Kumlien S, Calderini P, Brozzi A, Iori A, Cancrini G. The first report of Hepatozoon canis identified in Vulpes vulpes and ticks from Italy. Vector Borne Zoonotic Dis. 2010;10(9):855–9. https://doi.org/10.1089/vbz.2009.0182.
Najm NA, Meyer-Kayser E, Hoffmann L, Pfister K, Silaghi C. Hepatozoon canis in German red foxes (Vulpes vulpes) and their ticks: molecular characterization and the phylogenetic relationship to other Hepatozoon spp. Parasitol Res. 2014;113(7):2679–85. https://doi.org/10.1007/s00436-014-3923-8.
Díaz-Sánchez AA, Hofmann-Lehmann R, Meli ML, Roblejo-Arias L, Fonseca-Rodríguez O, Castillo AP, Cañizares EV, Rivero EL, Chilton NB, Corona-González B. Molecular detection and characterization of Hepatozoon canis in stray dogs from Cuba. Parasitol. 2021;80: 102200. https://doi.org/10.1016/j.parint.2020.102200.
Farkas R, Solymosi N, Takács N, Hornyák Á, Hornok S, Nachum-Biala Y, Baneth G. First molecular evidence of Hepatozoon canis infection in red foxes and golden jackals from Hungary. Parasit Vectors. 2014;7(1):1–7. https://doi.org/10.1186/1756-3305-7-303.
da Silva MR, Fornazari F, de Castro DL, Teixeira CR, Langoni H, O’Dwyer LH. Didelphis albiventris naturally infected with Hepatozoon canis in southeastern Brazil. Ticks Tick Borne Dis. 2017;8(6):878–81. https://doi.org/10.1016/j.ttbdis.2017.07.005.
Giannelli A, Latrofa MS, Nachum-Biala Y, Hodžić A, Greco G, Attanasi A, Annoscia G, Otranto D, Baneth G. Three different Hepatozoon species in domestic cats from southern Italy. Ticks Tick Borne Dis. 2017;8(5):721–4. https://doi.org/10.1016/j.ttbdis.2017.05.005.
Bisharat N, Koton Y, Oliver JD. Phylogeography of the marine pathogen, Vibrio vulnificus, revealed the ancestral scenarios of its evolution. Microbiol Open. 2020;9(9): e1103. https://doi.org/10.1002/mbo3.1103.
Ornelas JF, Licona-Vera Y, Vásquez-Aguilar AA. Genetic differentiation and fragmentation in response to climate change of the narrow endemic Psittacanthus auriculatus. Trop Conserv Sci. 2018;11:1940082918755513. https://doi.org/10.1177/1940082918755513.
Pérez-Crespo MJ, Ornelas JF, González-Rodríguez A, Ruiz-Sanchez E, Vásquez-Aguilar AA, Ramírez-Barahona S. Phylogeography and population differentiation in the Psittacanthus calyculatus (Loranthaceae) mistletoe: a complex scenario of climate–volcanism interaction along the Trans-Mexican Volcanic Belt. J Biogeogr. 2017;44(11):2501–14. https://doi.org/10.1111/jbi.13070.
Panzera Y, Sarute N, Iraola G, Hernández M, Pérez R. Molecular phylogeography of canine distemper virus: geographic origin and global spreading. Mol Phylogenet Evol. 2015;92:147. https://doi.org/10.1016/j.ympev.2015.06.015.
García-Varela M, Masper A, Crespo EA, Hernández-Orts JS. Genetic diversity and phylogeography of Corynosoma australe Johnston, 1937 (Acanthocephala: Polymorphidae), an endoparasite of otariids from the Americas in the northern and southern hemispheres. Parasitol Int. 2021;80: 102205. https://doi.org/10.1016/j.parint.2020.102205.
Eamsobhana P, Yong HS, Song SL, Gan XX, Prasartvit A, Tungtrongchitr A. Molecular phylogeography and genetic diversity of Angiostrongylus cantonensis and A. malaysiensis (Nematoda: Angiostrongylidae) based on 66-kDa protein gene. Parasitol Int. 2019;68(1):24–30. https://doi.org/10.1016/j.parint.2018.09.006.
Yusof R, Ahmed MA, Jelip J, Ngian HU, Mustakim S, Hussin HM, Fong MY, Mahmud R, Sitam FA, Japning JR, Snounou G. Phylogeographic evidence for 2 genetically distinct zoonotic Plasmodium knowlesi parasites, Malaysia. Emerg Infect Dis. 2016;22(8):1371–80. https://doi.org/10.3201/eid2208.151885.
Bertranpetit E, Jombart T, Paradis E, Pena H, Dubey J, Su C, Mercier A, Devillard S, Ajzenberg D. Phylogeography of Toxoplasma gondii points to a South American origin. Infect Genet Evol. 2017;1(48):150–5. https://doi.org/10.1016/j.meegid.2016.12.020.
R Development Core Team. R: a language and environment for statistical computing. Vienna: R Foundation for Statistical Computing; 2017.
Larkin MA, Blackshields G, Brown NP, Chenna R, McGettigan PA, McWilliam H, Valentin F, Wallace IM, Wilm A, Lopez R, Thompson JD. Clustal W and clustal X version 2.0. BIX. 2007;23(21):2947–8. https://doi.org/10.1093/bioinformatics/btm404.
Müller K, Müller J, Neinhuis C, Quandt D. PhyDE–phylogenetic data editor, v0. 995. Program distributed by the authors; 2006. http://www.phyde.de/download.html. Accessed 19 Nov 2020.
Clement M, Posada DC, Crandall KA. TCS: a computer program to estimate gene genealogies. Mol Ecol. 2000;9(10):1657–9. https://doi.org/10.1046/j.1365-294x.2000.01020.x.
Pfenninger M, Posada D. Phylogeographic history of the land snail Candidula unifasciata (Helicellinae, Stylommatophora): fragmentation, corridor migration, and secondary contact. Evolution. 2002;56(9):1776–88. https://doi.org/10.1111/j.0014-3820.2002.tb00191.x.
Corander J, Sirén J, Arjas E. Bayesian spatial modeling of genetic population structure. Comput Stat. 2008;23:111–29. https://doi.org/10.1007/s00180-007-0072-x.
Pons O, Petit RJ. Measwring and testing genetic differentiation with ordered versus unordered alleles. Genetics. 1996;144(3):1237–45. https://doi.org/10.1093/genetics/144.3.1237.
Excoffier L, Laval G, Schneider S. Arlequin (version 3.0): an integrated software package for population genetics data analysis. Evol Bioinform. 2005;1:47–50. https://doi.org/10.1177/117693430500100003.
Excoffier L, Smouse PE, Quattro JM. Analysis of molecular variance inferred from metric distances among DNA haplotypes: application to human mitochondrial DNA restriction data. Genetics. 1992;131(2):479–91.
Dupanloup I, Schneider S, Excoffier L. A simulated annealing approach to define the genetic structure of populations. Mol Ecol. 2002;11(12):2571–81. https://doi.org/10.1046/j.1365-294X.2002.01650.x.
Tajima F. Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics. 1989;123(3):585–95.
Fu YX. Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. Genetics. 1997;147(2):915–25. https://doi.org/10.1093/genetics/147.2.915.
Harpending HC. Signature of ancient population growth in a low-resolution mitochondrial DNA mismatch distribution. Hum Biol. 1994;60:591–600.
Schneider S, Excoffier L. Estimation of past demographic parameters from the distribution of pairwise differences when the mutation rates vary among sites: application to human mitochondrial DNA. Genetics. 1999;152(3):1079–89. https://doi.org/10.1093/genetics/152.3.1079.
Rogers AR, Harpending H. Population growth makes waves in the distribution of pairwise genetic differences. Mol Biol Evol. 1992;9(3):552–69. https://doi.org/10.1093/oxfordjournals.molbev.a040727.
Rozas J, Ferrer-Mata A, Sánchez-DelBarrio JC, Guirao-Rico S, Librado P, Ramos-Onsins SE, Sánchez-Gracia A. DnaSP 6: DNA sequence polymorphism analysis of large data sets. Mol Biol Evol. 2017;34(12):3299–302. https://doi.org/10.1093/molbev/msx248.
Ramos-Onsins SE, Rozas J. Statistical properties of new neutrality tests against population growth. Mol Biol Evol. 2002;19(12):2092–100. https://doi.org/10.1093/oxfordjournals.molbev.a004034.
Drummond AJ, Rambaut A, Shapiro BE, Pybus OG. Bayesian coalescent inference of past population dynamics from molecular sequences. Mol Biol Evol. 2005;22(5):1185–92. https://doi.org/10.1093/molbev/msi103.
Bouckaert R, Heled J, Kühnert D, Vaughan T, Wu CH, Xie D, Suchard MA, Rambaut A, Drummond AJ. BEAST 2: a software platform for Bayesian evolutionary analysis. PLoS Comput Biol. 2014;10(4):e1003537. https://doi.org/10.1371/journal.pcbi.1003537.
Darriba D, Taboada GL, Doallo R, Posada D. jModelTest 2: more models, new heuristics and parallel computing. Nat methods. 2012;9(8):772–772. https://doi.org/10.1038/nmeth.2109.
Morrison DA. How old are the extant lineages of Toxoplasma gondii? Parassitologia. 2005;47(2):205–14.
Rambaut A, Drummond AJ, Xie D, Baele G, Suchard MA. Posterior summarization in Bayesian phylogenetics using Tracer 1.7. Syst Biol. 2018;67(5):901. https://doi.org/10.1093/sysbio/syy032.
James SP. A new Leucocytozoon of dogs. BMJ. 1905;1(2320):1361.
Guo WP, Xie GC, Xue ZQ, Yu JJ, Jian R, Du LY, Li YN. Molecular detection of Hepatozoon canis in dogs and ticks in Shaanxi province, China. Comp Immunol Microbiol Infect Dis. 2020;72: 101514. https://doi.org/10.1016/j.cimid.2020.101514.
Penzhorn BL, Netherlands EC, Cook CA, Smit NJ, Vorster I, Harrison-White RF, Oosthuizen MC. Occurrence of Hepatozoon canis (Adeleorina: Hepatozoidae) and Anaplasma spp. (Rickettsiales: Anaplasmataceae) in black-backed jackals (Canis mesomelas) in South Africa. Parasite Vectors. 2018;11(1):1–7. https://doi.org/10.1186/s13071-018-2714-y.
Greay TL, Barbosa AD, Rees RL, Paparini A, Ryan UM, Oskam CL, Irwin PJ. An Australian dog diagnosed with an exotic tick-borne infection: should Australia still be considered free from Hepatozoon canis? Int J Parasitol. 2018;48(11):805–15. https://doi.org/10.1016/j.ijpara.2018.05.002.
Kwon SJ, Kim YH, Oh HH, Choi US. First case of canine infection with Hepatozoon canis (Apicomplexa: Haemogregarinidae) in the Republic of Korea. Korean J Parasitol. 2017;55(5):561–4. https://doi.org/10.3347/kjp.2017.55.5.561.
Helm CS, von Samson-Himmelstjerna G, Liesner JM, Kohn B, Müller E, Schaper R, Pachnicke S, Schulze C, Krücken J. Identical 18S rRNA haplotypes of Hepatozoon canis in dogs and foxes in Brandenburg, Germany. Ticks Tick Borne Dis. 2020;11(6): 101520. https://doi.org/10.1016/j.ttbdis.2020.101520.
Taylor JE, Pacheco MA, Bacon DJ, Beg MA, Machado RL, Fairhurst RM, Herrera S, Kim JY, Menard D, Póvoa MM, Villegas L. The evolutionary history of Plasmodium vivax as inferred from mitochondrial genomes: parasite genetic diversity in the Americas. Mol Biol Evol. 2013;30(9):2050–64. https://doi.org/10.1093/molbev/mst104.
Spotin A, Moslemzadeh HR, Mahami-Oskouei M, Ahmadpour E, Niyyati M, Hejazi SH, Memari F, Noori J. Phylogeography, genetic variability and structure of Acanthamoeba metapopulations in Iran inferred by 18S ribosomal RNA sequences: a systematic review and meta-analysis. Asian Pac J Trop Med. 2017;10(9):855–63. https://doi.org/10.1016/j.apjtm.2017.08.011.
Beesley NJ, Williams DJ, Paterson S, Hodgkinson J. Fasciola hepatica demonstrates high levels of genetic diversity, a lack of population structure and high gene flow: possible implications for drug resistance. Int J Parasitol. 2017;47(1):11–20. https://doi.org/10.1016/j.ijpara.2016.09.007.
Garcia-R JC, French N, Pita A, Velathanthiri N, Shrestha R, Hayman D. Local and global genetic diversity of protozoan parasites: spatial distribution of Cryptosporidium and Giardia genotypes. PLoS Negl Trop Dis. 2017;11(7):e0005736. https://doi.org/10.1371/journal.pntd.0005736.
Cole R, Viney M. The population genetics of parasitic nematodes of wild animals. Parasite Vectors. 2018;11(1):1–20. https://doi.org/10.1186/s13071-018-3137-5.
Tajima F. Evolutionary relationship of DNA sequences in finite populations. Genetics. 1983;105(2):437–60. https://doi.org/10.1093/genetics/105.2.437.
Ornelas JF, García JM, Ortiz-Rodriguez AE, Licona-Vera Y, Gándara E, Molina-Freaner F, Vásquez-Aguilar AA. Tracking host trees: the phylogeography of endemic Psittacanthus sonorae (Loranthaceae) mistletoe in the Sonoran Desert. J Heredity. 2019;110(2):229–46. https://doi.org/10.1093/jhered/esy065.
Giraudoux P, Raoul F, Pleydell D, Craig PS. Multidisciplinary studies, systems approaches and parasite eco-epidemiology: something old, something new. Parasite. 2008;15(3):469–76. https://doi.org/10.1051/parasite/2008153469.
Perri AR, Feuerborn TR, Frantz LA, Larson G, Malhi RS, Meltzer DJ, Witt KE. Dog domestication and the dual dispersal of people and dogs into the Americas. Proc Natl Acad Sci. 2021;118(6):1–8. https://doi.org/10.1073/pnas.2010083118.
Gomes LD, Moraes LA, Aguiar DC, Dias HL, Ribeiro AS, Rocha HP, Nunes MR, Gonçalves EC. Genetic diversity of Hepatozoon spp in Hydrochoerus hydrochaeris and Pecari tajacu from eastern Amazon. Ticks Tick Borne Dis. 2018. https://doi.org/10.1016/j.ttbdis.2017.11.005.
We thank all the contributors who have submitted sequence data to GenBank over the years. Our special thanks go to D. Hernandez-Rodriguez for valuable criticisms and suggestions and E. A. Lopez-Huicochea for the support and guidance with the BAPS analysis.
Open Access funding enabled and organized by Projekt DEAL. This work was supported by MDC Berlin (Helmholtz Consortium).
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
GenBank accession numbers for the sequences used on this study.
Statistical parsimony network of 76 rRNA 18S haplotypes found in 46 populations of H. canis worldwide. Black dots represent unsampled haplotypes. The frequency of each haplotype is represented by the size of the circle. Each haplotype is coded with a different color according to the continent. The numbers inside the haplotypes in the network indicate the number of individuals that share that haplotype. The four rectangles represent the four haplogroups reported in this study. Figure S2. Geographical distribution of host and statistical parsimony network of 76 rRNA 18S haplotypes found of H. canis in the world. a Distribution map of the hosts. Black dots represent sampled populations and pie charts represent hosts found in each sampling population. Section size of pie charts corresponds to the proportion of individuals with given hosts. b Haplotype network. Black dots represent unsampled haplotypes. The frequency of each haplotype is represented by the size of the circle. Each haplotype is coded with a different color according to the hosts. The numbers inside the haplotypes in the network indicate the number of individuals that share that haplotype. Table S1. Number of analyzed samples (n) for molecular marker (18S rRNA gene) and hosts found of H. canis.
About this article
Cite this article
Vásquez-Aguilar, A.A., Barbachano-Guerrero, A., Angulo, D.F. et al. Phylogeography and population differentiation in Hepatozoon canis (Apicomplexa: Hepatozoidae) reveal expansion and gene flow in world populations. Parasites Vectors 14, 467 (2021). https://doi.org/10.1186/s13071-021-04924-x
- Dog parasites
- Tick-borne pathogens
- 18S rRNA gene