No CrossRef data available.
Published online by Cambridge University Press: 17 October 2003
Babesia microti (Apicomplexa: Piroplasmida) has historically been considered a common parasite of Holarctic rodents. However, human babesiosis due to this species has generally been limited to the northeastern seaboard of the United States and Minnesota and Wisconsin. The absence of reports of B. microti babesiosis from sites where the agent is enzootic, such as in western Europe, remains unexplained. Previous work focusing on the 18S rDNA demonstrates little sequence diversity among samples from allopatric host populations across a wide geographical area. It may be that genetic diversity is underestimated due to sample size or the gene analysed. Accordingly, we collected blood or spleen samples from American or Eurasian animals with parasites that were morphologically consistent with B. microti, amplified the 18S rDNA and beta-tubulin gene, and conducted phylogenetic analysis. Surprisingly, what was considered to be ‘B. microti’ by microscopy appears to be a diverse species complex. We identify 3 distinct clades within this complex, including parasites from non-rodent hosts. Rodent parasites comprise 2 clades, one representing zoonotic isolates, and the other apparently maintained in microtine rodents, and therefore their morphological detection within animals from a site does not necessarily imply a risk to public health.
Human babesiosis due to Babesia microti (Apicomplexa: Piroplasmida) has mainly been reported from coastal New England and upper Midwestern sites in the US. Hundreds of cases have been reported from immunocompetent persons, with a 5% case fatality rate (Telford & Maguire, 1999). Although B. microti is considered to have a Holarctic distribution, few cases outside the US have been reported (van Peenen et al. 1977; Shih et al. 1997; Tsuji et al. 2001). European sites, in particular, would seem to be likely venues of transmission, because B. microti-infected rodents have frequently been found throughout Europe (Franca, 1910; Krampitz, 1979; Hussein, 1980; Healing, 1981; Telford & Spielman, 1993; Karbowiak et al. 1999; Bajer et al. 2001). The absence of reports of B. microti babesiosis in Europe has been suggested to be the consequence of the great degree of host specificity of the mouse tick, Ixodes trianguliceps, which appears to be the main enzootic vector (Hussein, 1980; Randolph, 1994; Young, 1970). However, the European vector of Lyme borreliosis, I. ricinus, is sympatric in many of these sites, feeds upon small rodents that would be infected by B. microti, and is capable of transmitting the piroplasm in the laboratory (Gray et al. 2002). Recent reports suggest that B. microti DNA may be detected within host-seeking I. ricinus (Duh, Petrovec & Avsic-Zupanc, 2001; Foppa et al. 2002), and that serological evidence of exposure to B. microti may be detected in tick-exposed people (Hunfeld et al. 2002). Although less than optimal physician awareness may partially explain the absence of recognized B. microti cases, we note that a different babesiosis is well-known by European physicians: B. divergens babesiosis in immunocompromised individuals usually constitutes a life-threatening situation, and is a differential diagnosis for acute febrile haemolytic syndromes (Gorenflot et al. 1998). These considerations make the absence of reports of clinical babesiosis due to B. microti in Europe even more enigmatic.
It may be that B. microti strains circulating in rodents vary in their capacity to cause pathology in humans. A recent report of B. microti babesiosis in a Japanese resident suggested that the parasite that caused human infection differed from the strains that were commonly found in rodents throughout Japan, and was found only in a small insular population of rodents (Tsuji et al. 2001). Alternatively, because B. microti has historically been identified by means of light microscopy of Giemsa-stained blood smears, reports of this species from diverse sites and even hosts such as shrews (Franca, 1910; Telford et al. 1990), cottontail rabbits (Spielman et al. 1981), macaques (as Entopolypoides macaci (Hawking, 1972)) or dogs (Camacho et al. 2001; Zahler et al. 2000 b) may not truly reflect parasite identity or diversity.
The pan-Holarctic distribution of parasites identified by microscopy as B. microti would suggest that molecular phylogenetic methods should detect significant variation due to typical processes of genetic divergence between allopatric populations separated by long spans of time. Recent DNA sequence analyses of B. microti, however, have yielded some unexpected results. Minor, if any, DNA sequence divergence was demonstrated between samples collected in Switzerland, Russia and the northeastern US (Duh et al. 2001; Foppa et al. 2002; Telford et al. 2002). But, only small portions of the 18S rDNA gene (approx 500 bp) were used in these analyses. The 18S rDNA gene is highly conserved over a wide range of taxa and is most often used for interspecific comparisons. Because mutations accumulate so slowly in this gene, large portions must be sequenced to accumulate ample signal to distinguish closely related species. On the other hand, we recently reported the presence of B. microti in Maine that, although found near endemic areas for human babesiosis, differs in its 18S rDNA sequence from human-derived isolates; the Maine sample is more similar to those that have been classified as B. microti from Alaskan microtine rodents (Goethert et al. 2003). To date, no study has attempted to analyse genetic diversity comprehensively among parasites that are morphologically identifiable as B. microti. Accordingly, to determine whether B. microti may genetically differ by site and host, we collected blood and spleen samples from animals with microscopy-based evidence of infection. DNA was extracted, and a major portion of the 18S rDNA amplified and sequenced. In addition, because the beta-tubulin gene has been successfully used for population-level comparisons of other parasitic protozoa (Caccio et al. 1999, 2000), we amplified and sequenced this target as a comparison. The resulting sequences were aligned and analysed by multiple phylogenetic algorithms.
Animals or ticks were generally collected during the course of epizootiological studies of Lyme disease over a number of years or as part of the Beringian Coevolution Project (Brooks & Hoberg, 2000; Conroy & Cook, 1999; Hoberg et al. 1999). Sites sampled include: Mis, Russia; Cordova, Alaska; Big Sky, Montana; Spooner, Wisconsin; Nantucket, Ipswich, Martha's Vineyard, Naushon and Cape Cod, Massachusetts; Reutiwis, Switzerland; Vinalhaven and Grand Lake Stream, Maine; Old Lyme, Connecticut; and Galicia, Spain. (Table 1). Rodents were snap trapped or live trapped; medium-sized mammals were captured in Tomahawk live traps. Citrate-anticoagulated blood and spleens were obtained from animals at most sites, but in a few instances, only blood (medium sized mammals) or spleens were sampled. The canine blood sample was obtained from a medical-veterinary diagnostic laboratory. All of the analyses were performed directly on field-collected materials, thereby avoiding the risk of contamination by serial blood passage of strains within the laboratory. Sequences from other taxa that were not directly sampled were from GenBank. Blood smears were obtained from most samples. These were stained with Giemsa's stain and examined under oil immersion.
Blood and spleens were extracted using either the Isoquick kit (Orca Research) or DNeasy Tissue kit (Qiagen Inc) following the instructions of the manufacturer. All samples were resuspended in 50 μl final volume. PCR was conducted using Taq polymerase (Qiagen) following the instructions of the manufacturer. Briefly, 12·5 pmol of each primer, 3 mM dUTP, 25 mM dNTPs, 1× Taq Buffer, 0·25 μl Taq polymerase and 2·5 μl of template were used in 25 μl reactions. An approximately 1300 bp piece of the 18S rDNA gene was amplified using primer PiroA (Armstrong et al. 1998) and primer B as described previously (Kjemtrup et al. 2000). Primers Babtub1 (GCCAAATTCTGGGAAGTCATCTCWGATGARCAYGG) and Babtub3 (TCACAAACACTTGACTTTGTGTTGTGTGGDATYCAYTC), that amplified 900 bp of the beta-tubulin gene were inferred from protein alignments of B. bovis with other apicomplexan parasites using CODEHOP (Rose et al. 1998). The amplification parameters were as follows: 94 °C 1 min, then 40 cycles of 94 °C 1 min, 63 °C 1 min and 72 °C 1·5 min with a step-down 1 degree each cycle until 55 °C. Final extension was 5 min at 72 °C. Once the initial sequencing was done, more specific primers were designed, BtubmicF (SGGTATWGATCCRGTAWGTA) and BtubmicR (ACRAAGAAACTKGARTTCTTAT). We practiced ‘safe PCR’ utilizing barrier pipette tips, dedicated rooms for pre- and post-PCR, and always included negative control (no template) samples with every set of reactions. Amplicons were visualized with agarose gel electrophoresis. Bands were excised and purified using Qiaquick gel purification columns. Purified DNA was quantified by running a small amount on a gel next to a low molecular weight ladder sample (Gibco) and sent to the University of Maine Sequencing Facility for sequence analysis (both strands were sequenced). Amplicons from multiple animals from each field site, if available, were analysed separately to test for consistency in the resulting sequence.
Sequences were aligned using ClustalX (BCM Search Launcher) and then adjusted by eye using GeneDoc (Nicholas & Nicholas, 1997). Maximum parsimony, neighbour-joining, and maximum likelihood methods were used for analysis (Swofford, 1998). For the 18S rDNA analysis a 1249 base-pair piece was used. Areas of ambiguity in the alignment were excluded from the analysis, 33 bps in total. Maximum parsimony analysis compared trees using equally weighted bases with a transition/transversion ratio of 2. Heuristic searches were conducted with TBR swapping. Bootstrap replicates, 500, were used to assess the robustness of the resulting topology. Neighbour-joining analysis was performed using the Kimura 2-parameter model so that we could compare our results with previously reported phylogenetic analyses of piroplasms (Allsopp et al. 1994; Kjemtrup et al. 2000; Zahler, Rinder & Gothe, 2000 a). For maximum likelihood analysis, the statistically best model for our data was determined using Modeltest (Posada & Crandall, 1998). This model, GTR+G+I (G=shape parameter of the gamma distribution; I=proportion of invariable sites) using parameters estimated from the data was then used for subsequent analysis. Analysis of the beta-tubulin gene used a 904 bp piece. Codon positions were assigned and verified by translation; the correct highly conserved protein sequence for this gene was obtained. Maximum parsimony, neighbour-joining, and maximum likelihood analyses (model GTR+G) were conducted in a manner similar to those for the 18S gene. Trees were evaluated with and without the highly variable 3rd codon. B. divergens and B. odocoilei were used as outgroups for all analyses, except for one neighbour-joining analysis that used Toxoplasma gondii and Crypthecodinium cohnii. Sequence divergence among clades and polymorphism within clades was calculated with the Jukes-Cantor correction using Dnasp (Rozas & Rozas, 2000).
GenBank Accession numbers: ‘T. annae’ AF188001; B. rodhaini AB049999; B. divergens U16370; B. odocoilei U16369; WA1 AF158700; B. gibsoni 1 AF158702; T. youngi AF245279; T. buffeli Z15106; T. parva L02366, B. gibsoni 2 AF175300; Toxoplasma gondii X68523; C. cohnii M64245.
All the parasites included in this analysis are morphologically similar to B. microti (Fig. 1) i.e., small ring forms with a white vacuole and extended as opposed to punctate chromatin; parasites were generally centrally located within the erythrocyte. Phylogenetic analysis of 18S rDNA sequences, by maximum parsimony or neighbour-joining methods establishes that all of our samples are closely related to B. microti, as represented by the well-characterized Harvard experimental GI strain from Nantucket Island, Massachusetts (Figs 2 and 3). All the parasites that were analysed appear to be monophyletic, with a highly statistically significant bootstrap value for this major node (97%). B. rodhaini, classically considered the most closely related species to B. microti, is excluded from this clade. Although the parasites that we considered to be B. microti by microscopy are monophyletic, within the large clade is what we conservatively interpret to be a species-complex with significant sequence diversity. The 18S gene is relatively conserved for this group of organisms. The divergence of these sequences from the GI sequence from Nantucket ranges from 0 to 5·4% (Table 2). Three distinct clades are also demonstrable by analysis of the beta-tubulin gene (Figs 2–4). Sequence divergence in the beta-tubulin gene is much greater than within the 18S rDNA, ranging from 0 to 16·2% (Table 2).
Fig. 1. Morphology of representative parasites on blood smears. (A) Vole Cordova, AK; (B) white-footed mouse Nantucket Island, MA; (C) Vole Mis, Russia; (D) skunk Martha's Vineyard, MA; (E) dog Galicia, Spain. Scale bars and numbers equal average size of host red blood cells in microns.
Fig. 2. Phylogenetic analyis of Babesia microti using the neighbour-joining algorithm. Shown are bootstrap consensus trees of (A) representative 18S rDNA sequences of Babesia and Theileria (one representative sequence from each B. microti clade is included); (B) 18S rDNA of B. microti clades; (C) beta-tubulin gene of B. microti clades with the 3rd codon included; (D) beta-tubulin gene of B. microti clades with 3rd codon excluded.
Fig. 3. Phylogenetic analysis of Babesia microti clades using maximum parsimony. Shown are bootstrap consensus trees of (A) 18S rDNA gene and (B) beta-tubulin gene with the 3rd codon included.
Fig. 4. Phylogenetic analysis of Babesia microti clades using maximum likelihood. (A) 18S rDNA gene and (B) beta-tubulin gene with the 3rd codon included.
Clade 1 includes all samples from humans as well as samples from ticks and rodents from Massachusetts and the rest of New England, Wisconsin, Switzerland and Russia. Clade 2 includes the sequences from non-rodent hosts: skunk, fox, dog and raccoon. Clade 3 includes samples from rodents from Maine, Montana and Alaska. Amplicons from a median of 4 animals (range 1 to 8) were sequenced from each location. Sequences within clade 1 are highly invariant. Both genes yielded little sequence diversity, the 18S gene showed 0·02% sequence diversity and the beta-tubulin gene showed 1·2% (Table 2). Clade 2 shows much greater intraclade diversity: 18S rDNA has 2·3% sequence diversity, whereas the beta-tubulin gene has 9·5% diversity. Although Clade 3 was less diverse in the overall sequence differences, 18S sequences were all identical and beta-tubulin showed 8·6% sequence diversity; it was the only clade in which population level polymorphism occurs i.e., more than 1 sequence was documented from different animals within a single field site (Figs 2–4). The presence of the 3 clades is supported in this analysis with the sequence diversity among the clades greater than the diversity within each clade (Table 2).
Although the existence of each clade was consistent regardless of the analytical method and these clades were always highly supported by bootstrap analysis, the relationship between the clades is uncertain. The 18S rDNA did not have enough phylogenetic signal to determine the branching pattern at the terminal nodes. The branching pattern connecting the 3 nodes was extremely variable and dependent on the method of analysis and the taxa included. We provide the maximum parsimony bootstrap consensus tree that presents the 3 clades as a polytomy (Fig. 3). The beta-tubulin gene tree corroborates the topology of the 18S rDNA tree (Figs 2–4) in that the same 3 major clades are formed. However, the beta-tubulin gene has greater sequence diversity and greater phylogenetic signal (Table 2). Most of the sequence diversity identified in the beta-tubulin gene has accumulated in the 3rd codon position. Analysis with this codon excluded still yields the same 3 clades but all sequence diversity within each clade disappears (Fig. 2).
The monophyly of this species-complex is not borne-out by maximum likelihood analysis. This analysis of the 18S gene divides the parasites analysed into 2 separate clades, one that includes all the samples from rodents (including B. rodhaini) and one that includes only samples from medium-sized mammals (Fig. 4). This branching order is supported by the beta-tubulin gene when the 3rd codon is included (Fig. 4). Unfortunately, we do not have beta-tubulin sequence from B. rodhaini so the placement of this key species cannot be assessed.
Contrary to expectations derived from morphological and host characteristics, B. microti is not a single invariant species across its distribution and among diverse hosts; instead it is a genetically diverse species complex. All of the sampled parasites are monophyletic by either parsimony or neighbour-joining analyses, which supports our suggestion that they represent a species complex, but a substantial degree of genetic divergence is apparent. Three distinct groups comprise the species complex. These groups are consistently demonstrated with high bootstrap support for the nodes regardless of the method of analysis or gene.
With parsimony or neighbour-joining, parasites from carnivores, surprisingly, clustered with those from rodents. Indeed, they did so more than did B. rodhaini, long considered the most closely related species to B. microti (Levine, 1971); indeed ATCC has listed B. rodhaini (ATCC 30222) as B. microti, and it has been designated a lethal B. microti strain in certain experimental studies (e.g. Hughes & Oz, 1995; Oz & Hughes, 1996). On the other hand, using a more objective method of relating sequence information, selecting a statistically robust model (Modeltest) a priori with subsequent maximum likelihood analysis, B. rodhaini clusters with Clades 1 and 3 to the exclusion of Clade 2, which contains the carnivore parasites. This conclusion is more biologically consistent in that Fahrenholz' Rule is not violated: phylogenetically related hosts should contain closely related parasites (Fahrenholz, 1913). However, there are many known exceptions to Fahrenholz' Rule, and without additional sequence comparisons, including other known rodent piroplasms such as Nuttallia danii of Middle Eastern gerbils (Killick-Kendrick, 1974; Tsur, Hadani & Pipano, 1960) we cannot conclude, at present, whether rodent piroplasms may be polyphyletic.
Paralogous sequences have been described for the 18S rDNA of apicomplexans (Corredor & Enea, 1994; McCutchan et al. 1995). It seems unlikely that this phenomenon accounts for the sequence diversity that we observed for this gene. The 18S rDNA sequences of B. bovis and B. bigemina contain limited sequence heterogeneity (Dalrymple, 1990; Reddy et al. 1991). The great sequence diversity that we found in our clades, such as the 38 bp difference between B. microti from Massachusetts and that from Alaska, are unlikely to be accounted for by this phenomenon. We sequenced amplicons from multiple animals sampled at each field site and found consistency in the sequence derived from each trapping location. Although the beta-tubulin gene has been noted to have polymorphic alleles (Caccio et al. 1999; Widmer et al. 1998), other apicomplexans tested revealed the presence of only one copy of the gene (Caccio, LaRosa & Pozio, 1997; Nagel & Boothroyd, 1988; Wesseling et al. 1989).
Clade 1 encompasses areas where human babesiosis is regularly reported (Massachusetts, Connecticut and Wisconsin), as well as some areas where it is rare or absent, such as VH-Maine, Russia and Switzerland. This clade has little genetic diversity for the 18S rDNA gene, with all sequences virtually identical, and the least sequence diversity for the beta-tubulin gene. Human-biting Ixodes (subgenus Ixodes) persulcatus-complex ticks – I. dammini, eastern US; I. ricinus, Switzerland; I. persulcatus, Russia – are present in these sites and, indeed, serve as the main vectors of Lyme disease. Thus, we cannot readily explain the absence of human babesiosis reports from Russia and Switzerland. Recent serosurveys in Switzerland and Germany, however, indicate that humans are exposed to infection but have not reported symptoms (Hunfeld et al. 2002), suggesting the possibility that local strains may vary in their capacity to cause overt disease, as recently suggested for Japan (Tsuji et al. 2001).
Clade 3 includes rodent parasites from Alaska, Montana and GLS-Maine, areas from which human cases have not been reported. Clade 3 may be associated with ticks of the subgenus Ixodiopsis, such as I. angustus, a widely distributed nidicolous species that very rarely bites humans. This tick was incriminated as the vector for Alaskan B. microti (Fay & Rausch, 1969), and is present in GLS-Maine (E. H. Lacombe, Maine Medical Center, personal communication). The putative vector for the Montana B. microti population may also be I. angustus, although I. spinipalpis was reported to be the vector for B. microti in Colorado voles (Burkot et al. 2000); this tick also rarely bites humans. Accordingly, the absence of reports of B. microti babesiosis in the Clade 3 sites may represent the host specificity of the vector. The capacity for Clade 3 parasites to infect and cause pathology in humans, however, remains unknown.
Clade 2 represents parasites from carnivores, and includes the recently described B. microti-like agent that causes fulminating disease in Spanish dogs, Theileria annae (Camacho et al. 2001; Zahler et al. 2000 b). We note that this parasite clusters within the B. microti species complex, and is only distantly related to true Theileria spp. such as T. parva. No evidence was presented (Zahler et al. 2000 b) for pre-erythrocytic, lymphocyte-infecting stages, nor for an absence of transovarial transmission in ticks (biological features distinguishing Babesia spp. from Theileria spp.), and we can only conclude that describing this agent as a Theileria was premature.
The emergence of human babesiosis, among other lines of evidence, has eroded the long-held dogma of piroplasm host specificity, once an important taxonomic character. Our analysis of B. microti suggests that Babesia may be more speciose than previously thought: many researchers have suggested that as modern biochemical or molecular phylogenetic methods were applied to piroplasm systematics, species that were mainly described on the basis of morphology and host associations would be synonymized (Killick-Kendrick, 1974; Levine, 1971, 1984; Telford et al. 1993). The amount of 18S rDNA sequence diversity that we demonstrate in our study is greater than that seen within the Theileria buffeli complex as well as within and between T. parva subspecies (Collins & Allsopp, 1999; Gubbels et al. 2000). Analysis of the ribosomal internal transcribed spacer sequences may allow the measurement of the amount of gene flow between B. microti populations, as it has for T. parva (Collins & Allsopp, 1999) to test the biological distinction of the apparent clades.
Because clades 1 and 3 are morphologically indistinguishable and are found to infect the same species of rodents, there is a degree of uncertainty as to the identity of rodent-infecting piroplasms that previously have been reported. The recent description of Theileria youngi (Kjemtrup, Robinson & Conrad, 2001) from woodrats would also contribute towards the difficulty of specifically identifying rodent piroplasms by morphology. We note that Franca, in his description of ‘Smithia microti’ from voles (1910), also provided descriptions of similar parasites in shrews and other small mammals within the same study sites. Because there is no certain way of assigning ‘S. microti’ to either clade, our finding renders ‘Babesia microti’ incertae sedis. But, because a change in nomenclature would not promote taxonomic stability, we suggest retaining this designation as a means of referring to the monophyletic group (as B. microti sensu lato). For public health purposes, Clade 1 could be referred to as B. microti sensu stricto. Although Clade 3 is genotypically distinct, a consensus has not yet emerged among protozoologists as to molecular criteria for species distinction. The conservative approach is to continue to complement molecular data with careful morphology and life-cycle studies prior to concluding novelty. Nearly 100 species of Babesia have been described; the rules of zoological nomenclature (International Trust for Zoological Nomenclature, 2000) dictate that an entity to be described as a new species must be demonstrated to differ from those that have previously been validly published.
Our work was supported by NIH grants AI 39002, AI 37993, AI 42402, and CDC U50/CCU 119560. H. K. G. was supported, in part, by an NIH Training Grant AI 07535. This is a contribution of the Beringian Coevolution Project (NSF 0196095 and 9981915), and was submitted in partial fulfillment of the requirements for a doctoral degree at the Harvard School of Public Health (H. K. G.). We thank Joseph Cook, Edward Korenberg, Tomas Camacho, the Maine Medical Center Lyme Disease Research Group, and Ivo Foppa for their assistance in providing samples. Stephen Rich and Benjamin Rosenthal provided invaluable advice on phylogenetic analysis.
Table 1. GenBank accession numbers of samples used in the analysis

Fig. 1. Morphology of representative parasites on blood smears. (A) Vole Cordova, AK; (B) white-footed mouse Nantucket Island, MA; (C) Vole Mis, Russia; (D) skunk Martha's Vineyard, MA; (E) dog Galicia, Spain. Scale bars and numbers equal average size of host red blood cells in microns.

Fig. 2. Phylogenetic analyis of Babesia microti using the neighbour-joining algorithm. Shown are bootstrap consensus trees of (A) representative 18S rDNA sequences of Babesia and Theileria (one representative sequence from each B. microti clade is included); (B) 18S rDNA of B. microti clades; (C) beta-tubulin gene of B. microti clades with the 3rd codon included; (D) beta-tubulin gene of B. microti clades with 3rd codon excluded.

Fig. 3. Phylogenetic analysis of Babesia microti clades using maximum parsimony. Shown are bootstrap consensus trees of (A) 18S rDNA gene and (B) beta-tubulin gene with the 3rd codon included.

Fig. 4. Phylogenetic analysis of Babesia microti clades using maximum likelihood. (A) 18S rDNA gene and (B) beta-tubulin gene with the 3rd codon included.
Table 2. Sequence polymorphism within each clade and sequence divergence between clades