Host biology and environmental variables differentially predict flea abundances for two rodent hosts in a plague-relevant system.

While rodents frequently host ectoparasites that can vector zoonotic diseases, often little is known about their ectoparasite communities, even in places where hosts frequently interact with humans. Yosemite National Park is an area of high human-wildlife interaction and high potential zoonotic disease transfer. Nonetheless, relatively few studies have surveyed the flea communities on mammalian hosts in this area, and even fewer have characterized the environmental and host factors that predict infestation. We focused on two species, the alpine chipmunk (Tamias alpinus) and the lodgepole chipmunk (T. speciosus), which inhabit Yosemite and surrounding areas and can host fleas that vector plague. Because these hosts are exhibiting differential responses to environmental change, it is valuable to establish baselines for their flea communities before further changes occur. We surveyed fleas on these chipmunk hosts during three years (2013-2015), including in the year of a plague epizootic (2015), and documented significant inter-host differences in flea communities and changes across years. Flea abundance was associated with host traits including sex and fecal glucocorticoid metabolite levels. The average number of fleas per individual and the proportion of individuals carrying fleas increased across years for T. speciosus but not for T. alpinus. To better understand these patterns, we constructed models to identify environmental predictors of flea abundance for the two most common flea species, Ceratophyllus ciliatus mononis and Eumolpianus eumolpi. Results showed host-dependent differences in environmental predictors of flea abundance for E. eumolpi and C. ciliatus mononis, with notable ties to ambient temperature variation and elevation. These results provide insight into factors affecting flea abundance on two chipmunk species, which may be linked to changing climate and possible future plague epizootics.

Introduction

Although small mammals can serve as reservoirs of ectoparasite-vectored pathogens (Gage et al., 1995; Mills and Childs, 1998), in many cases relatively little is known about the ectoparasite communities living on these hosts. This is true even for host species that frequently interact with humans and therefore pose a threat of zoonotic disease transfer. The Sierra Nevada Mountain Range of California is both a hub of small mammal diversity (Grinnell and Storer, 1924) and an area of heavy human-wildlife interaction, with Yosemite National Park alone visited by over 3.7 million people annually (https://www.nps.gov/yose/planyourvisit/visitation.htm). As a result, it is also an area of high potential zoonotic disease transfer (Adjemian et al., 2008). Most recently, an epizootic of plague, a flea-vectored bacterial disease caused by Yersinia pestis, was responsible for two human infections in 2015, in addition to the closure of numerous campgrounds and many documented rodent fatalities, largely California ground squirrels (Otospermophilus beecheyi) and golden mantled ground squirrels (Callospermophilus lateralis; Danforth et al., 2016). Despite these risk factors, relatively few studies have characterized the flea communities hosted by small mammal species in this locality (Jameson and Brennan, 1957; Stark, 1970; Adjemian et al., 2008; Smith et al., 2010; Fleer et al., 2011; Danforth et al., 2016). Such work is especially important in the context of climate change, which has increased average temperatures in the park by approximately 3 °C and has transformed vertebrate communities throughout the Sierra Nevada range (Moritz et al., 2008; Tingley et al., 2012).

Like most organisms, ectoparasites could be affected by anthropogenic environmental change (Colwell et al., 2012; Carlson et al., 2017). To better predict how these communities might change, it is critical to understand how they are influenced by both host biology and external environmental factors. Certain host traits are often predictive of flea abundance; for example, male mammalian hosts often carry higher numbers of ectoparasites (Zuk and McKean, 1996; Schalk and Forbes, 1997), and individuals in poor body condition often carry more fleas, though in some cases the opposite is true (Hawlena et al., 2005; Krasnov et al., 2005). Flea densities also tend to be positively correlated with host population densities (Arneberg et al., 1998; Krasnov et al., 2002; Young et al., 2015). While fleas exhibit variable and species-specific responses to temperature and humidity (Lang et al., 1996; Krasnov et al., 2001; Kreppel et al., 2016), it has been hypothesized that increased variability in climate regimes, which is projected in several climate change scenarios (Raisanen, 2002; Della-Marta et al., 2007; Fischer and Schär, 2009; Hansen et al., 2012; Vázquez et al., 2017), may generally benefit pathogens and parasites above hosts, possibly in part because their smaller body sizes and higher metabolic rates increase their rate of acclimatization (Rohr et al., 2010, 2013; Raffel et al., 2013). Because high flea densities are associated with increased rates of plague transmission (Krasnov et al., 2006; Pham et al., 2009; Tripp et al., 2009) and the geographic range of plague is thought to be mediated by vector ecology (Maher et al., 2010), understanding the ecological factors that impact flea abundance is increasingly important as climate change alters the environment.

Here we focus on two rodent species that inhabit the Sierra Nevada mountains, the alpine chipmunk (Tamias alpinus) and the lodgepole chipmunk (T. speciosus), and the fleas living on them. In the published literature, we could find only one account characterizing the flea species that parasitize T. alpinus, consisting of one flea found on a single individual (Fleer et al., 2011). While this high-elevation Sierra Nevada endemic is less frequently in contact with humans, its flea community deserves attention given the lack of past documentation and the changing spatial distribution of this host (Moritz et al., 2008). The elevational range of T. alpinus in Yosemite National Park has contracted significantly upwards over the last century, presumably due to warming temperatures (Mortiz et al., 2008). Tamias speciosus is known to carry flea species that vector plague and has frequently tested seropositive for Y. pestis (Eskey and Haas, 1939; Barnes and Kartman, 1960; Nelson, 1980; Gage and Kosoy, 2005; Adjemian et al., 2008; Holt et al., 2009; Smith et al., 2010; Fleer et al., 2011; Danforth et al., 2016; Dubyanskiy and Yeszhanov, 2016; Straub et al., 2017). This species co-occurs with T. alpinus but lives across a much broader range of elevations (Moritz et al., 2008) and can frequently be found in areas of high human use. We compare and contrast the flea communities carried by these hosts, examine interannual differences in flea abundances, and test which phenotypic (sex, mass, fecal glucocorticoid metabolites) and environmental traits (population densities, elevation, date, and ambient temperature regimes) are predictive of flea abundance for each of the two most dominant flea species on each host. In doing so, we provide a thorough account of the flea communities found on both host species and begin to tease apart the traits of hosts and environments that may contribute to site- and individual-level variations in ectoparasites of small mammals.

Methods

Study sites and species

The alpine chipmunk (Tamias alpinus, Ta) is a small (30–50 g) species found mainly at and above the treeline (∼3000 masl; Clawson et al., 1994). The lodgepole chipmunk (T. speciosus, Ts) is a larger (50–80 g) species that is found mainly at and below the treeline (∼1900–3200 masl; Best et al., 1994; Moritz et al., 2008). Between June and October of 2013–2015 we sampled populations at seven to ten sites per year (Fig. 1; Supplementary Data S1). Grids (∼2–5/site) of ∼40–80 Sherman live-traps were set out at each site and baited with peanut butter and oats. Traps were arranged in pairs with ∼10 m between pairs. Traps were opened at dawn, closed at dusk, and checked approximately every 4-6h during the day. Each captured chipmunk was identified to species, sexed, weighed, measured, and ear-tagged with uniquely numbered metal tags (1005–1, National Band and Tag Co. Newport, KY). Animals that were trapped in multiple years and had lost their ear-tags between years, preventing individual re-identification, were excluded from statistical analyses. Within each calendar year, only data from the first capture of each individual were included in analyses.

Fig. 1

Map showing study sites. Sites (see Supplementary Data S1 for more information) located in and around Yosemite National Park (green) were visited either in all three years (2013, 2014, and 2015; black), in two of the years (yellow), or in only one year (red). Yellow and black lines show significant roadways in the area. Lakes are shown in blue, including Mono Lake at top right. Inset shows Yosemite National Park (green) on a map of California. Site codes: AL: Arrowhead Lake; CL: Cathedral Lake (upper); GA: Glen Aulin; GL: Gaylor Lakes; HC: Hoffmann Creek; MA: Mammoth Lakes; ML: May Lake; PC: Porcupine Creek; SL: Saddlebag Lake; SLN: Saddlebag Lake, north-side (Greenstone and Steelhead Lakes); TM: Tuolumne Meadows. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Fecal glucocorticoid metabolite analyses

Glucocorticoids are important metabolic hormones that are often related to an organism's allostatic load (Sapolsky et al., 2000). Because parasites can impact host allostatic loads, we integrated fecal glucocorticoid metabolite levels, which typically reflect an average of circulating glucocorticoids over many hours (Touma and Palme, 2005; Palme, 2019), in our analyses of flea abundances for each host species. Feces were collected from each captured individual as part of another ongoing study and the methods for their collection and analysis are described in detail in Hammond et al. (2015, 2018). Briefly, fecal pellets were collected directly from each individual's capture trap and immediately frozen. Samples were dried to remove water, crushed, weighed, extracted with methanol, dried, and dissolved in the assay buffer before being assayed with a 5α-pregnane-3β,11β,21-triol-20-one enzyme immunoassay first described for use with laboratory mice (Touma and palme, 2005), which measures FGMs with a 5α-3β,11β-diol structure. This assay has been previously validated for both chipmunk species (Hammond et al., 2015).

Flea collection and identification

Each chipmunk was scruffed and combed with a metal-toothed flea comb five times down the dorsum, each back leg, and the tail. Chipmunks were not combed above a white plastic tray and were not anesthetized for this procedure. Flea combing was performed by multiple researchers. Fleas appeared in between the teeth of the comb and were collected with forceps or fingers by a second researcher. In some cases, fleas disappeared back into the pelage before they could be captured; in these cases, if another flea was later observed, to be conservative we assumed it was the same flea and did not double count that flea. Some fleas (∼20% of counted fleas) visibly jumped away from the collection area and were counted but not collected. All collected fleas were stored in 100% ethanol in cryogenic vials (1 vial per collection time per host), which were frozen within one month upon return to the University of California, Berkeley campus.

Collected fleas were identified to species as in Pigage et al. (2017). 122 flea specimens were cleared, dehydrated and mounted in Canada balsam on microscope slides using standard techniques (Lewis et al., 1988). They were then identified to species using several keys (Hubbard, 1947; Traub et al., 1983; Lewis et al., 1988; Lewis and Jameson, 2002). The majority of fleas were maintained in 100% ethanol, observed microscopically and identified using the mounted fleas as references. All fleas on microscope slides were added to the permanent collection of fleas at the Denver Museum of Nature and Science (accession numbers ZP.2000–2176). Unmounted fleas were turned over to the California Department of Public Health for potential disease testing. All procedures were approved by the University of California, Berkeley IACUC (protocol #2015-01-7084) and followed the American Society of Mammalogists guidelines for the use of wild mammals in research (Sikes et al., 2016).

Environmental data collection

Ambient temperatures at each trapping grid were recorded using iButton thermochron temperature loggers (model DS1921G). Loggers were deployed within 1 m of the substrate at approximately 75% of the trap stations on each grid. Loggers were programmed to collect hourly readings. From these data, average maximum, mean, and variation in temperature were calculated for each trapping grid on the day of capture for each individual. Due to correlations between these summary statistics and in order to reduce dimensionality, principal components analysis (PCA) was applied to these data to generate a smaller number of uncorrelated temperature variables that could be included in models (Table 1). The first principal component (which aligned positively with maximum, mean, and variance in temperature) and second principal component (which aligned most strongly, positively with variance in temperature) were used as fixed effects in models (see below).

Table 1

Principal components analysis results for ambient temperature data collected with iButtons.

Coefficient	PC1	PC2	PC3
Maximum Temperature	0.66	−0.15	−0.74
Mean Temperature	0.61	−0.46	0.64
Variance in Temperature	0.44	0.87	0.22
Eigenvalues	1.96	0.81	0.23
S.D.	1.40	0.90	0.48
Proportion of Variance	0.65	0.27	0.08

Mark-recapture data were used to generate estimates of chipmunk density (including both species, when appropriate) at each trapping grid. Relative population density was quantified using the Schnabel index, a method that is appropriate when populations have been sampled more than twice (Napolitano et al., 2008; Schnabel et al., 1938). This index is defined by the equation:where M is the total number of individuals captured at time t, C is the number of marked individuals in the population just before t, and R is the number of previously marked individuals captured at t. This value was calculated for each trapping grid for each year and then divided by that grid's area to generate an estimate of relative population density. Grid area was calculated by establishing a ∼20 m buffer around each trapping station, merging buffers for all trapping stations in the same grid, and calculating the resulting area. Density values were included as fixed effects in statistical models. In areas where both chipmunk species co-occur, we used combined density values that took both species into account.

Statistics

Interspecific comparisons and host traits— All statistics were implemented in R (R Core Team, 2017). Flea count data was non-normal based on Shapiro tests and visual inspection of the data, thus, non-parametric tests were used. Non-parametric Wilcoxon rank sum and Spearman's rank correlation rho tests were used to examine relationships between flea abundance and species, year, sex, mass, and glucocorticoids. Flea abundance was defined as the number of fleas per individual, whether or not the individual was infested or un-infested, and to calculate average abundances both infested and un-infested hosts were included (Bush et al., 1997). P-values were false-discovery rate adjusted to account for multiple testing (Benjamini and Hochberg, 1995).

Models of Flea Abundances and Environmental Factors—Poisson or negative binomial models are most commonly used for non-normal count data. Likelihood ratio tests and Dean's overdispersion tests were implemented in the Dcluster package; results indicated that a negative binomial model better fit our over-dispersed data in comparison to a Poisson model. Negative binomial mixed models, implemented in the glmmTMB package (Magnusson et al., 2017), were used to examine relationships between flea abundances and environmental traits while controlling for site, year, and certain host biological parameters. All continuous fixed effects were mean-rescaled prior to analysis to allow for comparison of effect sizes across variables. Pairwise plots and VIF values of fixed effects were assessed prior to model construction to test for multicollinearity, which was not found. Graphs (QQ plots, plots of residuals vs. predicted values and vs. fixed effects) and statistical tests (over-dispersion, uniformity, outliers, zero inflation), implemented in the R package DHARMa, were used to assess residuals to validate all models and ensure that they did not violate assumptions (Hartig, 2017; Zuur et al., 2009). A small number of outliers were removed from each model in order to meet assumptions (8 or fewer points with response variable values of over four standard deviations from the mean).

Separate models were constructed for each of the two most common flea species, Ceratophyllus ciliatus mononis and Eumolpianus eumolpi, and for overall flea counts. Models could not be fitted for other individual flea species due to low sample size. Fixed effects in models included interactions between host species and each of the following: elevation, and the first two principal components of temperature data (Table 1) collected on the focal individual's trapping grid on the day of sampling. Julian date was included as a continuous fixed effect. Based on the aforementioned analyses of relationships between host traits and flea abundances, an interaction term between sex and glucocorticoid levels was also included in the models. The significance of each fixed effect term was assessed using the Satterthwaite approximation in the ‘lmerTest’ package (Kuznetsova et al., 2017). Site and year were included as random effects in all models; significance of the random effect was determined using a likelihood ratio test (‘anova’ function).

Results

Flea communities and host differences

Over 1100 fleas from five genera and at least seven species, largely from the Family Ceratophyllidae, were collected from both chipmunk hosts at multiple sites (Table 2, Table 3). Ceratophyllus ciliatus mononis, Eumolpianus eumolpi, and E. eutamiadis were the three most commonly found fleas for both chipmunks, with the relative frequencies of these species differing between hosts: Eumolpianus eumolpi was the most common flea collected from T. alpinus while C. ciliatus mononis was the most common flea collected from T. speciosus (Table 2, Fig. 2C and D). Aetheca wagneri, an unknown Catallagia sp., Oropsylla idahoensis, and Oropsylla montana were also infrequently collected from both chipmunk species (Table 2). On average, the larger-bodied T. speciosus (N = 920 individuals combed) hosted higher numbers of fleas per individual than T. alpinus (N = 260 individuals combed; Wilcoxon rank sum test, W = 131960 p = 0.03; Fig. 2A–B). For T. speciosus but not T. alpinus, the average number of fleas per individual increased significantly from year to year (Spearman's rank correlation rho, T. alpinus: S = 3024200, r = −0.03. p = 0.73; T. speciosus: S = 112690000, r = 0.15. p = 0.0004; Fig. 2A–D). While a slightly different set of sites was visited in each year, this pattern held even if analyses were restricted to sites visited in all years, or to sites within specific ranges of elevations.

Table 2

Fleas collected from The “# Hosts” column lists the numbers of infested (I) and uninfested (U) individual chipmunk hosts from which fleas were collected in each year. Subsequent columns present the total counts of fleas of each species, followed in parenthesis by the prevalence of each flea species (proportion of total chipmunks tested that hosted at least one individual). Prevalence across all years is shown in the last two rows. Proportions do not sum to one as many animals had no fleas, while others carried multiple flea species. Ccm: Ceratophyllus ciliatus mononis; Eem: Eumolpianus eumolpi; Eet: Eumolpianus eutamiadis; Aw: Aetheca wagneri; Cat.: Catallagia sp.; Oi: Oropsylla idahoensis; Om: Oropsylla montana.

	Year	# Hosts	Ccm	Eem	Eet	Aw	Cat.	Oi	Om
T. alpinus	2013	I: 32U: 31	3 (0.05)	63 (0.46)	4 (0.05)	1 (0.02)	0 (0)	0 (0)	1 (0.02)
	2014	I: 43U: 60	25 (0.18)	55 (0.27)	5 (0.04)	0 (0)	1 (<0.01)	0 (0)	1 (<0.01)
	2015	I: 40U: 54	9 (0.09)	61 (0.33)	6 (0.04)	0 (0)	1 (0.01)	1 (0.01)	1 (0.01)
T. speciosus	2013	I: 66U: 92	96 (0.27)	31 (0.13)	20 (0.11)	0 (0)	0 (0)	1 (<0.01)	1 (<0.01)
	2014	I: 219U: 208	326 (0.43)	89 (0.13)	74 (0.11)	1 (<0.01)	0 (0)	1 (<0.01)	2 (<0.01)
	2015	I: 145U: 190	346 (0.44)	137 (0.24)	56 (0.12)	1 (<0.01)	2 (<0.01)	4 (<0.01)	3 (<0.01)
T. alpinus		I: 115U: 145	0.12	0.34	0.04	<0.01	<0.01	<0.01	0.01
T. speciosus		I: 475U: 445	0.40	0.17	0.12	<0.01	<0.01	<0.01	<0.01

Table 3

Summary of counts for each flea species at each study site. The number of hosts is the total number of both chipmunk species that were combed for fleas. See Table 2 caption for flea species abbreviations.

Elev. (m)	Site	# Hosts	Ccm	Eem	Eet	Aw	Cat	Oi	Om	Years Visited	Hosts Present
2850–3020	Arrowhead Lake	189	176	65	46	0	0	0	0	2013–14	Ts
2800–3060	Cathedral Lake (upper)	166	82	96	18	1	1	4	4	2013–15	Ta, Ts
2380–2435	Glen Aulin	11	0	0	0	0	0	0	0	2013	Ts
3095–3230	Gaylor Lakes	50	27	32	2	0	1	0	1	2014–15	Ta, Ts
2585–2660	Hoffmann Creek	12	5	2	3	0	0	0	0	2013	Ts
2570–2710	Mammoth Mountain	189	81	6	17	0	0	1	1	2013–15	Ts
2690–2995	May Lake	242	185	129	45	2	2	0	3	2013–15	Ta, Ts
2330–2460	Porcupine Creek	16	3	0	0	0	0	0	0	2013	Ts
2945–3230	Saddlebag Lake	90	52	53	0	0	0	2	0	2013–15	Ta, Ts
2955–3210	Saddlebag Lake - North	122	143	53	0	0	0	0	0	2013–15	Ta, Ts
2545–2690	Tuolumne Meadows	93	51	0	34	0	0	0	0	2013–14	Ts

Fig. 2

Patterns of flea abundance across years, hosts, and flea species. Overall average flea abundances (A–B) and abundances of each flea species (C–D) in each year for T. alpinus (A, C) and T. speciosus (B, D). Abundances of each flea species on hosts of each sex (Males: closed circles, Females: open circles) on T. alpinus (E) and T. speciosus (F).

For both hosts, a large proportion of captured individuals did not have fleas (Fig. 3), and the majority of all fleas collected were carried by a disproportionately small number of hosts. For example, for both chipmunk species across all years the top 10% most flea-parasitized individuals accounted for approximately half of the total number of fleas observed (Fig. 3). These dynamics also changed across years; for example, in 2013, 2014, and 2015 T. alpinus individuals with no fleas represented 49.1%, 58.8%, and 57.1% of the population respectively, while for T. speciosus these percentages were 58.2%, 47.9%, and 42.5% (Fig. 3). This indicates that the average number of fleas per individual on T. speciosus increased across years due to both a larger prevalence of infested individuals (Fig. 3B), and higher flea abundances per infested individual (Fig. 2B).

Fig. 3

Cumulative distribution plots divided by year for (A) For each species 2013 is shown in red, 2014 in teal, 2015 in pink. The x-axis represents each host individual, ordered from least to most flea infested, and the y-axis shows the cumulative proportion of total flea counts. The dotted line indicates individuals in the 90th percentile of flea abundances, illustrating that the top 10% most infected chipmunks usually account for close to 50% of all counted fleas. The proportion of individuals without fleas in each year is represented graphically as the proportion at which each colored line departs from the x-axis. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Host biological predictors of flea abundance

For both hosts, males had higher mean flea abundance than females (T. alpinus: W = 10530, p = 0.001; T. speciosus: W = 116030, p = 0.009), a pattern that was most strongly driven by E. eumolpi abundances (Fig. 2E–F). Fecal glucocorticoid metabolites were significantly, negatively correlated with flea abundances for females but not males of both host species (Spearman's rank correlation rho, T. alpinus females: S = 419090, r = −0.23, p = 0.02; T. alpinus males: S = 318590, r = −0.03, p = 0.84; T. speciosus females: S = 19199000, r = −0.15, p = 0.003; T. speciosus males: S = 11952000, r = −0.003, p = 0.95; Fig. 4). Body mass was significantly but weakly negatively correlated with flea abundances for T. speciosus but not for T. alpinus (Spearman's rank correlation rho, T. alpinus: S = 2692800, r = −0.059, p = 0.48; T. speciosus: S = 127690000, r = −0.098, p = 0.009). Generalized linear mixed models revealed that these relationships between host biological factors and flea abundances depended on the flea species. Aforementioned relationships between sex, glucocorticoid metabolite levels, and flea abundances appeared to be driven by E. eumolpi; these predictors did not explain significant variance in models of C. ciliatus mononis abundance (Table 4, Table 5). After controlling for sex, GCs also had a marginally significant (p = 0.05) relationship with E. eumolpi abundance (Table 5).

Fig. 4

Relationships between fecal glucocorticoid metabolite levels, sex, and flea abundance for (A) Points show the mean±S.E. number of fleas counted for female (white) and male (black) individuals within FGM categories (FGM values were rounded to the nearest 10). Lines of best fit (based on all raw data points) ±95% confidence intervals are overlaid for each sex.

Table 4

Negative binomial mixed model results testing which factors are most strongly predictive of C. ciliatus mononis abundances. Significant terms are bolded.

Coefficient	Estimate	S.E.	z value	p-value
(Intercept)	−2.07	0.37	−5.65	<0.00001
GC	−0.11	0.07	−1.50	0.13
Sex (Male)	0.14	0.12	1.22	0.222
Species (Ts)	1.67	0.35	4.70	<0.00001
PC1_temp	−0.08	0.17	−0.44	0.66
PC2_temp	0.02	0.32	0.07	0.94
Density	0.01	0.07	0.1	0.86
Date	−0.16	0.06	−2.52	0.01
Elevation	−0.23	0.37	−0.63	0.53
GC*Sex (Male)	0.17	0.12	1.45	0.15
Species (Ts)* PC1_temp	0.14	0.18	0.76	0.45
Species (Ts)* PC2_temp	−0.05	0.33	−0.15	0.88
Species (Ts)* Elevation	0.58	0.37	1.56	0.12

Table 5

Negative binomial mixed model results testing which factors are most strongly predictive of E. eumolpi abundances. Significant terms are bolded.

Coefficient	Estimate	S.E.	z value	p-value
(Intercept)	−1.24	0.42	−2.96	0.003
GC	−0.29	0.15	−1.92	0.05
Sex (Male)	0.48	0.15	3.23	0.001
Species (Ts)	−0.77	0.26	−2.99	0.003
PC1_temp	−0.13	0.13	−1.07	0.28
PC2_temp	−0.46	0.19	−2.45	0.01
Density	0.07	0.13	0.56	0.58
Date	0.23	0.10	2.24	0.03
Elevation	0.31	0.33	0.95	0.34
GC*Sex (Male)	0.40	0.18	2.17	0.03
Species (Ts)* PC1_temp	0.36	0.15	2.49	0.01
Species (Ts)* PC2_temp	0.33	0.21	1.53	0.13
Species (Ts)* Elevation	0.91	0.27	3.33	0.0009

Environmental predictors of flea abundance

Overall flea abundances (across all flea species) were predicted by ambient temperature and elevation in a host-species-dependent manner (Table 6). Variance in temperature appeared to be the most salient, temperature-related environmental predictor. Tamias speciosus trapped on days and in areas with higher PC2 scores (which aligned with higher variance in daily temperature) tended to have more fleas; the opposite pattern was true for T. alpinus (Table 6; Fig. 5A). Hosts of both species caught in areas with higher PC1 scores (which aligned with higher mean and maximum daily temperatures), on the other hand, had slightly fewer fleas on average, but this pattern was less obvious (Table 6). Tamias speciosus trapped at higher elevations had more fleas, while the opposite pattern was found for T. alpinus (Table 6; Fig. 5B).

Table 6

Negative binomial mixed model results testing which factors are most strongly predictive of overall flea abundances. Significant terms are bolded.

Coefficient	Estimate	S.E.	z value	p-value
(Intercept)	0.14	0.27	0.51	0.61
GC	−0.03	0.05	−0.53	0.59
Sex (Male)	0.33	0.10	3.33	0.0009
Species (Ts)	0.08	0.20	0.40	0.70
PC1_temp	−0.20	0.10	−1.99	0.046
PC2_temp	−0.45	0.16	−2.90	0.004
Density	−0.008	0.07	−0.11	0.91
Date	0.03	0.07	−0.43	0.67
Elevation	−0.31	0.21	−1.49	0.14
GC*Sex (Male)	0.10	0.09	1.02	0.31
Species (Ts)* PC1_temp	0.26	0.11	2.38	0.02
Species (Ts)* PC2_temp	0.46	0.17	2.73	0.006
Species (Ts)* Elevation	0.61	0.21	2.96	0.003

Fig. 5

Relationships between flea abundances and (A) the second principal component of temperature data; or (B) elevation for Points show the mean±S.E. number of fleas counted for a given study site in a given year. For each study site in each year, a mean±S.E. temperature or elevation value is shown. Lines of best fit (based on all raw data points) ±95% confidence intervals are overlaid for each species.

Flea-species-specific models (Table 4, Table 5) revealed species-dependent differences in which of these environmental factors were most predictive of flea abundance. Models did not identify any significant environmental predictors of C. ciliatus mononis abundance, except for day of the year (C. ciliatus mononis abundances were slightly higher later in the season; Table 4). Eumolpianus eumolpi abundance followed similar patterns as overall flea abundances: T. speciosus hosted more E. eumolpi at higher elevations and in areas with more variable temperatures, whereas the opposite patterns were found for T. alpinus (Table 5).

Host population densities were never predictive of flea abundances (Table 4, Table 5, Table 6). Eumolpianus eumolpi abundance significantly increased as the season progressed. Site explained significant variance in E. eumolpi and overall flea abundance, and year explained significant variance in overall flea abundances; no other random effects were significant in these or other models (Table 4, Table 5, Table 6).

Discussion

Survey of fleas

We identified at least seven flea species from five genera hosted by the focal chipmunk species. Only one published account of fleas parasitizing T. alpinus exists in the literature, which documents a flea species we did not find (a prairie dog flea, Oropsylla tuberculata tuberculata; Fleer et al., 2011). In contrast, all flea species we identified have been previously collected from T. speciosus (Barnes and Kartman, 1960; Best et al., 1994; Adjemian et al., 2008; Fleer et al., 2011). Eumolpianus eumolpi was most prevalent on T. alpinus, whereas C. ciliatus mononis was the most prevalent species on T. speciosus, and E. eutamiadis was found in lower frequencies on both hosts. Flea abundances and the proportion of individuals carrying fleas increased from year to year for T. speciosus but not for T. alpinus. High flea abundances are associated with increased rates of plague transmission (Krasnov et al., 2006; Pham et al., 2009; Tripp et al., 2009), thus, the increase in flea abundance and prevalence in T. speciosus from 2013 to 2015 has potential human health consequences and may have been related to the 2015 plague epizootic. All three of these flea species can vector plague, can infest other common rodent hosts like tree squirrels, ground squirrels and deer mice, and will bite humans (Eskey and Haas, 1939; Barnes and Kartman, 1960; Fagerlund et al., 2001; Eisen et al., 2009; Smith et al., 2010; Fleer et al., 2011; Danforth et al., 2016). In a recent study in the northern Sierra Nevada and the Southern Cascade mountains of California, C. ciliatus mononis and E. eumolpi were found to be the most frequently Y. pestis positive fleas (Smith et al., 2010). In the 2015 plague epizootic in Yosemite, five of 27 tested T. speciosus individuals were Y. pestis positive, as were a pool of C. ciliatus mononis fleas (Danforth et al., 2016). Tamias speciosus is also known to reservoir Borrelia spp., Rickettsia spp., and other zoonotic pathogens at relatively high frequencies (Adjemian et al., 2008; Fleer et al., 2011; Straub et al., 2017). Follow-up work related to this study will test sampled fleas for infection with Y. pestis and other pathogens.

These survey results can also serve as a baseline to which future surveys may be compared. Because T. alpinus has exhibited more drastic spatial, genetic, dietary, and morphological responses to climate change in comparison to T. speciosus, it may be worthwhile to examine how flea communities differentially change on these hosts in the upcoming decades (Moritz et al., 2008; Rubidge et al., 2012; Walsh et al., 2016). One caveat is that our method for collecting fleas did not make use of a white tray or sheet to facilitate capturing fleas that were dislodged from the host and the comb, meaning that a larger proportion of fleas may have escaped the collection area. Future studies should take methodological differences like these into consideration, and it would be valuable for more studies to quantify detection error or recovery differences in flea collection methods (e.g. Eads et al., 2013).

Predictors of flea abundances

Biological factors were predictive of flea abundances for both hosts. Males tended to have higher flea abundances than females, a pattern apparently driven by E. eumolpi abundances. This is a commonly documented pattern, likely due in part to sexual differences in space use related to polygynous mating systems (Zuk and McKean, 1996; Schalk and Forbes, 1997, but see Kiffner et al., 2013). We also found sex-specific relationships between glucocorticoids and flea abundances. For both hosts, females (but not males) with higher glucocorticoids had lower flea abundances. Models of E. eumolpi abundance also suggested a negative correlation between abundance and GCs. While it is sometimes hypothesized that animals with higher glucocorticoids may be immunocompromised and therefore may harbor higher ectoparasite loads (Owen et al., 2010), studies explicitly testing for relationships between ectoparasitism and glucocorticoids have found variable results (Poiani et al., 2000; Lobato et al., 2008; Monello et al., 2010; St. Juliana et al., 2014). Our results suggest that higher glucocorticoids are associated with reduced ectoparasite abundances. These findings underscore the importance of taking sex-related effects into consideration in studies of parasitism.

Environmental predictors of flea abundances diverged between host and flea species. Surprisingly, host population density, which is often associated with flea densities in other systems (Arneberg et al., 1998; Krasnov et al., 2002; Young et al., 2015), was never predictive of flea abundances. Moreover, for C. ciliatus mononis we could not identify any environmental predictors of flea abundance. This may be a result of the large number of predictor variables included in models, which likely reduced statistical power. When non-significant variables were eliminated in a backwards, stepwise manner, it eventually emerged that C. ciliatus mononis abundances increased at higher elevations. However, stepwise model fitting can be statistically problematic, and was therefore not formally implemented in this study (Mundry and Nunn, 2008). What is clear is that E. eumolpi abundances were better predicted by the environmental variables we tested than C. ciliatus mononis abundances. Because C. ciliatus mononis seems to preferentially infest T. speciosus, which is more of a generalist species inhabiting a broad range of habitats in comparison to T. alpinus, it is possible that specific environmental factors are simply less critical to determining C. ciliatus mononis abundances, as they are for its host. However, it is also possible that other environmental variables that we were not able to measure, including relative humidity or soil characteristics, may better predict C. ciliatus mononis abundances than temperature-related variables. Moreover, we did not sample T. speciosus at the lower end of its elevational range, and thus likely did not sample the entire environmental niche of C. ciliatus mononis. This means that we may have not captured relevant environmental variance in temperature for this flea species; possibly when including lower elevations with warmer temperatures, different patterns might emerge.

We did capture the majority of T. alpinus' elevational range in our surveys, and thus may have collected data on a larger proportion of E. eumolpi's niche (as this species was more common on T. alpinus), possibly facilitating the detection of environmental patterns. For E. eumolpi variance in daily temperatures and elevation were negatively associated with abundances on T. alpinus, and positively associated with abundances on T. speciosus. Variance in daily temperatures was also associated with overall flea abundances (across all species), such that in more variable environments there were higher abundances for T. speciosus, and lower for T. alpinus. Variance in daily temperature appeared to be the strongest temperature-related driver of flea abundances in general; indeed, models with only variance in temperature usually performed as well as models using principal component axes that summarized all three temperature variables (variance, mean, and maximum).

The climate variability hypothesis for disease-related declines hypothesizes that, given their smaller body size and faster metabolism, pathogens and parasites may be able to acclimatize more quickly to variability in climate in comparison to their hosts (Rohr et al., 2010; 2013). In some climate change scenarios variability in temperature is expected to increase (Raisanen, 2002; Della-Marta et al., 2007; Fischer and Schär, 2009; Vázquez et al., 2017). It is unclear to what extent patterns of daily variance in temperature align with overall climate variability at different spatial and temporal scales, and even whether daily variance in temperature is expected to increase or decrease with climate change (e.g. Michaels et al., 1998; Easterling et al., 2000). However, daily temperature variance is known to impact disease dynamics in other systems (e.g. Paaijmans et al., 2010), and the fact that flea abundances appear to respond to temperature variance in both host species is therefore potentially worrisome in the context of climate change. That said, many other studies have predicted declines in parasites with projected changes in climate (e.g. Carlson et al., 2017; Gehman et al., 2018). Because environmental parameters can have shared or opposing effects on host and parasite physiology, and due to the trade-offs between a host's quality as a food source and its resistance to parasitism, these dynamics are complicated to predict (Hawlena et al., 2005). The fact that year was consistently a significant predictor for flea abundances indicates that there are factors we did not measure – for example, snowpack, which decreased significantly from year to year (Berg and Hall, 2017), or humidity, which is known to impact flea development – that may account for a significant proportion of the variation in flea abundance on this species. Clearly, the factors predicting C. ciliatus mononis remain poorly understood. In future studies it would be valuable to integrate measures of relative humidity, which is expected to change with climate change and which is known to impact fleas in species-specific manners (Krasnov et al., 2001; Kreppel et al., 2016). Because fleas often inhabit host burrows and are impacted by burrow conditions, it would also be valuable to quantify the impacts of soil type on these species.

Conclusions

Our findings underscore the fact that environmental predictors of flea abundances differ depending on host and flea species. In some cases, variance in climate, which has been indicated in other studies of parasitism and disease (Rohr et al., 2010; 2013), may be a strong positive predictor of flea abundance. While we cannot draw causal connections based on our data and the number of years is limited, our surveys document that flea abundances on T. speciosus and the proportion of T. speciosus carrying fleas – both of which could facilitate disease transmission (Krasnov et al., 2006; Pham et al., 2009; Tripp et al., 2009) – increased in the years of and leading up to a plague epizootic, suggesting the possibility that these metrics should be further examined for their potential as warning signs of a plague event. Our results also document significantly different flea communities in two co-occurring congeners, record T. alpinus's flea communities in detail for the first time in the literature, and provide an important ecological baseline to which future surveys may be compared.