Edited by Stanley Falkow, Stanford University, Stanford, CA (received for review March 4, 2004)
Plague, the disease caused by the bacterium Yersinia pestis, has greatly impacted human civilization. Y. pestis is a successful global pathogen, with active foci on all continents except Australia and Antarctica. Because the Y. pestis genome is highly monomorphic, previous attempts to characterize the population genetic structure within a single focus have been largely unsuccessful. Here we report that highly mutable marker loci allow determination of Y. pestis population genetic structure and tracking of transmission patterns at two spatial scales within a single focus. In addition, we found that in vitro mutation rates for these loci are similar to those observed in vivo, which allowed us to develop a mutation-rate-based model to examine transmission mechanisms. Our model suggests there are two primary components of plague ecology: a rapid expansion phase for population growth and dispersal followed by a slower persistence phase. This pattern seems consistent across local, regional, and even global scales.
Plague, caused by the bacterium Yersinia pestis, has impacted humans for thousands of years (1). A relatively young species, Y. pestis most likely originated 1,500-20,000 years ago in Africa (2). Since then, there have been three major plague pandemics. The first pandemic was active in North Africa, Europe, central and southern Asia, and Arabia during the sixth and seventh centuries A.D. During the 14th to 17th centuries A.D., Europe experienced the second pandemic, which included the infamous Black Death (A.D. 1347-1351). The third pandemic is thought to have originated in China in 1855 and continues to this day. In the 1890s, due to international steamship traffic, plague was spread around the world by ship-borne rats. Despite its major impacts on human populations, plague originated as (and continues to be) primarily a disease of rodents and their associated fleas (1). It is established currently in mammal populations on all continents except Australia and Antarctica (3).
Plague arrived in North America at San Francisco around 1900 and moved eastward through native rodent populations (4). Because they are naive and nonadapted hosts, many North American rodents are highly susceptible to plague and experience massive population die-offs (epizootics) after exposure (5). These epizootic hosts are incapable of maintaining plague in the environment. Long-term maintenance of plague is through enzootic or reservoir species, which either are resistant to plague or have high annual recruitment that compensates for plague-caused mortality (6). Although plague maintenance in North America is poorly understood, kangaroo rats (Dipodomys spp.), deer mice (Peromyscus spp.), and grasshopper mice (Onychomys leucogaster) are possible reservoirs (5). Because these species are nocturnal and relatively inconspicuous to humans, plague activity in their populations is rarely observed. As a result, plague activity in North America is usually documented only when it causes epizootics in more conspicuous mammals such as ground squirrels (Spermophilus spp.) and prairie dogs (Cynomys spp.) (3).
Prairie dogs are large (0.5-1 kg), diurnal, burrow-dwelling ground squirrels that form dense colonies in grasslands and are highly susceptible to plague. Mortality rates for Gunnison's prairie dog (Cynomys gunnisoni), the species that occurs in our study area, are typically >99% during plague epizootics (6). Plague occurs throughout the geographic range of this species and is the only disease known to extirpate entire colonies in a short time period (weeks to several months) (7); plague was first documented in Gunnison's prairie dog populations in Arizona in 1938 (8). The specific factors that influence the interspecific transmission of plague from reservoir populations into Gunnison's prairie dog populations are poorly understood but may be triggered by environmental cues such as mild winters and moist springs (9, 10). After the transfer of plague from the reservoir to one or several prairie dogs, it quickly spreads throughout the entire population through a transmission cycle that includes Gunnison's prairie dogs and their two main fleas, Oropsylla hirsuta and Oropsylla tuberculata (3). Thus, plague populations increase rapidly during an epizootic in a prairie dog colony. After the prairie dogs die, plague-infected fleas, which survive in burrows for up to 1 year after epizootics (11), lose their host specificity (3) and are collected easily from the prairie dog burrow systems (11).
Although plague outbreaks have been well described and documented throughout history, there has been little attempt to describe genetic patterns associated with these outbreaks, because the plague genome is extremely homogeneous (2, 12). As a result, genetic tools with the necessary resolution to describe genetic variability, mutation rate, or pattern of spread of Y. pestis within a single outbreak have been unavailable previously. Recently, however, hypervariable regions known as variable-number tandem repeats (VNTRs), which occur throughout the genome of Y. pestis, have been used to distinguish unique genotypes and determine phylogenetic relationships among plague samples collected around the world and within the United States (12, 13). We used a high-resolution multiple-locus VNTR analysis (MLVA) system developed for Y. pestis and global positioning system and geographical information system technologies to address the following questions: (i) What is the regional (total area ≈ 350 km2) pattern of genetic variation in Y. pestis across multiple epizootic events? (ii) What is the local pattern (total area ≈ 100 hectare) of genetic variation in Y. pestis across a single epizootic event? (iii) Can in vitro data such as mutation rates provide insights into in vivo processes such as modeling of plague-transmission cycles?
Regional Study Area and Sampling. In the late spring and summer of 2001, we were actively monitoring northern Arizona prairie dog populations and observed catastrophic but localized die-offs across an ≈350-km2 area near Flagstaff, AZ. (Fig. 1). In the outbreak area, open grassland habitat suitable for prairie dog habitation is limited and has an irregular distribution within the dominant Ponderosa pine vegetation type (white vs. green areas on Fig. 2A ). The boundary between these two habitat types is relatively discrete and often defines the outer boundary of a prairie dog colony, because they will construct burrows directly adjacent to but not in forested areas (Fig. 3B ). Thus, prairie dog colonies in the outbreak area were located in islands of suitable habitat surrounded by unsuitable habitat (Fig. 2 A ), which limited intercolony dispersal. Distances between colonies also limited intercolony dispersal, because prairie dogs seldom disperse >10 km from their natal colony (14). Of the 99 colonies in this area, 49 experienced >99% mortality and were classified as die-offs (Fig. 1). We attempted to collect fleas from ≈50 prairie dog burrows in each of the 49 die-off colonies (for specific methodology, see Supporting Materials and Methods, which is published as supporting information on the PNAS web site).
Distribution and status of Gunnison's prairie dog colonies near Flagstaff. Colored circles indicate the status of colonies as of September 2001. The gray polygon identifies the extent (≈350 km2) of plague activity that affected many of these populations during the spring and summer of 2001. Ninety-nine prairie dog colonies were surveyed within this area. Of the 99 colonies, 49 (49%) experienced die-offs (>99% mortality) between May and September of 2001; Y. pestis was confirmed as the causative agent of 19 die-offs. Colonies were categorized based on the following criteria. Plague (n = 19): (i) evidence of recent (<1 month) prairie dog activity throughout the area (i.e., fresh scat) but no prairie dogs seen or heard or just a few (<1 per hectare) individuals scattered throughout the entire colony area; (ii) Y. pestis-positive fleas collected from at least one prairie dog burrow; and (iii) most burrow entrances open. Die-off (n = 30): same as plague criteria but no Y. pestis-positive fleas collected. Active (n = 19): (i) live prairie dogs observed throughout the area at normal densities; and (ii) most burrows entrances open. Inactive (n = 31): (i) no prairie dogs seen or heard anywhere in the area and only old scat, if any, present; and (ii) most burrow entrances closed.
Spatial (A) and phylogenetic (B) relationships among 39 regional Y. pestis DNA samples from north-central Arizona indicating significant clustering in both geographic and genetic space (for the specific methodology, see Supporting Materials and Methods). Y. pestis DNA was extracted from 39 flea pools collected from 19 different prairie dog colonies. (A) Map of regional study area illustrating the patchy distribution of grassland habitat (white areas) within the dominant forest habitat (green areas) as well as the location of the 19 prairie dog colonies from which samples were collected. Colonies were grouped and assigned a corresponding color code based on phylogenetic analyses of the samples collected from them. Because all samples from any given colony were assigned to the same phylogenetic group, just the locations and group status of the 19 colonies are indicated on the map. (B) Unrooted neighbor-joining tree for the 39 regional samples, which was created in paup (20) by using size data for the 43 Y. pestis MLVA markers and character weighting. Bootstrap values also were generated in paup by using 10,000 simulations. Of the 43 MLVA markers, 22 were polymorphic in the data set and the samples were divided into seven distinct groups, including: Hochderffer (HFR), San Francisco Peaks (SFP), Mormon Lake (ML), Flagstaff (FLG), Government Prairie (GVP), Spring Valley (SPV), and Sitgreaves Mountain (SMT). Samples within the San Francisco Peaks and Sitgreaves Mountain groups separated into two distinct clades based on colony of origin. These within-group clades are distinguished on the map and on the tree by squares and circles of the same color.
Phylogenetic and spatial analyses of the 2001 Ft. Valley plague outbreak. (A) Unrooted phylogenetic analysis based on cladistic principles and maximum-parsimony assumptions (20). We assumed that the dominant genotype (FV-1) was central to the other rarer types and that character-state changes were consistent with in vitro mutation rates and products (Table 1). This unrooted tree contains 24 mutational steps and has a consistency index of 0.80. Individual character-state changes are indicated by lowercase letters. Relative to the FV-1 genotype, these changes are: a, M19Δ-6; b, M27Δ+8; c, M58Δ+17; d, M19Δ-24; e, M34Δ+9; f, M76Δ-41; g, M31Δ+8; h, M79Δ-40; i, M31Δ-8; j, M19Δ-54; k, M34Δ-54; l, M19Δ-6; m, M34Δ-45; n, M34Δ-81; o, M34Δ-90; p, M19Δ+6; q, M34Δ-90; r, M34Δ+9, s, M34Δ+9; t, M19Δ+12; u, M27Δ-8; v, M19Δ-6; w, M75Δ-18; x, M22Δ+7. (B) Individual genotypes are represented by colored symbols and spatially mapped by using arcview. Genotypes observed only once are represented by squares and are numbered. More common genotypes are represented by colored circles and defined in the legend. The FV-1 isolate used to generate the Table 1 data was collected from the black-circled location.
Local Study Area and Sampling. To examine plague dynamics at smaller temporal and spatial scales, we intensely sampled an outbreak in a single prairie dog colony (≈100 hectares). In late May 2001, a prairie dog die-off was reported by residents of the rural Ft. Valley, AZ, community. When we examined the area, we found only remnants of a once-thriving prairie dog colony. Easterly burrows showed no evidence of annual excavation of burrow entrances after winter hibernation, suggesting that prairie dogs in this region died-off over the winter or in the preceding fall. In most of the colony, we observed only scattered individuals, but the westernmost sector still had high densities of prairie dogs in late May, which all died by mid-June. These observations suggest that the outbreak began in the fall of 2000 in the eastern sectors and progressed west after the emergence of prairie dogs from hibernation in the spring of 2001. We attempted to collect fleas from 944 prairie dog burrows in this colony and recorded the coordinates of each burrow with a differential global positioning system receiver, which provided detailed location data for environmental samples.
Analyses. DNA was extracted from all environmental samples (flea pools) and then subjected to Y. pestis-specific PCR to confirm the presence or absence of Y. pestis DNA. Genetic analysis was performed on all positive samples by using 43 VNTR loci (MLVA; Tables 2 and 3, which are published as supporting information on the PNAS web site), which represent some of the fastest-evolving chromosomal sequences known for Y. pestis (13). We used neighbor-joining analysis (15), with a simple matching coefficient, to examine phylogenetic relationships among regional samples and cladistic analysis (maximum parsimony) to examine phylogenetic relationships among samples from the Ft. Valley colony (for the specific methodology, see Supporting Materials and Methods).
Mutational Processes. We used an in vitro Y. pestis population, representing ≈21,000 generations, to define VNTR mutational processes, including mutational products and rates. This process was accomplished by serially transferring each of 96 clonal lineages, which were derived from a single initial colony, for 10 passages and then examining diversity among the 96 terminal (T10) colonies by using MLVA. The in vitro population was derived from the dominant genotype in the Ft. Valley outbreak and isolated from the eastern sector of the outbreak (i.e., the area of the suspected outbreak root; Fig. 3B ). Although in vivo Y. pestis populations may be affected by many factors not present in vitro, this study provided data for constructing evolutionary and epidemiological models to test this assertion.
Molecular Clock. To establish a molecular clock for plague-outbreak dynamics, we developed a plague-transmission model that includes the primary mammalian host and its flea vector, coupled with mutational rates. The model assumes that a single transmission cycle involves Y. pestis subpopulations within a single prairie dog and a single flea, which transmits the infection to the next prairie dog. Using values from the literature (for the specific methodology, see Supporting Materials and Methods), we included in the model best estimates of subpopulation size within the host and vector as well as VNTR mutation rates estimated in vitro (Table 1). The model predicts that 52 Y. pestis population doublings (effective generations) will occur within each transmission cycle and that there is a 6.8% chance of a VNTR mutation occurring and passing onward in any particular transmission cycle. We applied this model to data from the local Ft. Valley outbreak.
Regional Outbreak. Fleas from 19 of the 49 die-off colonies tested positive for plague; plague also was likely the causative agent of the other 30 die-offs (Fig. 1). Regional samples separated into seven distinct groups, which were strongly supported by branch-length and/or bootstrap data (Fig. 2B ). On average, branch lengths separating groups were three times longer than branch lengths within groups, and bootstrap values were >50% for six groups and ≥99% for four groups. Six groups were tightly clustered in geographic space; only the Sitgreaves Mountain group was not spatially grouped, and it had low bootstrap support as well as long within-group branch lengths.
Local Outbreak. Fleas were obtained from 352 of the 944 burrows sampled in the Ft. Valley outbreak; 149 burrows yielded plague-positive flea pools. MLVA analysis indicated two distinct genotypes each from 28 of the 149 burrows, for a total of 177 samples. Twenty-five unique MLVA genotypes were observed (Fig. 3A ). One genotype, FV-1, was dominant (71%) and widely distributed across the study area. It was particularly dominant in the eastern sectors, from which the outbreak is thought to have originated (Fig. 3B ). Fifteen genotypes were observed only once and were widely distributed across the area. Nine minor genotypes were observed multiple times and showed a strong tendency for spatial clustering. The topology, based on parsimony analysis, approaches that of a star phylogeny, with most genotypes only a single mutational step from the dominant FV-1 genotype (Fig. 3A ). Rooting such a phylogenetic topology is problematic, but the overall local dominance of FV-1, its presence in regionally adjacent locations (Fig. 2, Flagstaff group), and its dominance at the earliest local outbreak sites suggest it is the evolutionary root for this Y. pestis population. Only seven two-step branches are derived from the central FV-1 genotype, and these are spatially consistent with their proposed phylogenetic relationships. For example, the single FV-20 genotype is located immediately adjacent to several individual FV-4 genotypes, and the two FV-8 genotypes are located immediately adjacent to FV-2 genotypes (Fig. 3B ).
In Vitro Y. pestis Population. Twenty-six mutational events were observed across 10 different loci within this population (Table 1). Overall, single-repeat mutations were most common, and insertions were only slightly more common than deletions. The fastest mutating loci in vitro were also among those that varied the most in the Ft. Valley outbreak (Fig. 3A ) and regionally (data not shown). Slightly different mutation rates for loci with the same number of mutations are caused by missing data. Combined, the mutation rate for all 43 loci examined was 1.3 × 10-3 mutations per generation.
Transmission Model. Given the most parsimonious topology for Ft. Valley genetic patterns (Fig. 3A ), 24 mutations occurred across the 177 observations in this outbreak. Our transmission model predicted 355 transmissions cycles for the 177 data points, or 2.0 transmission cycles per sample point. From empirical field data and burrow density at the Ft. Valley site (for the specific methodology, see Supporting Materials and Methods), we estimated that a total of 472 individual prairie dogs lived in this region, all of which died in the plague outbreak. If the entire population was included within our sampling, it would have represented 2.7 transmission cycles per sample point and 32 mutational events would have been expected. It is reasonable to believe the maximum-parsimony assumption is slightly underestimating mutational events and, also, not all prairie dog-flea transmission events are represented among our 177 sample points. If so, the already close predicted and expected values would be even closer.
By using MLVA, we defined Y. pestis population structure on a landscape scale, which provides unique insights into plague ecology. First, plague is probably regionally widespread in mammalian reservoir populations. This pattern is difficult to examine, however, because the reservoirs are cryptic. Second, the lack of robust deeper phylogenetic topology (star phylogeny) and the clustering of similar genotypes in geographic space suggest the initial spread of plague across the landscape was probably quite rapid and possibly involved a single, common genotype, which later differentiated within local, dispersal-limited reservoir populations. Third, because plague-infected prairie dog colonies were interspersed with healthy colonies (Fig. 1), interspecific transmission events from the reservoir to prairie dogs probably occurred not once but multiple times. This pattern suggests these events are independent but coordinately triggered by environmental cues (9, 10). Fourth, in areas in which prairie dog colonies are located close to each other (<10 km) and separated by suitable grassland habitat (e.g., Fig. 2 A , within groups Hochderffer, Government Prairie, Mormon Lake, and Spring Valley), intercolony dispersal of plague is probably through infected prairie dogs. Indeed, samples from these four discrete grassland areas form distinct genetic groupings that are strongly supported by bootstrap data (Fig. 2B ). Fifth, the presence of unsuitable, forested habitat between colonies located close to each other (<10 km) likely prevents intercolony dispersal of plague through prairie dogs (e.g., Fig. 2 A , between groups Government Prairie and Spring Valley). Sixth, the presence of similar genotypes in prairie dog colonies that are separated by long distances (>10 km) and unsuitable habitat (e.g., Fig. 2 A , within groups Sitgreaves Mountain, Flagstaff, and San Francisco Peaks) is probably because of intercolony dispersal of plague-infected fleas through domestic dogs, coyotes, or other predator/scavengers (16, 17).
The fragmented population structure of the outbreak host in this area, along with interspersed reservoir species habitat, allows a dynamic balance between the amplification of the pathogen and its long-term persistence, which seems ideal for a plague focus. Amplification of the pathogen can occur in a limited fashion that does not drive the susceptible host to extinction. At the same time, these amplification events also increase the likelihood of interspecific transmission back into the reservoir population, permitting survival of the pathogen between outbreaks in the susceptible host.
Despite the rapid progression of localized outbreaks, our findings demonstrate that Y. pestis population sizes during a single outbreak are large enough that mutational events are observed in multiple VNTR loci. These mutational events can provide insights into the dispersal of plague during an outbreak. The genotype diversity patterns in the outbreak described here (Fig. 3) are consistent with rapid wave-like dispersal by a single fit genotype (FV-1) followed by continued, localized Y. pestis population growth. The fitness of FV-1, which could be caused by either stochastic or adaptive processes (18), led to a dominant overall genotype but with diverse, locally derived genotype clusters. Overall, the phylogenetic relationships are spatially consistent. Either dispersal is slow enough that evolutionarily derived relationships are spatially conserved or dispersal is rapid but spares many of the prairie dogs during the initial plague progression. The differentiated population then could be associated with a second wave of prairie dog deaths; differentiated populations also may be associated with growth inside surviving, infected fleas.
Whereas the occurrence of VNTR mutations can provide insights into plague dispersal during an outbreak, understanding the rate at which these mutation events occur provides additional insights into transmission processes. The strong association between observed and expected VNTR mutation events indicates our model is approximating the natural transmission processes. Overall, the flea-prairie dog transmission model reasonably explains the Ft. Valley plague outbreak and illustrates that the generation of VNTR variation in vivo is not so different from that observed in vitro. However, if only the frequency and distribution of genotype FV-1 are considered, the model seems inadequate. The model would require FV-1 to be widely dispersed in only a few weeks and no mutations to occur during the dispersal. The probability of 124 successive transmissions without a mutation would be 1.7 × 10-4 in our model. Thus, we believe the initial dispersal of the FV-1 contagion was not strictly a burrow-to-burrow progression but rather involved either a smaller number of mammalian hosts distributing plague-infected fleas widely or another completely different transmission mechanism. A small number of initially infective animals may have spread infected fleas across the entire outbreak area, causing rapid flea-borne transmission. Alternatively, plague can spread without vectors through respiratory droplets or feeding on infected carcasses. These more rapid transmission mechanisms would require smaller Y. pestis population sizes, reducing the likelihood of mutations. Regardless of the actual mechanism, our data would be consistent with an initial rapid dispersal by a few animals, followed by a slower, localized, vector-mediated process.
At both the landscape and local scales, the genetic data suggest very rapid plague transmission followed by slower localized differentiation, which is evidenced by long branch lengths accompanied by a lack of structure in the deeper branches of our phylogenetic hypotheses (Figs. 2 and 3). A “star phylogeny” is consistent with a common ancestor from which multiple lineages are derived simultaneously (19). In our data, these multiple lineages are spatially dispersed across our study areas at both the regional and local scales. At the local scale, the likely ancestral genotype is observable (FV-1), and its rapid dispersal occurred only weeks before complete extinction of the Ft. Valley prairie dog colony. In contrast, the regional common ancestor can only be postulated, but its wide dispersal must have been in the last century, because plague was first documented in Arizona in 1938 (8) and in North America in 1899 (4). It is interesting to note that this pattern of rapid and wide dispersal followed by localized differentiation is also observed at the global scale. During each of the three major plague pandemics that occurred over the last 1,500 years, there was rapid dispersal followed by regional establishment of plague foci that have existed, in some cases, for millennia (2).
We postulate that when plague encounters a susceptible rodent species that is both plague-naive and highly dense (perhaps social), there is rapid dispersal across the landscape. A natural consequence is a rapid decline in susceptible host density, which would extinguish the epidemic. In some cases, the rapid dispersal is followed by slower transmission cycles in mostly resistant or low-density hosts as the disease is established into a stable reservoir for future emergence. Whereas the pathogen would go extinct without the reservoir, it would also be very rare and limited in distribution without the explosive expansion in susceptible hosts. This dynamic balance explains the emergence and subsequent success of plague at the global and regional scales, and even in localized outbreaks.
We thank B. Davis, J. Gallie, D. Krpata, C. Levy, P. Service, S. Shuster, T. Theimer, and W. Van Pelt. Funding was provided by the U.S. Department of Energy, Arizona Game and Fish Department, National Fish and Wildlife Foundation, and Kaibab National Forest.
↵ § To whom correspondence should be addressed. E-mail: paul.keim{at}nau.edu.
↵ † J.M.G. and D.M.W. contributed equally to this work.
This paper was submitted directly (Track II) to the PNAS office.
Abbreviations: VNTR, variable-number tandem repeat; MLVA, multiple-locus VNTR analysis.
