Local factors associated with on‐host flea distributions on prairie dog colonies

Abstract Outbreaks of plague, a flea‐vectored bacterial disease, occur periodically in prairie dog populations in the western United States. In order to understand the conditions that are conducive to plague outbreaks and potentially predict spatial and temporal variations in risk, it is important to understand the factors associated with flea abundance and distribution that may lead to plague outbreaks. We collected and identified 20,041 fleas from 6,542 individual prairie dogs of four different species over a 4‐year period along a latitudinal gradient from Texas to Montana. We assessed local climate and other factors associated with flea prevalence and abundance, as well as the incidence of plague outbreaks. Oropsylla hirsuta, a prairie dog specialist flea, and Pulex simulans, a generalist flea species, were the most common fleas found on our pairs. High elevation pairs in Wyoming and Utah had distinct flea communities compared with the rest of the study pairs. The incidence of prairie dogs with Yersinia pestis detections in fleas was low (n = 64 prairie dogs with positive fleas out of 5,024 samples from 4,218 individual prairie dogs). The results of our regression models indicate that many factors are associated with the presence of fleas. In general, flea abundance (number of fleas on hosts) is higher during plague outbreaks, lower when prairie dogs are more abundant, and reaches peak levels when climate and weather variables are at intermediate levels. Changing climate conditions will likely affect aspects of both flea and host communities, including population densities and species composition, which may lead to changes in plague dynamics. Our results support the hypothesis that local conditions, including host, vector, and environmental factors, influence the likelihood of plague outbreaks, and that predicting changes to plague dynamics under climate change scenarios will have to consider both host and vector responses to local factors.


| INTRODUC TI ON
Sylvatic plague, caused by the bacteria Yersinia pestis, manifests in the western United States as episodic outbreaks in prairie dogs, squirrels, and small rodents, with occasional spillover to human populations (Gage, 1998). Prairie dogs (Cynomys spp.) are highly susceptible to plague and colonies often suffer >90% mortality during plague outbreaks (Cully & Williams, 2001). By decimating prairie dog populations, plague can have devastating effects on species such as black-footed ferrets (Mustela nigripes), which depend on prairie dogs for their prey (Antolin et al., 2002). At present, the factors driving plague outbreaks are difficult to ascertain. A number of studies have indicated that small rodent composition (Stapp et al., 2009) and local weather conditions (Collinge et al., 2005;Hubbart, Jachowski, & Eads, 2011;Savage, Reich, Hartley, Stapp, & Antolin, 2011) may play a role in contributing to plague epizootics, potentially through their effects on flea distribution and abundance. However, these factors likely vary geographically in the strength of their effects. As the primary vector of plague, fleas play an important role in the dynamics of the disease (Bacot & Martin, 1914). In order to understand the conditions that are conducive to plague outbreaks and potentially predict spatial and temporal variations in risk, it is important to understand the factors associated with flea abundance, prevalence, and distribution that may lead to plague outbreaks.
Traditionally, the hypothesized mechanism for flea transmission of plague bacteria was the result of "blocking" (Bacot & Martin, 1914), whereby Y. pestis bacteria form a blockage in the midgut of fleas, causing them to increase the number of feeding attempts and regurgitate infectious material, increasing the likelihood of plague transmission. Lorange, Race, Sebbane, and Hinnebusch (2005) developed models that suggested, with this mechanism of transmission, high flea abundance was necessary for driving plague epizootics due to the low competence of fleas as vectors and the short lifespan of blocked fleas. Eisen et al. (2006) introduced the idea of early-phase transmission to explain the rapid spread of plague through host populations, as an alternative hypothesis to blocking.
Evidence of early-phase transmission, which does not require blocking of fleas, was observed in Oropsylla montana (Eisen et al., 2006), and modeling has demonstrated that early-phase transmission can drive plague epizootics (Buhnerkempe et al., 2011). Previous research has demonstrated that the ability of fleas to block decreases at higher temperatures (Cavanaugh, 1971;Hinnebusch, Fischer, & Schwan, 1998;Kartman, 1969), therefore potentially explaining the relationship between decreased temperature and frequency of plague outbreaks (Cavanaugh, 1971;Collinge et al., 2005;Enscore et al., 2002;Pollitzer, 1954). However, in contrast, early-phase transmission in Xenopsylla cheopis was delayed or inhibited at cooler temperatures, but did not seem to be affected by warmer temperatures (Schotthoefer et al., 2011). These results indicate that local climate conditions may affect flea biology and the likelihood of plague outbreaks in prairie ecosystems.
Numerous factors likely affect on-host flea prevalence (number of hosts with fleas) and abundance (number of fleas on hosts).
Previous research has indicated that on-host flea abundance has been associated with local weather conditions, although the relationships are seemingly complex. Increasing flea abundance has been associated with dry summers in a New Mexico population of black-tailed prairie dogs (BTPD, C. ludovicianus;Eads, Biggins, Long, Gage, & Antolin, 2016), wetter springs and cooler summers on great gerbils (Rhombomys opimus) in Kazakhstan (Stenseth et al., 2006), and warmer autumns in the previous year in Mongolian gerbils (Meriones unguiculatus; Xu et al., 2015), while decreasing flea abundance was associated with increasing prior year growing season and winter precipitation in South Dakota BTPD (Eads & Hoogland, 2016).
Several studies have examined the relationship between the incidence of plague and climatic variables. Some studies have proposed that hot weather desiccates flea larvae and reduces overall flea abundance (Snäll, O'Hara, Ray, & Collinge, 2008), therefore, reducing the likelihood of plague, while others have proposed that drier weather leaves hosts in poorer condition resulting in increased flea abundance and increased potential for plague outbreaks . Precipitation also may affect the likelihood of plague, with increased precipitation leading to better forage conditions and increased host densities, which may create favorable conditions for plague outbreaks (Enscore et al., 2002). For example, Stenseth et al. (2006) found that the incidence of Y. pestis in great gerbils in Central Asia increased with warmer springs and wetter summers. Eads and Hoogland (2017) observed that flea abundance on Gunnison's prairie dogs (C. gunnisoni, GPD) in Arizona and white-tailed prairie dogs (C. leucurus, WTPD) in Colorado varied inversely with increasing precipitation during the prior year growing season and suggested that this may explain increased risk of plague after drought. Enscore et al. (2002) found that the numbers of human plague cases in New Mexico and Arizona decreased with increasing summer temperatures and increased when spring precipitation in previous years was high. However, greater summer rainfall, but not previous year precipitation, was associated with plague events in BTPD in Colorado (Savage et al., 2011). At last, Eisen et al. (2012) suggested that flea diversity, which was positively related to rainfall and negatively related to temperature, may play a role in plague outbreaks in Uganda.
In addition to climate and weather variables, flea and host characteristics may influence flea abundance. Oropsylla hirsuta and O. tuberculata, flea species strongly associated with prairie dog colonies and plague outbreaks, exhibit seasonal trends (Mize & Britten, 2016).
The abundance of O. tuberculata is often highest in early spring (Cully, Barnes, Quan, & Maupin, 1997;Salkeld & Stapp, 2008;Tripp, Gage, Montenieri, & Antolin, 2009), while O. hirsuta numbers tend to peak in mid or late summer (Salkeld & Stapp, 2008;Tripp et al., 2009). Eads et al. (2013) noted that overall prevalence of fleas also fluctuated seasonally. Flea abundance on BTPD has been found to be higher on adult males than on juvenile and/or adult female prairie dogs in Colorado . However, in South Dakota, no difference in flea abundance was found on male and female BTPD, although juveniles had fewer fleas than adults (Eads & Hoogland, 2016).
In summary, the relationships between weather factors, flea abundance, and host characteristics (body condition or abundance) remain difficult to ascertain and may not be consistent between different host species, fleas, and locations. We collected fleas from four species of prairie dogs over a 4-year period (2013-2016) on paired plots located along a latitudinal gradient from Montana to Texas and assessed local climate and other factors associated with flea prevalence and abundance (total fleas, O. hirsuta and P. simulans counts), as well as the incidence of plague outbreaks. Elucidating the relationship between flea species composition, abundance and distribution, the occurrence of plague, and environmental factors will provide information to management agencies responsible for controlling plague, as well as help predict how plague dynamics may change under future climate conditions.

| Study areas and design
We analyzed flea abundance and prevalence on prairie dogs from 23 paired plots that had been included in a large-scale sylvatic plague vaccine (SPV) trial (see further details in Rocke et al., 2017). A pair consisted of one plot treated with SPV-laden baits and one plot that received placebo baits. In brief, our study included 11 paired plots on BTPD colonies, one pair on a GPD colony, four pairs on WTPD, and seven pairs on Utah prairie dog (UPD, Cynomys parvidens) colonies, sampled over a 3-year period, 2013-2015 ( Figure 1). Fifteen of the pairs were resampled in 2016. For plots to be included in the study, we required that no pesticides had been used for at least one year prior to the beginning of the study. Most pairs had never been dusted; however, two pairs (4 plots) (Charles M Russell, Montana-CMR) were dusted in 1997, and one pair (Wind Cave, South Dakota-WCSD) was dusted in 2011. Tripp et al. (2016) observed a waning effect of dusting on flea abundance on sites in Colorado after 10-12 months.
SPV-laden baits were distributed on one randomly assigned plot in each pair and placebo baits were distributed on the other. Within 1-4 weeks, prairie dogs were trapped and sampled on the same days at both plots within a pair by the same personnel using Tomahawk live traps. Traps were baited with sweet feed and oats. Sex, age, weight, and foot length were recorded for each animal sampled, and fleas and other samples were collected from each animal upon first capture during a sampling period. Sampling dates varied by year, but occurred between June and November (Appendix 1), and a sampling period consisted of a minimum of 3 days of trapping (on consecutive days unless prevented by field conditions) (Appendix 2). During 2013, prairie dogs were sampled both pre-and post-baiting, with approximately a month between sampling periods. In 2014-2015, each plot was sampled once (with the exception of the GPD pair and one of the WTPD pairs in Utah which were sampled twice, once in July and once in August).
For flea sampling, animals were anesthetized with isoflurane in an induction chamber and immediately combed for fleas that were collected into tubes of sterile saline or ethanol. After sample collection, animals were allowed to fully recover from anesthesia and then released at the location of capture (NWHC Animal Care and Use Committee, Protocol #EP130214). Fleas collected from live prairie dogs were stored at −20°C until identification. To remove surface contamination that might alter PCR results, fleas were rinsed with 70% ethanol containing 0.2% iodine and rinsed twice in sterile water (La Scola, Fournier, Brouqui, & Raoult, 2001;Raoult et al., 2001;Kernif et al., 2014). Fleas were then counted and identified to species using published taxonomic references (Furman & Catts, 1982;Hubbard, 1947;Stark, 1970), and pooled by species and sex, up to 10 individuals per pool from a single animal.
To determine plague status of our study plots, spleen and liver tissues from prairie dog carcasses, as well as flea pools collected from live and dead animals, were tested for the presence of Y. pestis using real-time PCR (see Rocke et al., 2017 for details). In brief, tissue DNA was extracted using the Wizard SV Genomic DNA Purification System (Promega; Madison, WI), and PCR was performed for the caf1 gene sequence located on the Y. pestis pG8786 plasmid (Genbank accession NC_006323) and the pla gene located on the Y. pestis pPCP1 plasmid (Genbank accession AL109969). Tissues from carcasses were also cultured on blood agar plates, and suspect colonies were DNA samples were then screened for the pim gene, which resides on the pPCP1 plasmid of Y. pestis, using real-time PCR, as we found this gene to be more sensitive than the pla gene for fleas. Suspect positive fleas were then confirmed using caf1 gene as described above.
A plot was classified as plague positive if at least one carcass or one flea pool from a live prairie dog tested positive for Y. pestis by culture or PCR. Once a plot was classified as "plague positive," it was considered plague positive in subsequent years.

| Statistical analyses
We do not formally account for detection error in our analyses; therefore, we have no repeated measures from which to estimate detection (see Eads et al., 2013 for methods with repeat detections).
The amount of time spent combing for fleas was not standardized between pairs (although we assume effort between paired plots was similar); therefore, our counts should be considered relative indices.
We used logistic regression and negative binomial regression with LASSO (least absolute shrinkage and selection operator; Tibshirani, 1996) priors in a Bayesian framework to assess factors related to total flea abundance, O. hirsuta abundance, and Pulex simulans abundance. This allowed us to examine factors related to both the number of fleas on a host (abundance) and the presence/absence of fleas on a host (prevalence). The LASSO allows for model selection and regularization to take place automatically, while improving predictive accuracy and providing interpretable parameter estimates (Tibshirani, 1996). Parameter estimates for variables that are not important predictors shrink to zero. We demonstrate fit of the data to a negative binomial distribution graphically (Appendix 3).
We evaluated factors associated with prevalence of fleas and overall flea abundance, including variables related to climate and seasonal weather patterns (NDVI [normalized difference vegetation index-a measure of vegetation greenness, see Table 1], number of days with a temperature over 85°F [Collinge et al., 2005], amount of precipitation during the prior year's growing season, winter precipitation amount, and day of sampling), plot-level variables (catch per unit effort as an index of relative abundance of prairie dogs [see Rocke et al., 2017  for BTPD (n = 22 plots), UPD (n = 14), and WTPD (n = 8). Three of seven UPD pairs were located on the same prairie dog colony and significant movement was noted between the plots on these pairs, so treatment effects at these pairs may be compromised. GPD were collected at one pair of plots only, therefore, accounting for location was not necessary. For O. hirsuta and P. simulans abundance, we used the same covariates as described above; however, negative binomial regression was only appropriate for BTPD, due to the relatively few numbers of hosts with >1 flea observed for other species. Therefore, for all other species, we report results for logistic regression only.
We estimated factors associated with initial plague detection on placebo plots using logistic regression, including all plot-level factors previously mentioned and flea species diversity as predictor variables. Data were summarized for each plot by year combination, and because we were interested in determining factors related to initial plague outbreaks (the first year plague was detected in fleas or carcasses), plot by year combinations with continuing plague outbreaks were excluded. We were interested in evaluating conditions that could predict plague outbreaks; therefore, we included catch per For all models, logistic and negative binomial, we calculated effect sizes by exponentiating the parameter estimates. Effects estimated from parameters of a negative binomial indicate the change in the expected number of fleas on a host given a 1 unit increase in the covariate (or a change in the covariate from one categorical level to another); we represent these effects as a percent change. Effects for coefficients from a logistic regression are odds ratios and indicate the odds of flea presence given a 1 unit increase in the covariate.
Odds ratios of 1 indicate no change in the probability of presence given an associated increase in the covariate value, odds ratios<1 indicate that the probability of presence is less likely given an associated increase in the covariate value, and odds ratios >1 indicate an increase in the probability of presence.
All models were run in R (R Core Team, 2017). Models were run for a total of 60,000 iterations with a burn in of 20,000. All models were run in rjags (Plummer, 2013). We used Gelman-Rubin diagnostics to assess convergence (Gelman & Rubin, 1992). We used area under the curve statistics in the "pROC" package (Xavier et al., 2011) to evaluate goodness of fit for our logistic regression models.
For negative binomial models, we used posterior predictive checks on the parameters of the negative binomial distribution by simulating 1,000 data sets from the estimated model parameters and comparing to the observed data (Gelman, Carlin, Stern, & Rubin, 2004). To display a visual representation of flea community structure, we used nonmetric multidimensional scaling (NMDS; an ordination technique) in R (package "vegan"; Oksanen et al., 2017). We assessed the fit of the NMDS by reporting a "stress" value (Oksanen et al., 2017).

| Climate and weather variables
Multiple weather and climate-related variables were associated with flea abundance and presence/absence. Total flea abundance increased 92% (95% Credible Interval [C.I.] 9-217%) on UPD pairs over every 10-day period (Figure 2). The odds of a prairie dog having F I G U R E 2 Effect sizes (represented as average percent change in flea abundance) as a function of (a) climate and environmental variables (precip=precipitation, see Table 1 (Table 3), while P. simulans abundance decreased 66% (95% C.I. 50-78%) ( Table 3). On UPD pairs, day of sampling was associated with an increase in the odds of O. hirsuta presence (1.01, 95% C.I. 1.12-1.23) but not P. simulans.
With a 10% increase in NDVI, total flea abundance decreased by 89% (95% C.I. 80-95%) on BTPD pairs, but these differences were not significant for individual flea species (Figures 2 and 3). ing winter precipitation. These results were largely driven by declines F I G U R E 3 Odds ratios from results of logistic regression relating climate and weather variables (precip=precipitation, see Table 1 for description of parameters) to (a) flea presence (all flea species), (b) Oropsylla hirsuta presence, and c) Pulex simulans presence for black-tailed prairie dog (Cynomys ludovicianus; black dots), Gunnison's prairie dog (Cynomys gunnisoni; gray dots), Utah prairie dog (Cynomys parvidens; black squares), white-tailed prairie dog (Cynomys leucurus; gray squares). *Indicates odds ratios that do not overlap one  (Table 3).

| Plot-level characteristics
Plot-level variables also indicated associations with flea abundance and presence/absence. The presence of plague was positively associated with flea abundance for all species except GPD (Figure 2). In addition, O. hirsuta and P. simulans abundance was associated with plague detection for BTPD (Table 3) for BTPD plots, while average flea abundance declined 97% (95% C.I.

F I G U R E 4
Odds ratios from results of logistic regression relating plot-level variables (see Table 1 for description of parameters) to (a) flea presence (all flea species), (b) Oropsylla hirsuta presence, and (c) Pulex simulans presence for black-tailed prairie dog (Cynomys ludovicianus; black dots), Gunnison's prairie dog (Cynomys gunnisoni; gray dots), Utah prairie dog (Cynomys parvidens; black squares), white-tailed prairie dog (Cynomys leucurus; gray squares). * Indicates odds ratios that do not overlap one

| Individual characteristics of prairie dogs
Flea abundance increased by 6% (95% C.I. 3-9%) for every 1 unit increase in body condition (g/mm) for WTPD, but no other species

| Flea species composition
Results of our NMDS reflect results of our regression analyses ( Figure 6). High elevation pairs in Wyoming and Utah (Appendix 7) had distinct flea communities compared with the rest of the study F I G U R E 5 Odds ratios from results of logistic regression relating individual characteristics of prairie dogs (see Table 1 for description of parameters) to (a) total flea abundance, (b) Oropsylla hirsuta abundance, and (c) Pulex simulans abundance for black-tailed prairie dog (Cynomys ludovicianus; black dots), Gunnison's prairie dog (Cynomys gunnisoni; gray dots), Utah prairie dog (Cynomys parvidens; black squares), whitetailed prairie dog (Cynomys leucurus; gray squares). * Indicates odds ratios that do not overlap one

| Plague detections
The    (Lorange et al., 2005). Abundant fleas are also likely to result in more efficient early-phase transmission (Eisen et al., 2006), whereas flea presence merely indicates that at least one flea is present on the host. Therefore, a large number of hosts with low flea abundance (i.e., high flea prevalence) may not result in plague outbreaks. We did not find any relationship between local weather conditions and the likelihood of plague outbreaks.
However, our sample sizes were small and geographically varied. It is likely that conditions under which plague epizootics occur are as locally-specific as the conditions that result in high flea abundance.
Winter precipitation increases resulted in decreases in flea abundance for all species except WTPD, while overall flea abundance decreased with increasing number of days with temperatures above 85°F on WTPD plots but increased with warmer temperatures on BTPD plots. These results may reflect differences in species composition or local weather conditions. Salkeld and Stapp (2008) found increasing O. hirsuta and decreasing O. tuberculata populations in prairie dog burrows with increasing temperature, which may explain some of these differences. Soil type is an additional factor that may influence flea abundance. Savage et al. (2011) found that soils with increased soil moisture holding capacity were related to the incidence of plague epizootics; however, Salkeld and Stapp (2008) found no relationship between O. hirsuta abundance and soil type. Despite the overall increase in flea abundance on WTPD pairs, P. simulans abundance decreased with increasing winter precipitation on these plots (as well as BTPD plots), indicating that the increases in flea abundance on WTPD plots with winter precipitation are driven by other flea species. In addition, it is likely that flea abundance is optimal in a mid-range of conditions. On WTPD plots, average winter precipitation was roughly 2 cm greater than on BTPD plots, indicating that greater precipitation on these wet plots led to decreases in flea abundance in contrast to plots that received lesser amounts of precipitation.
Several possible mechanisms may explain the effects of local weather on flea abundance and plague outbreaks. Local weather factors may affect flea abundance directly by influencing flea mortality or reproductive rates (Krasnov, Khokholova, Fielden, & Burdelova, 2001) or indirectly by changing forage availability and affecting host body condition and/or grooming behavior . Flea abundance on BTPD plots declined with increasing NDVI, supporting the hypothesis that better forage conditions on BTPD plots may lead to lower flea abundance. In our data set, body condition was not strongly associated with flea abundance of most species, except WTPD, and for this species, animals with better body condition tended to have more fleas, possibly due to their larger size. It is possible that body condition may increase in importance during drought conditions which were not observed on our plots during 2013-2016.
One study noted that BTPD that survived plague outbreaks tended to have better body condition than prior to the die-off, potentially due to less competition for forage (Pauli, Buskirk, Williams, & Edwards, 2006). This observation may explain the relationship we observed between flea abundance on plague-positive plots and body condition in our study. We did not have direct before-after comparisons; however, there was no difference in average body condition between animals captured on plots with plague detections (mean body condition = 14.0, SD = 4.4), compared to plots with no plague detection (mean body condition = 14.5, SD = 4.6).
Previous research has indicated that trends in flea abundance during a season are highly dependent on the flea community  may be the result of fleas congregating on surviving animals . It is also possible that the combination of high flea abundance and high prairie dog abundance leads directly to plague epizootics and is therefore rarely observed (i.e., once these conditions are met there is a plague outbreak which drives prairie dog densities lower). Lastly, the combination of high densities of animals in summers that are drier than average could lead to poorer condition of animals overall, leading to higher susceptibility to plague the following summer, particularly if forage availability (as indicated by NDVI) is low.
Changing climate conditions will likely affect aspects of both flea and host communities, including population densities and species composition, and lead to changes in plague dynamics. We expected that plague outbreaks would be less likely when temperatures were high (Snäll et al., 2008); however, we did not observe this relationship in our data set. High temperatures have been observed to interfere with transmission of Y. pestis in rat fleas (Cavanaugh, 1971) and may lead to decreases in plague outbreaks with warming temperatures. In an alternative manner, hosts and plague vectors may expand their range northward, thus shifting the range of Y. pestis outbreaks (Nakazawa et al., 2007). A meta-analysis recently pro-

CO N FLI C T O F I NTE R E S T
None declared.

AUTH O R CO NTR I B UTI O N S
T.R., D.T., and R.R. conceived the ideas and design of the study; R.R. and T.R. designed the methodology; R.A. and D.T. collected the data; R.R. analyzed the data and R.R., T.R., D.T., and R.A. contributed to analysis interpretations; R.R. led the writing of the manuscript. All authors contributed critically to the drafts and gave final approval for publication.