|Home | About | Journals | Submit | Contact Us | Français|
Understanding the spread of infectious diseases is crucial for implementing effective control measures. For this, it is important to obtain information on the contemporary population structure of a disease agent and to infer the evolutionary processes that may have shaped it. Here, we investigate on a continental scale the population structure of Borrelia burgdorferi, the causative agent of Lyme borreliosis (LB), a tick-borne disease, in North America. We test the hypothesis that the observed population structure is congruent with recent population expansions and that these were preceded by bottlenecks mostly likely caused by the near extirpation in the 1900s of hosts required for sustaining tick populations. Multilocus sequence typing and complementary population analytical tools were used to evaluate B. burgdorferi samples collected in the Northeastern, Upper Midwestern, and Far-Western United States and Canada. The spatial distribution of sequence types (STs) and inferred population boundaries suggest that the current populations are geographically separated. One major population boundary separated western B. burgdorferi populations transmitted by Ixodes pacificus in California from Eastern populations transmitted by I. scapularis; the other divided Midwestern and Northeastern populations. However, populations from all three regions were genetically closely related. Together, our findings suggest that although the contemporary populations of North American B. burgdorferi now comprise three geographically separated subpopulations with no or limited gene flow among them, they arose from a common ancestral population. A comparative analysis of the B. burgdorferi outer surface protein C (ospC) gene revealed novel linkages and provides additional insights into the genetic characteristics of strains.
Populations of infectious agents are shaped by evolutionary and demographic processes, as well as by the population dynamics of their hosts and (in the case of vector-borne pathogens) arthropod vectors. These processes can leave signatures in the pathogens' genomes. Their frequency and distribution in space and time can be exploited via strain typing for important utilitarian purposes, such as associating specific genotypes with specific ecological niches and determining the geographic distribution of genotypes or of phenotypic characteristics, such as pathogenicity. By investigating the evolutionary history of pathogens, we may be able to elucidate how environmental drivers have in the past shaped patterns of their spread or dispersal and infer how they may do so in the future. Here, by studying Borrelia burgdorferi, the causative agent of Lyme borreliosis (LB) in North America, we investigate the population structure of a vector-borne agent on a continental scale and test the hypothesis that it reflects recent environmentally driven demographic events (population bottlenecks) experienced by vector populations and their reproductive hosts.
LB is a tick-borne bacterial infection which, if not diagnosed and treated early, may develop into a debilitating multisystemic disorder (48). The LB group of spirochetes currently comprises 18 named species, but B. burgdorferi sensu stricto (hereinafter called B. burgdorferi) is the species regularly associated with human disease in North America (34). In North America, B. burgdorferi is transmitted primarily by the blacklegged tick, Ixodes scapularis, east of the Rocky Mountains and by the western blacklegged tick, Ixodes pacificus, in the Far-Western United States, and a large variety of vertebrates serve as reservoir hosts (e.g., see reference 50). In the United States, LB risk currently occurs primarily in three regions: the Northeast, where >80% of reported cases in the United States occur, the Upper Midwest, where 11% of cases occur, and the Pacific coastal region (California and Washington) (2). In Canada, LB is also an emerging infectious disease that became a national reportable disease in 2010. Risk of exposure to vector ticks has been identified in the Maritimes and southeastern Quebec, southern central Canada (Ontario and Manitoba), and southern British Columbia (38–41).
Geographic variation in LB incidence arises in part from environmental factors (climate, natural habitats, and associated anthropogenic changes) that affect the occurrence, abundance, and activity of tick vectors and hosts, the efficiency of transmission cycles, regional variation in the frequencies of more- or less-pathogenic genotypes, and/or the behavior of humans (40, 50). Thus, the spatial distribution of human cases is driven by a combination of factors that vary regionally (17, 21, 22, 31) and that may have been influenced by recent historical events impacting the distribution of ticks and their reproductive hosts (reviewed in references 3 and 46).
The historical distribution of B. burgdorferi in North America can perhaps be inferred from early records of deer and vector distributions (3, 18, 36). It has been suggested that white-tailed deer had been driven almost to extinction by 1900 through both unmanaged exploitation and habitat loss, which would also have affected the availability of other hosts for ticks (7, 12, 24, 36). Deer are key hosts for adult I. scapularis, and their decline would be expected to have driven I. scapularis populations to a low level, which in turn probably produced severe bottlenecks in B. burgdorferi populations (reviewed in references 8, 28, 43, and 46). Reforestation caused reexpansion of deer populations (as well as that of other woodland hosts), followed by expansions of both vector-tick and B. burgdorferi populations (3, 46). Several recent studies reporting on contemporary I. scapularis and B. burgdorferi populations found that their current distribution in the Northeast and Midwest is discontinuous (4, 13, 14, 26). We hypothesize that these recent and past demographic events have left genetic signatures in the genome of B. burgdorferi, which we sought to investigate in the present study.
Multilocus sequence typing (MLST) allows the characterization of diverse bacterial populations (1, 32, 47) by analyzing the genetic variation of multiple loci encoding proteins essential for cell maintenance, or “housekeeping.” Previous studies using MLST suggested that North American B. burgdorferi populations were subdivided between the Northeast and Upper Midwest (26, 41). Other authors, using different genetic markers, such as the plasmid-encoded outer surface protein C (ospC) or the 16S-23S rRNA intergenic spacer region (IGS), suggested that the population genetic structure of B. burgdorferi isolates from the Northeast and Upper Midwest overlaps that from California (4, 5, 22). It has, however, been proposed that the genetic variation in ospC does not reflect the organisms' evolutionary history but, as a result of recombination and/or horizontal plasmid transfer, it instead reflects the evolution of the locus (4).
Here, we seek to elucidate the evolutionary history of B. burgdorferi in North America and its impact on the contemporary population structure using MLST (26, 33, 38, 41). Expecting that multiple chromosomal genetic markers will better reflect the organism's history, we determine whether the genetic structure of B. burgdorferi in North America fits a pattern consistent with long-term and/or more recent environmental changes. For this, we analyzed strains from the Northeast (New England and southeastern Canada), the Upper Midwest (including south-central Canada), and the Far West (California) in North America to investigate potential population subdivisions in an explicitly spatial context. Parallel to the MLST analyses, we typed the ospC locus to understand how variation at this locus may relate to the population divisions as determined by MLST.
Overall, of the 295 samples analyzed in this study, 162 from Canada were described previously by Ogden and coauthors (41), 78 from the Northeast and Upper Midwest by Hoen and coauthors (26), and 7 from California by Margos and coauthors (33). Two GenBank submissions, i.e., strains WI91-23 (genome project ID 28627) and CA11.2A (genome project ID 28629) (45), were included in MLST only, not the spatial analysis. An additional 25 samples from the Upper Midwest and 21 from California were analyzed by MLST for the purpose of this study. Although the number of samples from the Western coastal region is relatively small, the sample set includes representatives of the most common strains found in this region. A summary table of samples (Table 1) is given above, and a complete list of the samples evaluated and their geographic origins is given in Table S1 in the supplemental material.
I. scapularis ticks (mainly adults) from Canada were collected from 2005 to 2007 by means of a passive surveillance program from companion animals and humans at veterinary clinics or medical clinics as described previously (41). The collections comprised ticks from resident populations and others that probably had dispersed on migratory birds from locations in the United States where the I. scapularis and B. burgdorferi populations are established (10, 42). Questing adult I. scapularis ticks were sampled in Upper Midwestern sites between 2004 and 2008 by cloth dragging (25). Tick collections from California consisted of questing I. pacificus nymphs that were collected in 2004 from 78 dense woodlands in Mendocino County (16, 22). We are aware that adult ticks may not pose a risk for spirochete transmission, but the pathogens they harbor in many of our samples are still representative of the local Borrelia populations and constitute an important component of the overall population genetic pool of B. burgdorferi in North America. Although some ticks collected in Canada and their associated Borrelia organisms may not be from locally established transmission cycles, they potentially represent the propagule pool of Borrelia that may eventually become established locally. These samples were therefore included in analyses that explored bacterial relationships. However, Borrelia samples collected in Canada from regions that are not known to have resident tick populations were excluded from the spatial population analysis (n = 191; see Table S1 in the supplemental material). Samples were stored either in 70% or 95% ethanol until processed for DNA purification using Qiagen DNeasy blood and tissue purification kits (Qiagen) (16, 25, 41).
Nested PCRs for the eight housekeeping genes (clpA, clpX, nifS, pepX, pyrG, recG, rplB, and uvrA) were carried out using HotStarTaq (Qiagen, Germany) and a touchdown PCR protocol for the primary reaction (35, 38). For the secondary reaction, an annealing temperature of 50°C was chosen. Denaturation and elongation were run at temperatures of 95°C and 72°C for 30 s and 1 min, respectively.
PCR fragments were sequenced in forward and reverse directions (Qiagen) and manually compared using DNASTAR (Lasergene). Sequences that contained two base peaks at the same position in forward and reverse sequences were considered mixed infections and, therefore, were not included in our analyses. Besides the 293 samples we processed, data for two samples were obtained from GenBank (45).
Sequences for the individual genes were compared to sequences previously deposited in the borrelia.mlst.net database, and unique sequences were given new, consecutive allele numbers. Allele numbers for all eight loci make up the allelic profile that determines the sequence type (ST) for each strain. Novel allelic profiles were given consecutive ST numbers.
The plasmid-encoded ospC gene was amplified in a majority of samples using primers and conditions following Bunikis and coauthors (6).
To investigate the genetic population structure of B. burgdorferi, we used six analytical tools available for MLST data. The first three assess contemporary population structure, whereas the last three assess the population structure with respect to the relatedness of strains; both permit inferences of past events. (i) The geographical distribution of STs was visualized using ArcGIS version 10 and the website http://www.spatialepidemiology.net; (ii) a permutation test was performed (Margos et al. ) to assess to what extent the spatial distribution of STs corresponds to geographic location; (iii) the geographic occurrence of genetic boundaries was identified by “wombling” (11) to establish whether the genetic data fit the notion of genetically separated populations in North America. The relatedness of strains from the populations in different geographic locations was explored using (iv) eBurst and goeBurst (19, 20), (v) Bayesian analysis of population structure (BAPS) (9), and (vi) phylogenetic analysis (23). These analyses provided information on several hierarchical levels, as follows: (i) the spatial distribution of STs and population boundaries provide information about the contemporary population structure; (ii) the information contained in the allelic profiles corresponds to the degree of relatedness between STs, and the eBURST and goeBurst analyses can be used to infer the pattern of descent; and (iii) the BAPS analysis and gene genealogies use the information available in the sequences and provide information about the historical past of the populations. Detailed information on the software and settings used can be found in the supplemental material.
All housekeeping gene sequences have been deposited to the B. burgdorferi MLST database hosted at Imperial College London, United Kingdom, and are available at http://borrelia.mlst.net. ospC sequences have been submitted to GenBank under accession numbers JQ951094 to JQ951225.
The 295 B. burgdorferi samples analyzed were resolved into 85 sequence types (STs) (Table 1; also see Table S1 in the supplemental material). In general, the distribution of STs corresponded to geographic location. None of the STs determined for the 28 samples from California were present in the samples from the other geographic regions (Upper Midwest or Northeast) (Fig. 1; also see Fig. S1 in the supplemental material). Only 5 of the 65 STs found in the Northeast and Upper Midwest occurred in both regions (STs 12, 29, 19, 221, and 222). This differential spatial distribution of ST frequencies in the subpopulations was corroborated by a permutation test using the allelic profiles of B. burgdorferi (see Table S1). This test provided significant evidence (dissimilarity value, 5.5; P = 0.0001) (also see Fig. S2 in the supplemental material) that the geographic distribution of B. burgdorferi STs found in the Northeast, the Upper Midwest, and California differed. The test also provided significant support for spatial differences when Canadian samples from regions without established I. scapularis populations (n = 191; dissimilarity value, 5.5; P = 0.001; data not shown) were omitted.
Wombling of 191 samples from regions with known resident I. scapularis populations revealed two population boundaries that divided North American B. burgdorferi populations into three subpopulations. A pronounced boundary was located between longitudes 93°W and 110°W, which coincided with the Great Plains and the Rocky Mountains and separated the Californian and Upper Midwestern populations (Fig. 2). The boundary was also identified with more stringent and more relaxed settings, resulting in slightly differing geographic ranges (data not shown). A more accurate identification of the division of B. burgdorferi subpopulations between the Upper Midwest and California will require the investigation of more samples from both geographic regions. A second potential boundary was located between longitude 83°W (western Ohio) and 88°W (western Indiana), thereby separating the Northeastern and Upper Midwestern populations. Figure 2 shows the population boundaries determined using wombsoft (h-value of 1.4, and pB = 0.4). The boundary between the Northeastern and Upper Midwestern sites became more pronounced when only samples from these two regions were considered (data not shown).
Despite the contemporary geographic structure of North American B. burgdorferi populations, the eBURST/goeBURST, BAPS, and phylogenetic analyses revealed that these populations of B. burgdorferi are genetically closely related (Fig. 3 and and4).4). Application of the eBURST algorithm divided the strains into 17 clonal complexes. Founder strains were identified with reasonable bootstrap support (approximately 60%) in four clonal complexes, whereas founders were less certain in the other complexes (Fig. 3 and Table 2). In the goeBURST diagram, the groupings of strains were more “bushy” than in previous analyses and, consistent with previous data, the clonal complexes consisted of strains from the different geographic regions, especially the Northeast and Upper Midwest (Fig. 3) (26). STs from California were mostly found to be triple-locus variants (TLV) of strains from the Northeast and Upper Midwest, except for ST403 and ST2. ST403 was the first single-locus variant (SLV) for ST1, an abundant ST in the Northeast, and it connects ST1 with ST2. These findings further support the genetic relatedness of the subpopulations. Similarly, the Californian ST13 was found to be an SLV of ST12, which occurs in the Northeast and Upper Midwest.
To corroborate the findings made by eBURST, we used BAPS because it infers population structuring. The BAPS analysis divided the strains into 10 groups (Table 3 and Fig. 3), which mirrored the relationship assigned by the eBURST analysis. In Fig. 3, the BAPS groups are projected onto the goeBURST diagram. Small differences were discerned by the two analyses. Specifically, ST55 (a TLV of ST30) was assigned to BAPS group 7, whereas ST30 and ST314 were assigned to group 5. ST396 was assigned to group 2, while in goeBURST, it formed a double-locus variant (DLV) with ST3, a member of BAPS group 7. These differences, however, were also reflected in the phylogeny, wherein samples assigned to groups 5 and 7 were split into different phylogenetic clusters (Fig. 4).
STs that demonstrated significant admixture between groups in the BAPS analysis are given in Table 4. These STs shared genetic information from several groups, suggesting they had undergone recombination. This included STs that manifested disagreement between goeBURST and BAPS: ST30 and ST314 both yielded significant admixture.
The phylogenetic analyses revealed that strains from all three regions were represented in the major lineages (Fig. 4). Although some of the strains from California were closely related to strains from other regions, no single ST from California was found in either of the other regions. In fact, two STs from California, 398 and 399, were distant from neighboring clades. These two STs had unique alleles for most of the genes, except rplB (see Table S1 in the supplemental material).
The major ospC groups that had been determined for many of the STs are shown adjacent to the tree in Fig. 4. Many of the major ospC groups matched the ST clusters (Fig. 4; also see Table S1). Some STs, however, were matched with more than one major ospC group, for example, ST12 was found to carry group M ospC and group O ospC. Other major ospC groups, such as group A, were matched locally with different STs, e.g., ST1 in the Northeast, ST2 and ST5 in California, or ST55 in the Upper Midwest. Several major ospC groups fell into different phylogenetic clusters (e.g., N or H), presumably as a result of horizontal transfer.
In this study, we analyzed North American B. burgdorferi populations by MLST in an explicitly spatial context. Samples from the geographic regions most affected by LB, that is, the Northeast and Upper Midwest, and from a region in northern coastal California with a high prevalence of tick infection (29), were included in our analyses. All samples were derived from ixodid ticks (nymphs and adults), collected either by dragging vegetation (Northeast, Upper Midwest, and California) or by means of passive surveillance (southern Canada) between 2004 and 2008.
The geographic distribution of STs suggests that present day populations are separated along geographical lines: the Northeast, Upper Midwest, and California. The permutation test provided evidence that the distribution of STs was significantly different among the three regions, and we hypothesize that this probably reflects the local expansion of regional refuge populations in the Upper Midwest and Northeast during a time much more recent than events reflected in the phylogenetic tree.
Furthermore, the inference of genetic boundaries by wombling analysis between North American B. burgdorferi populations provided additional support for the present-day population structure. A western population boundary separated B. burgdorferi populations east and west of the Great Plains and the Rocky Mountains, landscape features that are credible natural barriers to B. burgdorferi gene flow. The finding of closely related strains, i.e., SLV east and west of the Rocky Mountains, suggests some population connectivity despite these landscape features. Although mutations in housekeeping genes accumulate slowly and the genes evolve over long periods, the uplifting of the Rocky Mountains dates back to about 50 to 100 million years ago (http://www.geo.arizona.edu/geo5xx/geo527/Rockies/uplift.html), which is likely to be further in the past than the mutation(s) found in an SLV. It is possible that the following hypotheses may alone or collectively explain gene flow among these geographic regions: (i) occasional Borrelia transport occurred (and may be still occurring) across these barriers (for example, by birds); (ii) the now divided subpopulations of B. burgdorferi constituted an admixed population when deer (36) and, possibly, ticks and B. burgdorferi had a much wider distribution range than at present; and (iii) gene flow occurred in an ancestral admixed population that might have been located further south during Pleistocene glacial periods from 2.5 million years ago to about 18,000 years ago (see below) (27).
The location of the boundary between the Northeastern and Upper Midwestern B. burgdorferi subpopulations is roughly consistent with a subdivision described for B. burgdorferi ospC at a longitude of 83°W (4). This boundary coincides with a region where a previous study described an absence of I. scapularis (13). However, a model for risk of Lyme borreliosis in North America based on the density of host-seeking nymphs predicted higher nymphal densities in northern Indiana and Ohio than were observed during the project (14), which suggests that the division of B. burgdorferi populations found between the Northeast and Upper Midwest may be the result of a population contraction as a consequence of vector population decline. The collapses of white-tailed deer, vector tick, and presumably, Borrelia populations that occurred from 500 to 110 years ago (46) are consistent with our findings and those of others (6, 26, 44) that local B. burgdorferi populations expanded out of refuge populations in the Northeast and Upper Midwest. If so, this scenario may explain the clonality and linkage disequilibrium between plasmid-encoded and chromosomal genes described in the initial population studies on B. burgdorferi in the Northeast (6, 44). Our study confirms and extends results from previous studies showing that on a wider geographic scale, the initially reported linkage disequilibrium between ospC and chromosomal markers was not as pronounced (5, 49). However, B. burgdorferi populations, especially at the edges of their distributional ranges, are likely to show an element of clonality that is related to their parasitic life style, i.e., the bacterium resides only in the tick or its vertebrate host. For recombination or horizontal gene transfer to occur, close physical contact among bacteria is required and depends upon the frequency of infection in vertebrate hosts or tick vectors with multiple strains.
Although the geographic distribution of STs provided information about the contemporary population structure of B. burgdorferi in North America, BAPS, eBURST, goeBURST, and phylogeny yielded information about its putative ancestral population structure. We hypothesize that the signals observed using these analytical tools reflect an ancestral population structure that existed before the recent population bottleneck caused by the decline of deer and the consequent effects on tick vector and Borrelia populations (3, 13, 14, 18, 36). Our findings suggest that before these events, B. burgdorferi strains formed overlapping populations and/or that gene flow between populations was more pronounced than today. For example, the denser “forest” shown in the goeBURST diagram and the additional clonal complexes obtained compared to the ones described before (26) underpin the genetic relatedness of strains from the different regions. In addition, the BAPS analysis largely agreed with the relationships inferred from goeBURST and the phylogenetic tree. Together, these data suggest that to understand the evolutionary history of B. burgdorferi as a species, samples from all three regions must be considered.
Bootstrap support for some clonal complex founder assignments was slightly higher than 60%. This was probably due to a lack of SLV, which may have been caused by loss of genetic diversity as a result of the severe bottleneck. However, SLV of STs from the Upper Midwestern and Northeastern sites were found in California (e.g., ST12/ST13 and ST403/ST1), suggesting that additional strains like this may be found during more exhaustive analyses of samples from these regions. An increased number of SLV or DLV will improve the bootstrap support for clonal complex founders and will increase our understanding of the direction of gene flow. Previous studies suggested that gene flow occurred in an east-to-west direction in North America (5, 26). This conclusion is consistent with our current and previous eBURST results that most founders of clonal complexes are found in the Northeast (26), with the fact that B. burgdorferi is absent in Asia and therefore would not have populated North America via the Bering Strait, and with the hypothesis that B. burgdorferi originated in Europe (33). SLV and DLV were found more frequently in Northeastern and Upper Midwestern sites than in either the Northeast and California or the Upper Midwest and California, intimating a closer relationship between subpopulations of Northeastern and Upper Midwestern spirochetes. Moreover, the analysis using BAPS found significant admixtures between strains from the Northeast and Upper Midwest, providing support for an overlap of these B. burgdorferi subpopulations in the past. We speculate that the current population structure reflects the strong bottleneck that B. burgdorferi has undergone due to loss of forest habitats (at least in the Northeast and Upper Midwest) and that this bottleneck caused a severe loss of diversity in the Northeast.
The finding of identical STs in the Upper Midwest and Northeast suggests that either these populations once overlapped, reflecting a random distribution of these STs in both regions at the onset of the bottleneck, or alternatively, there is currently a limited gene flow between the two regions. Differentiating between these possibilities requires the use of more sensitive methods, such as the analysis of genome-wide single nucleotide polymorphisms or ecological studies to determine the host associations of these STs. In this respect, it is interesting to note that although migratory birds are considered important to transport ticks and, possibly, tick-borne pathogens over long distances (39), their movements, which respond (directly or indirectly) to seasonal changes in temperature-driven resource availability, are mostly north-south in direction rather than east-west, and this behavior probably evolved or emerged due to increasing availability of northern habitat during interglacials (37). An alternative potential explanation for east-to-west spread of LB spirochetes might be spread via rodent hosts or their ancestors. This would require the concomitant spread of ticks and/or the presence of established vector tick populations. It is conceivable that Ixodes ticks may have had a much wider distribution range when its main reproductive host, the deer, was more widespread (36). For some rodent hosts, phylogeographic patterns similar to that seen here for B. burgdorferi have been shown, i.e., Northeastern, Midwestern, and Western populations, with the Midwest and Californian clades being linked geographically (and via interbreeding at hybridization zones) by Rocky Mountain-Great Plains clades (15, 51). Such a pattern is consistent with the phylogeographies of Hanta virus and its deer mouse reservoirs (15). We recognize, however, that gene flow might have occurred in ancestral populations located further south in North America.
In agreement with previous studies (4, 49), the phylogenetic analysis revealed that strains that carried the same ospC genotype were matched regionally with different MLST sequence types. Although the samples from California represented the most common major ospC groups present in questing I. pacificus nymphs in this region, ospC type H3 was described as the most prevalent in questing I. pacificus nymphs in northwestern California (22). Curiously, the three strains carrying ospC type H3 (ST398 and ST399) were genetically the most divergent in our analysis, suggesting that ST398, ST399, and perhaps, related strains probably represent a long-term separated population in California. Ecologically, ospC type H3 strains have been associated with hardwood-dominated forests where the western gray squirrel (Sciurus griseus) is the primary reservoir of B. burgdorferi (30). Epidemiologically, these strains have not been associated with human patients, and it had been speculated that their high prevalence might be one reason for the low human LB incidence in northern California (22).
In conclusion, this study has revealed that B. burgdorferi in North America comprises three separate but genetically related populations. These are separated by geographic barriers, with California populations being the most isolated, while Northeastern and Upper Midwestern populations overlap and could have been isolated more recently by land use changes. Nevertheless, both the geographic occurrence and the evidence for close relatedness of the three populations also correspond with known phylogeographic patterns among rodent hosts for the vector tick and bacterium, perhaps influencing current B. burgdorferi phylogeographic structure. Further studies of B. burgdorferi phylogeography, coupled with similar studies on vertebrate hosts and tick vectors, are required to enhance our understanding of how environmental drivers in the past have shaped the population and geographic structure of B. burgdorferi. Our study provides optimism that mining signatures in the B. burgdorferi genome may assist in prediction of the effects of current and future environmental changes on the risk of vector-borne diseases.
We are grateful to M. Hurn for support with statistical analyses, S. Manel for help interpreting genetic boundaries, D. Aanensen for maintaining the MLST database, and Natalia Fedorova for critical comments on the manuscript. We are indebted to Klaus Kurtenbach (died March 2009) for his inspirational ideas.
Published ahead of print 22 June 2012
Supplemental material for this article may be found at http://aem.asm.org/.