# Drivers of flea abundance in wild rodents across local and regional scales in the Chihuahuan Desert, northwestern Mexico

Paulina Pontifes, Adriana Fernández-González, Gabriel García-Peña, Benjamin Roche, Gerardo Suzán

► To cite this version:

Paulina Pontifes, Adriana Fernández-González, Gabriel García-Peña, Benjamin Roche, Gerardo Suzán. Drivers of flea abundance in wild rodents across local and regional scales in the Chihuahuan Desert, northwestern Mexico. Ecosphere, 2022, 13 (4), pp.e4013. 10.1002/ecs2.4013. hal-04978794

HAL Id: hal-04978794

https://hal.science/hal-04978794v1

Submitted on 5 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Drivers of flea abundance in wild rodents across local and regional scales in the Chihuahuan Desert, northwestern Mexico

Paulina A. Pontifes<fnref n="1" />,<fnref n="2" /> ID

Adriana Fernández-González<fnref n="2" />

<figure><img src="image_2.png" /></figure>

Gabriel E. García-Peña<fnref n="3" /> ID

Benjamin Roche<fnref n="4" /> ID

Gerardo Suzán<fnref n="2" /> ID

## Correspondence

Gerardo Suzán
Email: gerardosuz@gmail.com

## Funding information

Programa de Apoyo a Proyectos de Investigación e Innovación Tecnológica (PAPIIT), Grant/Award Number: IN225219; National Autonomous University of Mexico; Consejo Nacional de Ciencia y Tecnología

Handling Editor: Andrea Swei

## Abstract

The broad distribution of macroparasites and their thriving populations are matters of health and economic concern. Macroparasites cause damage both directly through their feeding habits, which impact host fitness, and indirectly through the transmission of various infectious diseases of relevance to human and domestic animal health and wildlife conservation. Because the impacts of macroparasites on host health and the risk of disease transmission are directly related to their abundance, understanding the drivers of macroparasite burden is of relevance. Various host traits and environmental factors have been associated with differences in macroparasite abundance. In addition to these variables, spatial scale is increasingly incorporated to understand how these drivers vary across space. However, variation in the relative importance of host traits and environmental factors as predictors of abundance at different scales is not well understood. To further clarify the relationship between scale and drivers of macroparasite abundance, we investigated the effects of host traits and environmental factors on flea abundance in rodents of the Chihuahuan Desert in northwestern Mexico on three levels: within a single site, between sampling sites with different vegetation types, and across the region. This partitioning allowed us to compare drivers at both local and regional scales. Fleas provide a natural model to assess the interplay between host and environmental variables across scales because their life cycles alternate between on-host and off-host environments and their hosts have varying ranges of distribution. We sampled 1311 fleas from 674 rodent individuals of 14 different species across 40 sampling plots between 2012 and 2013. Using generalized linear mixed models, we found that flea abundance was associated with different combinations of host traits such as size and sex. The specific combination of predictive variables differed across species, while the effects on flea abundance showed context and scale dependency, although this could only be tested at the full level of analysis on the most abundant species, *Dipodomys merriami*.

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

© 2022 The Authors. Ecosphere published by Wiley Periodicals LLC on behalf of The Ecological Society of America.

<figure><img src="image_3.png" /></figure>

<footnote n="1">¹Posgrado en Ciencias Biológicas, Facultad de Ciencias, Universidad Nacional Autónoma de México, Ciudad Universitaria, Ciudad de México, Mexico</footnote>

<footnote n="2">²Departamento de Etología, Fauna Silvestre y Animales de Laboratorio, Facultad de Medicina Veterinaria y Zootecnia, Universidad Nacional Autónoma de México, Ciudad Universitaria, Ciudad de México, Mexico</footnote>

<footnote n="3">³Centro de Ciencias de la Complejidad C3, Universidad Nacional Autónoma de México, Ciudad Universitaria, Ciudad de México, Mexico</footnote>

<footnote n="4">⁴MIVEGEC (Infectious Diseases and Vectors: Ecology, Genetics, Evolution and Control), IRD (Research Institute for Sustainable Development), CNRS (National Center for Scientific Research), University of Montpellier, Montpellier, France</footnote>
Sampling season was the only variable consistently significant across scales, reflecting the far-reaching effects of large-scale, interannual environmental fluctuations. These results emphasize that integrating spatial scale can strengthen study design for monitoring macroparasite burden.

### KEYWORDS

Chihuahuan Desert, flea abundance, macroparasites, mixed models, rodents, spatial scale

## INTRODUCTION

Macroparasites can reduce host fitness directly through their feeding habits and indirectly through the transmission of various infectious diseases of relevance to human and domestic animal' health and wildlife conservation (Wall & Shearer, 1997). The negative impact of macroparasites on their hosts is related to their abundance: Higher burdens result in more consumption of host resources and increased risk of disease transmission in the case of vector macroparasites (Bethony et al., 2006; Eisen & Gage, 2012). Therefore, understanding the drivers behind macroparasite abundance has important epidemiological implications.

Various environmental factors and host traits have been extensively studied as drivers of macroparasite abundance in vertebrates, with a few patterns reported across different systems. For example, temperature and humidity can modify macroparasite development and mortality (Wilson et al., 2002). Host sex influences macroparasite abundance through differences in host behavior, size, and immune-modulating effects of sex hormones (Zuk & McKean, 1996), with males typically being more parasitized than females (Skorping & Jensen, 2004). However, even widely observed patterns are not universal and many studies report variation in the role of host traits and environmental variables between different host-parasite systems and even within the same systems depending on location (Kiffner et al., 2013, 2014). Furthermore, other studies have found that the effects and relative importance of host traits and environmental variables on macroparasite abundance might be inconsistent across space (Cardon et al., 2011; Young & Maccoll, 2016). Thus, while extensive research has assessed the role of host traits and environmental variables on macroparasite abundance, less is understood about how the effect and relative importance of these drivers change when considered across spatial scales.

Studies performed at a local scale capture small-scale variability in the host's biotic and abiotic environment. However, hosts are subjected to spatiotemporal variation in environmental conditions across their distribution range (Penczykowski et al., 2016). This variation may

influence macroparasite abundance directly through differences in environmental conditions such as soil humidity or temperature across sites (Krasnov et al., 2001, 2002a) or indirectly through effects on resource availability or host behavior (Khokhlova, 2004; Ostfeld et al., 2006). Furthermore, spatiotemporal variation in individual-level (sex, age, size, and reproductive status) and population-level (density, age distribution, and sex ratio) host traits will also differ across the host's landscape. Thus, drivers of macroparasite abundance show dependence on the host's spatiotemporal context.

Incorporating spatial scale into the analysis of drivers of macroparasite abundance can help clarify their effects across scales and assess their relative importance across levels of ecological organization. For example, Young et al. (2015) found that environmental variables and host traits driving flea abundance in small mammals of the East African savanna had a higher predictive power across species but were not significant at the individual level. Linardi and Krasnov (2013) found that at lower hierarchical levels (between individuals), flea and mite abundance was affected by host and parasite traits and environmental factors (although effects differed between flea and mites), whereas at the higher levels (communities across a landscape), host traits and environmental variables drove variation in ectoparasite abundance.

In this study, we assess whether host traits and environmental variables have the same relative importance and effect on macroparasite abundance across spatial scales. We use fleas in rodent communities in a natural reserve within the Chihuahuan Desert in northwestern Mexico as a study system. Fleas (order Siphonaptera) provide a natural model to assess the interplay between host and environmental factors at different scales, as their life cycle alternates between on-host and off-host environment, requiring them to cope with the host's immune and behavioral responses, as well as with varying degrees of environmental exposure (dependent on the host's burrowing/nesting habits and range) (Krasnov, 2008). The broad distribution of some rodent species of the Chihuahuan Desert will allow us to assess how variation in host traits and environmental characteristics contributes to flea abundance at a local and a regional scale,
across different levels of host traits: individuals within a single site, between sites, and across host species within a region. Furthermore, by comparing flea abundance in sympatric rodents with a range of burrowing habits (fossorial, semi-fossorial, or shallow), we can explore the role of variation in microclimate conditions on flea abundance within and between host species. These results will enhance our understanding of host and environmental effects on flea abundance across spatial scales and highlight the importance of incorporating spatial context to accurately assess the effects of drivers of macroparasite abundance.

# METHODS

## Study region

The study was conducted in northwestern Chihuahua, Mexico, in the Janos Biosphere Reserve (JBR), a nature reserve located in the Chihuahuan Desert. The reserve covers an extent of 5264.9 km², with mosaics of grassland and shrubland vegetation interspersed with patches of agricultural land and human settlements (CONANP, 2013). The dominant climate in this region is temperate and semiarid, with an annual average of 381 mm of rain, 77% of which falls between April and August (CONANP, 2013). Temperature varies seasonally, with annual fluctuations of over 14°C.

## Rodent sampling

We sampled nocturnal rodents between 2012 and 2013 over three sampling sessions (May and October 2012, and May 2013), using Sherman traps (H.B. Sherman 8 × 8 × 23 cm, Tallahassee, FL) at four sampling areas (either MV, RO, EC, or PV; see Figure 1). At each sampling area, ten 7 × 7 grids with a 10-m spacing between traps were set. Sampling plots were located approximately 700–900 m apart from each other. Rodent movement between grids (assessed by the presence of individual rodents in multiple grids within a sampling area) was monitored to ensure sampling plot independence. Traps were baited with a mixture of oats and vanilla extract and set for three consecutive nights at each trapping grid, yielding 147 trap nights per site. After use, each trap was cleaned with hospital-grade detergent.

Captured rodents were identified to species level using taxonomic keys (Anderson, 1972; Reid, 2006). Body mass, length, and sex of host were recorded. Reproductive status was established upon observation of descended testes for males and perforated vagina or pregnancy

<figure><img src="image_4.png" /><figcaption>**FIGURE 1** Map of the study sites within the Janos Biosphere Reserve (red polygon outlines the study region) in the Chihuahuan Desert. Sampling was conducted across a 1000-km² area, considering a total of 40 sampling locations distributed across four sampling areas, to represent the local and regional ranges of ecological conditions of the study region. Sampling areas correspond as follows: Monte Verde (MV), Rancho Ojitos (RO), Pancho Villa (PV), and El Cuervo (EC)</figcaption></figure>

(determined by abdominal palpation) or lactation signs for females (Gurnell & Flowerdew, 2006). Sampled rodents were ear-tagged to avoid resampling. Most animals were released at the point of capture, although some were euthanized for morphological voucher specimens to verify identification. Procedures for trapping and handling were approved by the Animal Care Committee of the National Autonomous University of Mexico (UNAM) and by the Secretariat of Environment and Natural Resources (license number FAUT-0250) and followed the standards set by the American Society of Mammologists (Sikes & Gannon, 2011).

## Flea sampling

Each captured rodent was placed in a plastic chamber with a cotton ball dosed with isoflurane. This method anesthetizes both the host and its fleas, which dislodge from their hosts (Himsworth et al., 2020). Animals in the chamber were monitored to remove the lid and the cotton ball as soon as motor activity nearly ceased, in order to reduce the risk of death following anesthesia. Each anesthetized animal was held stretched and thoroughly combed within the chamber for 2 min with a standardized number of passes to collect fleas that had not fallen
off during anesthesia. Fleas were collected from the container by hand and stored separately in microtubes containing 70% ethanol and kept at -80°C until identification. The plastic chamber was cleaned with water and detergent after each use. Fleas were identified to species level using a dissecting scope (SZx12 Olympus, Melville, NY) and taxonomic keys (Acosta & Morrone, 2003; Hubbard, 1974). Given the purposes of the research, all analyses were conducted using total flea burden.

## Environmental characteristics and sampling season

Microclimate conditions in different vegetation types, even within the same region, can lead to differences in soil humidity and air temperature, parameters that have been shown to affect flea development and mortality (Krasnov et al., 2001, 2002a, 2002b). To explore the potential role of differences in microclimate on flea abundance, we recorded vegetation type at each sampling plot as either grassland or shrubland (considering the classification of previous analyses at the same sampling plots as part of ongoing research in the study area [see Rubio et al., 2015 for details]). Sampling session and year were combined as a categorical variable with three levels, corresponding to spring (May 2012 or May 2013) or autumn (October 2012) fieldwork sessions, to account for temporal variation between sampling sessions.

## Statistical analysis

We examined three levels of analysis: across individuals within a single site, between sites, and across species within the region (see Figure 1 for a map of the study region and see Appendix S1: Figure S1 for a visual representation of the analysis at each scale). Flea abundance was considered as the total number of individual fleas per host, expressed as either the absolute number of fleas or the mean number at the corresponding level (site or species). We considered the number of fleas per individual to represent the success of the fleas on the host once established (rationale discussed in Appendix S1). Although we report flea abundance for all host species (Appendix S1: Tables S1 and S2), statistical analyses were only conducted for host species where sample size *n* > 30, with at least three individuals in each of the levels of the factor variables. In addition to the previous criteria, only species that were present in at least half of the sampling locations in each sampling area were considered for analysis between sites.

All analyses described in the following sections were conducted using R v. 4.0.0 (R Core Team, 2020). Mixed

models were implemented with the package glmmTMB (Brooks et al., 2017). All statistical analyses were restricted to nonpregnant adults to avoid confounding effects of weight gain and loss associated with pregnancy and growth (Appendix S1: Table S3). Separate models were run for each species at the local and between-site levels. Fixed effects for each level are described in Appendix S1: Table S4. Collinearity between explanatory variables in the final models was assessed using a variance inflation factor test (VIF < 2). Continuous fixed effects were mean-centered prior to analysis. The intraclass correlation coefficient (ICC) was calculated as an indicator of variation in flea abundance due to differences within and between sampling plots, as a proxy to assess the role of unspecified environmental variation associated with local conditions of the sampling plots, either within a single site or across the host's landscape.

We used Nakagawa's and collaborators' $R^2$ as an estimator of the proportion of variance accounted for by the final models (Nakagawa et al., 2017), implemented in the package "MuMln" (Barton, 2013). This metric distinguishes between variance due to fixed effects ($R^2_{marginal}$) and variance conditional both on fixed and on random effects ($R^2_{conditional}$). To estimate the relative importance of the different variables on flea parasitism across scales, we used Akaike's information criterion weights (AICw), considering only the subset of models with $\Delta AIC_c < 2$ for model averaging and standardized predictor variables (Schielzeth, 2010). Model averaging to obtain the relative importance of predictors was performed using the MuMln package. Note that model inference is not based on averaged coefficients but rather on a single competitive model (as determined by model selection) (Cade, 2015). The full set of models for each level of analysis is presented in Appendix S1: Tables S5-S7.

## Across individuals within a locality

The number of fleas per host was modeled with either a negative binomial or a Poisson error structure (Appendix S1: Table S6). Only sampling locations from a single area (RO in Figure 1) with the same habitat type across sampling plots (shrubland) were considered for analysis at this level. Fixed effects included sex, reproductive status, body size, weight, and sampling season, with sampling plot as a random effect. Previous to statistical modeling, we assessed the correlation between morphometric variables (weight and body length) (Appendix S1: Figure S2). Model comparison was performed by backward stepwise elimination of nonsignificant terms (*p* < 0.05) from a maximal model that considered all terms and plausible biological interactions (Appendix S1: Table S5). The significance of the variables
and their interactions was evaluated using likelihood ratio tests (LRTs). Further support for variables in the final model was provided by inspecting their relative importance according to model averaging results (Appendix S1), although model coefficients in result tables are presented for a single competitive model (lowest AIC value). Model diagnostics and checks for overdispersion and zero inflation were conducted using the package DHARMa (Hartig, 2016) to ensure final models did not violate any assumptions (Zuur et al., 2010).

## Across sites

The response variable at this level of analysis was the mean abundance of fleas per sampling location. Fixed effects included habitat type at location (either shrubland or grassland) and the morphometric variables explored at the previous level, scaled appropriately. Note that data from the two sampling sessions in 2012 were pooled after checking for differences in mean abundance of fleas (see Appendix S1 for details), so the categorical variable sampling season variable has two levels. Finally, terms to represent conspecific host abundance or abundance of demographic subgroups (males, females, and reproductive or nonreproductive individuals) were also included. While host abundance has been associated with effects of macroparasite abundance on a theoretical and empirical level in some systems (Anderson & May, 1978; Stanko et al., 2002), we wanted to assess the role of specific subgroups on flea abundance, as demographic structure might be key to identifying drivers of abundance within a population (see, e.g., Perkins et al., 2003). Separate sets of models were run for each demographic subgroup. Additional methods, results, and further description of the rationale for analysis at this level are included in Appendix S1. We included sampling location as a random effect to assess the contribution of within-location and between-location variation. Model comparison for each set was conducted as described in the previous section. It is important to mention that although we only considered first capture individuals for our analysis, we monitored recaptures to assess the movement of individuals between sampling plots.

## Across host species within region

To assess drivers of flea abundance across species within the region, we used average flea abundance as a response variable. This response variable was modeled with a Gaussian distribution and identity link function. Fixed effects included sampling season, body size, and mass,

while sampling area was considered a random factor. Only results for the best-fit model are shown. Additionally, we assessed the role of host identity on flea abundance to evaluate whether certain host species were associated with higher flea abundance. For these models, we used total flea counts on individuals as a response variable, modeled with negative binomial distribution. Only species with at least *n* > 30 individuals were considered for analysis at this level. Fixed effects were included to control for variation associated with sampling season, vegetation type, and sex-related biases. Morphometric variables were excluded due to high collinearity (VIF > 2) with host identity. Additionally, sampling area and location were considered random effects to assess the proportion of variance explained by spatial differences within and between sampling sites within areas.

# RESULTS

## Fieldwork

We captured and sampled a total of 674 rodents belonging to 14 species across three families (Cricetidae, Heteromyidae, and Sciuridae) (Appendix S1: Table S1). A total of 1311 fleas were collected from sampled rodents (spring 2012: 98 fleas; autumn 2012: 400 fleas; and spring 2013: 813 fleas). Summary information and details of the fleas found are presented in Appendix S1: Table S2. Three host species, *Dipodomys merriami*, *D. spectabilis*, and *Onychomys arenicola*, represented 66% of total individuals sampled, 80% of individuals with at least one flea, and between 37% and 92% of total individuals per sampling locations. We did not trap any same rodents on different sampling grids throughout our fieldwork, indicating that the separation between our sampling plots was an adequate representation of the maximum movement distance of the sampled species.

## Across individuals within a locality

Three species met the criteria for analysis at this level: *D. merriami*, *D. spectabilis*, and *O. arenicola*. The negative binomial distribution offered the best fit for *D. merriami* and *D. spectabilis* flea abundance data, while *O. arenicola* was best modeled by a Poisson distribution (Appendix S1: Table S6). Sampling season was the only significant and important variable across the three species, with the strongest effects observed during the sampling season corresponding to May 2013 in two of the three species analyzed (Table 1). Different effects of host traits on flea abundance were observed in *D. merriami* and *O. arenicola*; in the latter, flea abundance was higher in
individuals of smaller body size (−0.02 [−0.04, −0.001 log units]). An interaction between sex and size was detected in *D. merriami*, with larger-than-average females presenting higher flea abundance than their male counterparts (Figure 2). This interaction was significant, as assessed by a LRT ($\chi^2 = 7.75$, df = 1, *p* = 0.005), with a relative importance greater than 0.5 (Table 1). However, the effect size on the response variable was small (0.03 [0.01, 0.05] log units). Interestingly, while the ICC value for *D. merriami* and *O. arenicola* was low (<0.001), indicating larger within-grid variance than between-grid variance, ICC for *D. spectabilis* had a value of 0.57, representing a moderate correlation between observations in the same grid. The relevance of the random effects in *D. spectabilis* is further supported by the total variance explained by the models, where $R^2_c$ accounted for 46% of the variation in flea abundance, in contrast to the proportion accounted for by $R^2_m$ (27%). The marginal variance explained by the models for *D. merriami* (36%) or *O. arenicola* (65%) did not change for $R^2_c$.

## Across sites

Analysis at this level was only conducted for one species, *D. merriami*, as it met the selection criteria outlined in the Methods section. Mean flea abundance in this species across sites was negatively associated with male abundance, while the third sampling season had a markedly

positive effect (Table 2). An interaction between vegetation and body size, where shrubland sites with larger-than-average individuals were associated with a higher mean flea abundance (compared to grassland sites with larger-than-average individuals), was significant (LRT: $\chi^2 = 11.89$, df = 1, *p* = 0.001). All fixed effects had high relative importance. The among-grid variance was larger than within-grid variance (Table 2). Fixed effects, according to the best-fit model, accounted for 79% of variation in the response variable, which increased when considering the conditional $R^2$ (86%).

## Across host species within region

Sampling season was identified as the single most significant and important variable associated with variation in mean flea abundance across species. Specifically, the last sampling season increased mean flea abundance across species by 2.79 (1.01, 4.58) (Table 3). Sampling season accounted for 19% of variation in mean flea abundance at this level of analysis, a moderate result that did not increase when considering the effect of the sampling site. However, the within-site variance was very high ($\sigma^2 = 6.7$). Models assessing the role of host identity on flea abundance showed that specific host species were associated with higher abundance (Appendix S1: Table S8). Sampling season had a marked effect on flea abundance across species (Figure 3).

<table><caption>TABLE 1 Variables explaining flea abundance across individuals of *Dipodomys merriami*, *D. spectabilis*, and *Onychomys arenicola* at a local scale, according to generalized linear mixed model results</caption><thead><tr><th>Host species</th><th>Fixed effect</th><th>Estimate</th><th>SE</th><th>p</th><th>RI</th><th>σ<sup>2</sup></th><th>τ<sub>00</sub></th><th>ICC</th><th>R<sup>2</sup><sub>m</sub></th><th>R<sup>2</sup><sub>c</sub></th></tr></thead><tbody><tr><td rowspan="5">D. merriami</td><td>Intercept</td><td>-1.97</td><td>0.46</td><td><strong>&lt;0.001</strong></td><td>1.00</td><td>0.82</td><td>&lt;0.001</td><td>&lt;0.001</td><td>0.36</td><td>0.36</td></tr><tr><td>Season 2</td><td>0.22</td><td>0.59</td><td>0.71</td><td>1.00</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Season 3</td><td>3.05</td><td>0.48</td><td><strong>&lt;0.001</strong></td><td>1.00</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Size</td><td>-0.02</td><td>0.01</td><td><strong>0.03</strong></td><td>0.83</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Sex (F)</td><td>-0.07</td><td>0.19</td><td>0.71</td><td>0.83</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="4">D. spectabilis</td><td>Size : Sex (F)</td><td>0.03</td><td>0.01</td><td><strong>0.005</strong></td><td>0.83</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Intercept</td><td>-0.70</td><td>0.59</td><td>0.23</td><td>1.00</td><td>0.51</td><td>0.67</td><td>0.57</td><td>0.27</td><td>0.46</td></tr><tr><td>Season 2</td><td>1.72</td><td>0.51</td><td><strong>0.001</strong></td><td>1.00</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Season 3</td><td>2.64</td><td>0.52</td><td><strong>&lt;0.001</strong></td><td>1.00</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="4">O. arenicola</td><td>Intercept</td><td>1.21</td><td>0.27</td><td><strong>&lt;0.001</strong></td><td>1.00</td><td>0.31</td><td>&lt;0.001</td><td>&lt;0.001</td><td>0.65</td><td>0.65</td></tr><tr><td>Season 2</td><td>-0.82</td><td>0.36</td><td><strong>0.02</strong></td><td>1.00</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Season 3</td><td>0.40</td><td>0.33</td><td>0.49</td><td>1.00</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Size</td><td>-0.02</td><td>0.01</td><td><strong>0.03</strong></td><td>1.00</td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

Note: Parameter estimates and standard errors (SE) shown correspond to the single best fixed-effects model as measured by the lowest AIC and likelihood ratio test comparison. Other models within 2 ΔAIC are reported in Appendix S1: Table S7.

Abbreviations: $\sigma^2$, within-group variance; $\tau_{00}$, between-group variance; ICC, intraclass correlation coefficient; $R^2_c$, conditional $R^2$; $R^2_m$ = marginal $R^2$; RI, the relative importance of the variable across models.
<figure><img src="image_5.png" /><figcaption>**FIGURE 2** Plots showing model predictions for flea abundance in individuals of *Dipodomys merriami*, *D. spectabilis*, and *Onychomys arenicola* at a local scale (single site). Variables shown correspond to those in the top model for each species, as determined by model comparison. Each plot shows the marginal effects with the nonfocal variables held constant. Plots were made using the package *ggeffects* (v1.1.1; Lüdecke, 2018). *Note*: The plot showing the interaction between size and sex in *D. merriami* (M for male and F for female) shows predicted results only for the third sampling season</figcaption></figure>

DISCUSSION

In this study, we assessed how spatial scale affects the role of host and environmental factors as drivers of flea abundance. We found that the effect and relevance of host traits differed across species and across scales, highlighting a dependence on the host's environmental context (Tables 1 and 2). We also found evidence for consistent effects of large-scale factors (sampling season) on flea abundance across species (Table 3) and across scales. While specific host species were associated with higher flea abundance (Appendix S1: Table S8), variation within and between sampling sites indicates an important role of local-scale variability, although their contribution toward accounting for flea abundance depended on the host species and the scale of analysis (Tables 1 and 2). Overall, these results underline that drivers of flea abundance, particularly those associated with host traits, exhibit variation across scales.

Across individuals within a locality

Host traits at a local scale (across individuals within a locality) were predictive of flea abundance in two of the three species analyzed. Specifically, in *D. merriami*, larger-than-average females presented higher flea abundance than larger-than-average males (Figure 2). Interestingly, larger size and male hosts are more frequently associated with higher parasite burdens in rodents and other vertebrates (Eads et al., 2020; Johnson & Hoverman, 2014). However, given the lack of size dimorphism in *D. merriami* and no effects of reproductive status, the interaction between sex and size in this species might be related to behavioral differences between males and females. Indeed, patterns of movement in this species, which are wider in males, have been proposed as a

<table><caption>**TABLE 2** Variables explaining mean flea abundance in *Dipodomys merriami* across sites, according to generalized linear mixed model results</caption><thead><tr><th>Host species</th><th>Fixed effect</th><th>Estimate</th><th>SE</th><th>p</th><th>RI</th><th>σ<sup>2</sup></th><th>τ<sub>00</sub></th><th>ICC</th><th>R<sup>2</sup><sub>m</sub></th><th>R<sup>2</sup><sub>c</sub></th></tr></thead><tbody><tr><td rowspan="6">D. merriami</td><td>Intercept</td><td>1.53</td><td>0.38</td><td><strong>&lt;0.001</strong></td><td>1.00</td><td>0.04</td><td>0.13</td><td>0.79</td><td>0.79</td><td>0.86</td></tr><tr><td>Season 3</td><td>2.49</td><td>0.25</td><td><strong>&lt;0.001</strong></td><td>1.00</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Size</td><td>-0.06</td><td>0.04</td><td>0.07</td><td>1.00</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Vegetation (shrubland)</td><td>-0.54</td><td>0.33</td><td>0.10</td><td>1.00</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Male abundance</td><td>-0.22</td><td>0.05</td><td><strong>&lt;0.001</strong></td><td>1.00</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Size : vegetation (shrubland)</td><td>0.16</td><td>0.04</td><td><strong>&lt;0.001</strong></td><td>1.00</td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

*Note*: Parameter estimates and standard errors (SE) shown correspond to the single best fixed-effects model as measured by the lowest AIC and likelihood ratio test comparison. Other models within 2 ΔAIC are reported in Appendix S1: Table S7.

Abbreviations: σ², within-group variance; τ₀₀, between-group variance; ICC, intraclass correlation coefficient; RI, the relative importance of the variable across models.
<table><caption>TABLE 3 Variables explaining mean flea abundance across species, according to generalized linear mixed model results</caption><thead><tr><th>Host species</th><th>Fixed effect</th><th>Estimate</th><th>SE</th><th>p</th><th>RI</th><th>σ<sup>2</sup></th><th>τ<sub>00</sub></th><th>R<sup>2</sup><sub>m</sub></th><th>R<sup>2</sup><sub>c</sub></th></tr></thead><tbody><tr><td rowspan="3">Across species</td><td>Intercept</td><td>1.03</td><td>0.64</td><td>0.11</td><td>1.00</td><td>6.7</td><td>&lt;0.001</td><td>0.19</td><td>0.19</td></tr><tr><td>Season 2</td><td>0.50</td><td>0.91</td><td>0.59</td><td>1.00</td><td></td><td></td><td></td><td></td></tr><tr><td>Season 3</td><td>2.79</td><td>0.91</td><td><strong>0.002</strong></td><td>1.00</td><td></td><td></td><td></td><td></td></tr></tbody></table>

Note: Parameter estimates and standard errors (SE) shown correspond to the single best fixed-effects model as measured by the lowest AIC. Other models within 2 ΔAIC are reported in Appendix S1: Table S7.

Abbreviation: σ², within-group variance; τ₀₀, between-group variance; ICC, intraclass correlation coefficient; RI, the relative importance of the variable across models.

<figure><img src="image_6.png" /><figcaption>FIGURE 3 Model predictions for flea abundance across sampling seasons in host species with *n* > 30 individuals. The third sampling season was associated with higher flea abundance across all of the analyzed species. Between-species comparison also shows that model predictions indicate a lower overall flea abundance for species such as *Chaetodipus penicillatus* and *Dipodomys merriami*, with more variable burdens for *D. spectabilis*, *Onychomys arenicola*, and *Peromyscus maniculatus*. 2012_1 = spring 2012; 2012_2 = autumn 2012; and 2013_1 = spring 2013. Each plot shows the marginal effects with the nonfocal variables held constant. The predicted effects were estimated using the ggpredict and the ggplot functions in the ggeffects package (Lüdecke, 2018)</figcaption></figure>

mechanism to escape parasitism pressure (Behrends et al., 1986). The size was also associated with variation in flea abundance in *O. arenicola*, with larger individuals harboring fewer fleas than smaller ones, a pattern hypothesized to correspond to better defenses against parasites in larger individuals (Kiffner et al., 2013; Sheldon & Verhulst, 1996).

Sampling season was the only variable associated with flea abundance in *D. spectabilis*. Indeed, sampling season had a large effect on individual-level flea abundance across all three species analyzed (Table 1), with the third sampling season, corresponding to spring 2013,

showing the most positive effect, whereas spring 2012 had an overall lower mean flea abundance (Figure 3), indicating a strong effect of annual variation over seasonal variation. The first period of our fieldwork coincided with the end of the most severe drought on record in northwestern Mexico (Murray-Tortarolo & Jaramillo, 2019), generating large-scale conditions of low humidity and high temperatures, which have been observed to decrease larval survival significantly and induce desiccation in adult fleas (Krasnov et al., 2002a, 2002b). Thus, under drought conditions, we would expect to see lower flea abundance, as observed in this study, particularly in rodent species with shallow burrows, which are more exposed to aboveground conditions.

Variation between sampling plots encompasses differences in variables such as substrate type or vegetation cover, which in turn affect parameters such as soil humidity and temperature that are relevant to flea development and survival (Krasnov, 2008). Although our study did not measure these directly, our models did indicate a role for variation between sampling plots at a local scale in *D. spectabilis* (Table 1) whose ICC shows a moderate correlation between observations from the same sampling plot. This implies that conditions within plots are more similar than conditions in other plots, even within the same sampling area. Such small-scale effects in plot-to-plot variation might be more relevant for flea abundance of fossorial or semi-fossorial species such as *D. spectabilis*, as small variation in these conditions could influence the burrow's microenvironment (Kay & Whitford, 1978).

## Across sites

The interaction between size and vegetation type observed at this level indicates that, while body size is a consistent predictor of flea abundance in *D. merriami* at both local and across-site levels, its effect is context-dependent: Mean flea abundance was higher at sites with larger-than-average specimens only at sites with shrubland, with the reverse pattern at grassland sites (Table 2). Increased near-surface air temperature at nighttime
has been reported in shrubland as compared to grassland vegetation in the Chihuahuan Desert (D'Odorico et al., 2010), resulting in local warming effects that could increase flea growth or reproductive parameters (Krasnov, 2008). This effect could be particularly marked on hosts with shallow burrows, such as *D. merriami*, where there are no significant differences between burrow and ambient atmosphere (Burda et al., 2007). Interestingly, while between-site variance was not relevant in *D. merriami* at a local scale (which included only sites with shrubland vegetation) (Table 1), across-site models for *D. merriami* (considering multiple locations with either shrubland or grassland vegetation) show that flea abundance at this level is affected by variation between sites (Table 2). Sampling season was also associated with flea abundance, showing the same effects as at the local scale. Finally, mean flea abundance in *D. merriami* across sites was negatively associated with male host abundance in this species (Table 2). While host abundance and macro-parasite abundance have long been known to correlate (Anderson & May, 1978), the demographic structure of host species populations is seldom considered, despite a potential role as drivers of parasite abundance; for example, large and sexually mature males of *Apodemus flavicollis* were found to drive *Ixodes ricinus* tick abundance (Perkins et al., 2003). Although longitudinal data would be required to understand the dynamics of rodent populations, our results suggest that higher abundance of male *D. merriami* could be associated with lower flea abundance, in line with other results that have found that specific demographic subgroups can drive ectoparasite abundance.

## Across host species within the region

No host traits were found to be associated with mean flea abundance across species within the study region (Table 3). Our results at this scale differ from those obtained from a similar study with rodents in Africa (Young et al., 2015), where the authors found that body mass accounted for a large proportion of the variation in mean flea burden across species. Additionally, the small amount of variation explained by models at this scale contrasts with the moderate-to-high variation accounted for by models at the previous scales (Tables 1 and 2). Furthermore, model results indicate high within-sampling-area variance (Table 3), suggesting that heterogeneities in the environmental context of the host, or in individual or across-site host traits, are important determinants of flea abundance, which are not represented by a high-level pooling of data at a regional scale.

In terms of host identity effects, three species were associated with higher flea abundance: *D. spectabilis*, *O. arenicola*, and *Peromyscus maniculatus* (Appendix S1: Table S8). Different host species differ in behavioral characteristics such as territoriality, which modifies movement and contact patterns, and burrowing habits, which fundamentally define flea exposure to the environment, which is, in turn, affected by variation in the microenvironmental conditions of the host's surroundings (Lareschi & Krasnov, 2010). Despite the diversity of life-history traits, sampling season had significant and consistent effects across all species studied. Interestingly, although there is a marked effect of sampling season on flea abundance across species, abundance in rodent species with shallow or simple burrows such as *D. merriami* and *P. maniculatus* reaches lower values during the sampling season associated with direct and post-drought effects (sampling seasons 1 and 2) according to model predictions (Figure 3). It is also interesting to note that, despite reports of higher flea abundance in spring and summer months as compared to fall or winter months in the ecoregion (López-Pérez et al., 2018), flea abundance during the second sampling season (autumn 2012) was four times higher than abundance in the first sampling season (spring 2012). Interannual season comparison shows that flea abundance during the third sampling season (corresponding to spring 2013) was eight times higher than sampling during spring 2012. While we note that the extraordinary climatic conditions associated with drought might not be representative of typical year-to-year seasonal variation in the study system, the significance of sampling season effects highlights that temporal variation, particularly if driven by large-scale interannual fluctuations, has far-reaching effects on flea abundance across all species included in the analysis.

## CONCLUSIONS

In agreement with current knowledge of drivers of flea abundance, we found that both host and environmental variables, with specific combinations differing among host species, drive flea abundance. However, we found that drivers of flea abundance also varied across space, highlighting the context dependency of host traits at local spatial scales, and the far-reaching effects of large-scale annual fluctuations. Indeed, although significant, the effect size of host traits on flea abundance was small, particularly at the local scale, and most noticeable on the third sampling season for some variable combinations (Figure 2). Thus, through its strong influence on flea abundance, environmental variation (associated with drought conditions in the case of our study) could be
modifying the effect of host traits on flea abundance. Long-term empirical data would be essential to determine whether the effects we observed vary in magnitude across years, particularly because our data were collected during nonstandard conditions. However, this does not preclude interest in the result that large-scale environmental fluctuations can have significant effects on flea abundance regardless of the scale of analysis and the species (with the important caveat that none of the species we analyzed have deep burrows that could buffer against drastic changes in aboveground conditions). Indeed, changes in weather patterns and in the frequency of extreme weather events worldwide due to climate change are expected to modify both macroparasite abundance and disease transmission (see, e.g., Eads et al., 2016).

Although the generality of our findings to other macroparasite systems would require further empirical research accounting for differences in transmission mode, dispersal capabilities, and a more comprehensive range of habitats, our results indicate that investigations of macroparasite abundance need to consider both large-scale fluctuations in environmental conditions and context dependence of effects of host traits across scales to accurately assess the relative importance of the factors that affect flea abundance and even disease dynamics. In this regard, a study by Ben Ari et al. (2011) analyzed the limitations of assuming scale independence and linearity in drivers of plague dynamics at large scales and found that effects at a given scale cannot be accurately extrapolated from effects observed at smaller scales, further highlighting the importance of understanding drivers across scales.

While our results might overlook specific host-flea interactions, focusing on total burden can help to elucidate the most significant drivers of abundance within and between spatial scales, and is of relevance not only for fleas, given the presence of macroparasites with varying degrees of host specificity in any given host assemblage. Indeed, generalist fleas represented 20% of the sampled fleas in our system, while fleas that not only show family-level specificity but can also parasitize sympatric rodents represented 44% of samples (Appendix S1: Table S2).

Our results highlight potential lines of research, but further exploration considering longitudinal study design, broader parasite groups, and geographic areas would be desirable to continue advancing our understanding of the effect of scale on the drivers of macroparasite abundance. Understanding these links will further expand our capacity to monitor the abundance of macroparasites and mitigate their ecosystem and health impacts.

## ACKNOWLEDGMENTS

This work was supported by a PhD scholarship from the CONACYT (Consejo Nacional de Ciencia y Tecnología) to Paulina A. Pontifes. The authors also acknowledge the National Autonomous University of Mexico (UNAM) for the support provided through a grant of the PAPIIT Program (Programa de Apoyo a Proyectos de Investigación e Innovación Tecnológica UNAM-PAPIIT), project IN225219. Paulina A. Pontifes would like to thank the Biological Sciences Graduate Program at the UNAM (Posgrado en Ciencias Biológicas) and acknowledge that this manuscript fulfills the publication requirement for obtaining a PhD degree in accordance with this program.

## CONFLICT OF INTEREST

The authors declare no conflict of interest.

## AUTHOR CONTRIBUTIONS

Adriana Fernández-González and Gerardo Suzán designed the field studies and were involved in the collection of field data. Adriana Fernández-González conducted flea identification and sample processing. Gabriel E. García Peña, Benjamin Roche, and Paulina A. Pontifes designed the statistical approach. Paulina A. Pontifes led the writing of the manuscript. Benjamin Roche and Gerardo Suzán contributed equally. All coauthors contributed comments to the manuscript.

## DATA AVAILABILITY STATEMENT

Data (Peña et al., 2021) are available from Figshare: https://doi.org/10.6084/m9.figshare.17303933.

## ORCID

Paulina A. Pontifes https://orcid.org/0000-0002-7334-8070

Gabriel E. García-Peña https://orcid.org/0000-0002-3515-6908

Benjamin Roche https://orcid.org/0000-0001-7975-4232

Gerardo Suzán https://orcid.org/0000-0003-2508-6376

## REFERENCES

Acosta, R., and J. J. Morrone. 2003. "Clave Ilustrada Para la Identificación de los Taxones Supraespecíficos de Siphonaptera de México." *Acta Zoologica Mexicana* 89: 39-53.

Anderson, S. 1972. "Mammals of Chihuahua Taxonomy and Distribution." *Bulletin of the American Museum of Natural History* 148: 149-410.

Anderson, R. M., and R. M. May. 1978. "Regulation and Stability of Host-Parasite Population Interactions: I. Regulatory Processes." *The Journal of Animal Ecology* 47: 219-47.

Barton, K. 2013. "MuMIn: Multimodel Inference." R Package Version 1.9.13. http://www.cran.r-project.org//package-MuMIn
Behrends, P., M. Daly, and M. I. Wilson. 1986. "Range Use Patterns and Spatial Relationships of Merriam's Kangaroo Rats (*Dipodomys merriami*)." *Behaviour* 96: 187-209.

Ben Ari, T., S. Neerinckx, K. L. Gage, K. Kreppel, A. Laudisoit, H. Leirs, and N. C. Stenseth. 2011. "Plague and Climate: Scales Matter." *PLoS Pathogens* 7: e1002160.

Bethony, J., S. Brooker, M. Albonico, S. M. Geiger, A. Loukas, D. Diemert, and P. J. Hotez. 2006. "Soil-Transmitted Helminth Infections: Ascariasis, Trichuriasis, and Hookworm." *The Lancet* 367: 1521-32.

Brooks, M., K. Kristensen, K. J. van Benthem, A. Magnusson, C. W. Berg, A. Nielsen, H. J. Skaug, M. Mächler, and B. M. Bolker. 2017. "glmmTMB Balances Speed and Flexibility among Packages for Zero-Inflated Generalized Linear Mixed Modeling." *The R Journal* 9: 378-400.

Burda, H., R. Šumbera, and S. Begall. 2007. "Microclimate in Burrows of Subterranean Rodents—Revisited." In *Subterranean Rodents*, edited by S. Begall, H. Burda, and C. E. Schleich, 21-33. Heidelberg, Berlin: Springer.

Cade, B. S. 2015. "Model Averaging and Muddled Multimodel Inferences." *Ecology* 96: 2370-82.

Cardon, M., G. Loot, G. Grenouillet, and S. Blanchet. 2011. "Host Characteristics and Environmental Factors Differentially Drive the Burden and Pathogenicity of an Ectoparasite: A Multilevel Causal Analysis." *Journal of Animal Ecology* 80: 657-67.

CONANP. 2013. *Programa de manejo de la Reserva de la Biosfera de Janos*. Mexico City: Comisión Nacional de Áreas Naturales Protegidas. https://simec.conanp.gob.mx/pdf_libro_27_13pm.pdf

D'Odorico, P., J. D. Fuentes, W. T. Pockman, S. L. Collins, Y. He, J. S. Medeiros, S. DeWekker, and M. E. Litvak. 2010. "Positive Feedback between Microclimate and Shrub Encroachment in the Northern Chihuahuan Desert." *Ecosphere* 1(6): 1-11.

Eads, D. A., D. E. Biggins, D. H. Long, K. L. Gage, and M. F. Antolin. 2016. "Droughts May Increase Susceptibility of Prairie Dogs to Fleas: Incongruity with Hypothesized Mechanisms of Plague Cycles in Rodents." *Journal of Mammalogy* 97: 1044-53.

Eads, D. A., D. E. Biggins, and K. L. Gage. 2020. "Ecology and Management of Plague in Diverse Communities of Rodents and Fleas." *Vector Borne and Zoonotic Diseases* 20: 888-96.

Eisen, R. J., and K. L. Gage. 2012. "Transmission of Flea-Borne Zoonotic Agents." *Annual Review of Entomology* 57: 61-82.

Gurnell, J., and J. R. Flowerdew. 2006. *Live Trapping Small Mammals: A Practical Guide*, Fourth edition, revised ed. Reading: Mammal Society.

Hartig, F. 2016. "DHARMA: Residual Diagnostics for Hierarchical (Multi-Level/Mixed) Regression Models." R Package Version 0.3.3.0. http://florianhartig.github.io/DHARMA

Himsworth, C. G., K. A. Byers, C. Fernando, L. Speerin, M. J. Lee, and J. E. Hill. 2020. "When the Sum of the Parts Tells You More Than the Whole: The Advantage of Using Metagenomics to Characterize Bartonella Spp. Infections in Norway Rats (*Rattus rattus*) and Their Fleas." *Frontiers in Veterinary Science* 7: 584724.

Hubbard, C. A. 1974. *Fleas of Western North America*. Iowa City, IA: Iowa State College Press.

Johnson, P. T. J., and J. Hoverman. 2014. "Heterogeneous Hosts: How Variation in Host Size, Behaviour and Immunity Affects Parasite Aggregation." *Journal of Animal Ecology* 83: 1103-12.

Kay, F. R., and W. G. Whitford. 1978. "The Burrow Environment of the Banner-Tailed Kangaroo Rat, *Dipodomys spectabilis*, in Southcentral New Mexico." *The American Midland Naturalist Journal* 99: 270-9.

Khokhlova, I. S. 2004. "Immune Response to Fleas in a Wild Desert Rodent: Effect of Parasite Species, Parasite Burden, Sex of Host and Host Parasitological Experience." *Journal of Experimental Biology* 207: 2725-33.

Kiffner, C., M. Stanko, S. Morand, I. S. Khokhlova, G. I. Shenbrot, A. Laudisoit, H. Leirs, H. Hawlena, and B. R. Krasnov. 2013. "Sex-Biased Parasitism Is Not Universal: Evidence from Rodent-Flea Associations from Three Biomes." *Oecologia* 173: 1009-22.

Kiffner, C., M. Stanko, S. Morand, I. S. Khokhlova, G. I. Shenbrot, A. Laudisoit, H. Leirs, H. Hawlena, and B. R. Krasnov. 2014. "Variable Effects of Host Characteristics on Species Richness of Flea Infracommunities in Rodents from Three Continents." *Parasitology Research* 113: 2777-88.

Krasnov, B. R. 2008. *Functional and Evolutionary Ecology of Fleas: A Model for Ecological Parasitology*. Cambridge: Cambridge University Press.

Krasnov, B. R., I. S. Khokhlova, L. J. Fielden, and N. V. Burdelova. 2001. "Effect of Air Temperature and Humidity on the Survival of Pre-Imaginal Stages of Two Flea Species (Siphonaptera: Pulicidae)." *Journal of Medical Entomology* 38: 629-37.

Krasnov, B. R., I. S. Khokhlova, L. J. Fielden, and N. V. Burdelova. 2002a. "The Effect of Substrate on Survival and Development of Two Species of Desert Fleas (Siphonaptera: Pulicidae)." *Parasite* 9: 135-42.

Krasnov, B. R., I. S. Khokhlova, L. J. Fielden, and N. V. Burdelova. 2002b. "Time of Survival under Starvation in Two Flea Species (Siphonaptera: Pulicidae) at Different Air Temperatures and Relative Humidities." *Journal of Vector Ecology* 27: 70-81.

Lareschi, M., and B. R. Krasnov. 2010. "Determinants of Ectoparasite Assemblage Structure on Rodent Hosts from South American Marshlands: The Effect of Host Species, Locality and Season." *Medical and Veterinary Entomology* 24: 284-92.

Linardi, P. M., and B. R. Krasnov. 2013. "Patterns of Diversity and Abundance of Fleas and Mites in the Neotropics: Host-Related, Parasite-Related and Environment-Related Factors." *Medical and Veterinary Entomology* 27: 49-58.

López-Pérez, A. M., K. Gage, A. V. Rubio, J. Montenieri, L. Orozco, and G. Suzán. 2018. "Drivers of Flea (Siphonaptera) Community Structure in Sympatric Wild Carnivores in Northwestern Mexico." *Journal of Vector Ecology* 43: 15-25.

Lüdecke, D. 2018. "ggeffects: Tidy Data Frames of Marginal Effects from Regression Models." *Journal of Open-Source Software* 3: 772.

Murray-Tortarolo, G. N., and V. J. Jaramillo. 2019. "The Impact of Extreme Weather Events on Livestock Populations: The Case of the 2011 Drought in Mexico." *Climate Change* 153: 79-89.

Nakagawa, S., P. C. D. Johnson, and H. Schielzeth. 2017. "The Coefficient of Determination R² and Intra-Class Correlation Coefficient from Generalized Linear Mixed-Effects Models Revisited and Expanded." *Journal of the Royal Society Interface* 14: 20170213.

Ostfeld, R. S., C. D. Canham, K. Oggenfuss, R. J. Winchcombe, and F. Keesing. 2006. "Climate, Deer, Rodents, and Acorns as Determinants of Variation in Lyme-Disease Risk." *PLoS Biology* 4: e145.
Peña, G. E. G., P. Pontifes, A. Fernández-González, B. Roche, and G. Suzán. 2021. "Flea Abundance in Rodents in NW Mexico." Figshare. Dataset. https://doi.org/10.6084/m9.figshare.17303933.v1.

Penczykowski, R. M., A. L. Laine, and B. Koskella. 2016. "Understanding the Ecology and Evolution of Host-Parasite Interactions across Scales." *Evolutionary Applications* 9: 37–52.

Perkins, S. E., I. M. Cattadori, V. Tagliapietra, A. P. Rizzoli, and P. J. Hudson. 2003. "Empirical Evidence for Key Hosts in Persistence of a Tick-Borne Disease." *International Journal for Parasitology* 33: 909–17.

R Core Team. 2020. *R: A Language and Environment for Statistical Computing*. Vienna: R Foundation for Statistical Computing. https://www.R-project.org/

Reid, F. A. 2006. *A Field Guide to Mammals of North America*, Fourth ed. New York: Houghton Mifflin.

Rubio, A. V., A. L. Vigueras-Galván, T. Schountz, K. Moreno-Torres, R. List, and R. E. Sarmiento-Silva. 2015. "Abundance of Hantavirus Hosts in a Landscape with Black-Tailed Prairie Dog Colonies in Northwestern Mexico." *Mammalian Biology* 80: 491–5.

Schielzeth, H. 2010. "Simple Means to Improve the Interpretability of Regression Coefficients: Interpretation of Regression Coefficients." *Methods in Ecology and Evolution* 1: 103–13.

Sheldon, B. C., and S. Verhulst. 1996. "Ecological Immunology: Costly Parasite Defenses and Trade-Offs in Evolutionary Ecology." *Trends in Ecology & Evolution* 11: 317–21.

Sikes, R. S., and W. L. Gannon. 2011. "Guidelines of the American Society of Mammalogists for the Use of Wild Mammals in Research." *Journal of Mammalogy* 92: 235–53.

Skorping, A., and K. H. Jensen. 2004. "Disease Dynamics: All Caused by Males?" *Trends in Ecology & Evolution* 19: 219–20.

Stanko, M., D. Miklisová, J. G. de Bellocq, and S. Morand. 2002. "Mammal Density and Patterns of Ectoparasite Species Richness and Abundance." *Oecologia* 131: 289–95.

Wall, R., and D. Shearer. 1997. *Veterinary Entomology*. London: Chapman & Hall.

Wilson, K., O. N. Bjørnstad, A. P. Dobson, S. Merler, G. Poglayen, A. F. Read, and A. Skorping. 2002. "Heterogeneities in Macroparasite Infections: Patterns and Processes." In *The Ecology of Wildlife Diseases*, edited by P. J. Hudson, B. T. Rizzoli, H. Grenfell, and A. P. Dobson, 6–44. Oxford: Oxford University Press.

Young, R. E., and A. D. C. Maccoll. 2016. "Spatial and Temporal Variation in Macroparasite Communities of Three-Spined Stickleback." *Parasitology* 144: 436–49.

Young, H. S., R. Dirzo, D. J. McCauley, B. Agwanda, L. Cattaneo, K. Dittmar, R. P. Eckerlin, et al. 2015. "Drivers of Intensity and Prevalence of Flea Parasitism on Small Mammals in East African Savanna Ecosystems." *Journal of Parasitology* 101: 327–35.

Zuk, M., and K. A. McKean. 1996. "Sex Differences in Parasite Infections: Patterns and Processes." *International Journal for Parasitology* 26: 1009–24.

Zuur, A. F., E. N. Ieno, and C. S. Elphick. 2010. "A Protocol for Data Exploration to Avoid Common Statistical Problems: Data Exploration." *Methods in Ecology and Evolution* 1: 3–14.

## SUPPORTING INFORMATION

Additional supporting information may be found in the online version of the article at the publisher's website.

**How to cite this article:** Pontifes, Paulina A., Adriana Fernández-González, Gabriel E. García-Peña, Benjamin Roche, and Gerardo Suzán. 2022. "Drivers of Flea Abundance in Wild Rodents across Local and Regional Scales in the Chihuahuan Desert, Northwestern Mexico." *Ecosphere* 13(4): e4013. https://doi.org/10.1002/ecs2.4013