Effect of host and environment-related factors on fleas of the pichi, an armadillo from Argentina

: The pichi ( Zaedyus pichiy ; Cingulata: Chlamyphoridae) is an armadillo whose ectoparasite fauna is composed of ticks and fleas. Fleas were collected from 218 pichis in southern Mendoza, Argentina, in summer and winter of 2015 and 2016. Prevalences were analyzed and differences in the intensities of the total number of fleas related to host (age, sex, weight, size and physical condition) and environment-related (seasonality and year) factors evaluated. Phthiropsylla agenoris was the only species found. Intensities of fleas were higher in 2015, in juveniles, and in males. Individuals with poor physical condition were more parasitized than those with good or normal body condition. The main explanatory variable was sampling year. This factor was directly associated with precipitation. The extreme conditions and heavy rains during the El Niño event in 2015/2016 led to environmental changes that seem to have severely affected the life cycle of fleas.


INTRODUCTION
There is ample evidence of the ecological impact of recent climate change, from polar terrestrial to tropical marine environments (Walther et al. 2002). When parasites are considered in the climate change literature, most studies focus on virulent pathogens that could become dominant in a changing climate, raising human health concerns. The majority of parasites have, however, no direct effect on human health, and the potential negative impacts of climate change on most wildlife parasites are still untested (Cizauskas et al. 2017 and references herein). A changing climate alters the availability of parasite niche space, driving a combination of habitat loss and range shifts, and potentially decreasing population growth and reproductive rates, all of which may encourage primary extinctions (Cizauskas et al. 2017). In this sense, the habitat of an ectoparasite should not be just a particular host, but a particular host in a particular habitat because of its sensitivity to factors of the off-host environment. Flea species composition in a habitat is not only determined by host species composition, but also by the environmental parameters of that habitat. These parameters determine the conditions of the burrow or the nest of the host (temperature, humidity, substrate, nest material) and thus affect flea assemblage (Krasnov 2008).
The pichi (Zaedyus pichiy; Xenarthra: Chlamyphoridae) is a small (body mass approximately 1 kg) armadillo endemic to arid and semiarid lands with firm sandy soils of Argentina and Chile (Superina 2008). Its range covers the horizontal strip from the province of La Rioja to southern Buenos Aires province in Argentina, as well as eastern Chile, south to the Straits of Magellan (Superina & Abba 2014). These semi-fossorial, solitary and predominantly diurnal mammals are opportunistic omnivores that primarily feed on insects (Superina 2008, Superina et al. 2009a. They are the only xenarthrans known to hibernate during winter and to enter torpor in warmer seasons (Superina & Boily 2007). One tick and 4 flea species have been found associated with this armadillo: Amblyomma pseudoconcolor (Ixodidae), Malacopsylla grossiventris, Phthiropsylla agenoris (both Malacopsyllidae), Tunga perforans, Hectopsylla (Hectopsylla) broscus (both Tungidae) (Mauri & Navone 1993, Superina et al. 2004, Ezquiaga et al. 2015. Future climate scenarios for the pichi's range include projected temperature increases that are, in general, more intense during the summer months. Precipitation levels are expected to decline especially in winter (Barros et al. 2015). These changes in environmental conditions could lead to a modification of species assemblages, range shifts, and local extinctions of fleas parasitizing pichis.

MATERIALS AND METHODS
The physical condition of pichis was assessed by an iterative outlier test on an ordinary-leastsquares (OLS) regression of log10-transformed weight and body length (Ezquiaga et al. 2014). Outliers in each iteration were discarded to fit the OLS regression in the next step. Outliers with positive externally studentized residuals were assigned to a good condition, whereas those with negative externally studentized residuals were assigned to a poor condition. The procedures finished when outliers were no longer detected. The remaining observations were then assigned to a normal condition. Flea prevalence, mean abundance (MA) and mean intensity (MI) were calculated as described by Bush et al. (1997).
Prevalences were analyzed by means of contingency tables using the log-likelihood ratio Chi-square statistic (Zar 2009) and applying the loglm function of MASS 7.3-45 R package (Venables & Ripley 2002). The first step was to evaluate multidimensional contingency tables for mutual independence (i.e., no interactions) among the variables year, season, age, sex and physical condition. Partial independence was tested whenever mutual independence was significant. Based on the results of the multidimensional analyses, all variables were tested by means of bidimensional contingency tables.
Generalized Linear Models (GLM) were used to evaluate differences in the intensities of the total number of fleas related to host factors and temporal dynamics. The number of fleas per individual host (count data) was the response variable, and year, age (categorical ordinal sequential), season, sex, physical condition (categorical nominal), body length and weight (continuous) were explanatory variables. A logarithmic link function was used because all predicted values had to be positive. We obtained the 127 possible combinations of 7 explanatory variables from 7 to 1. We selected the models that fit better than the null model by means of the Akaike Information Criterion (AIC), using the log-likelihood ratio test (Zuur et al. 2009).For this purpose, we calculated the weighted AIC (wAIC) (Akaike 1978, Burnham & Anderson 2002, which can be considered the conditional probability of each model of being the best (Wagenmakers & Farrell 2004). In addition, when the selected model had a wAIC lower than 0.5 or 2 or more competing models had similar wAIC values, an averaged model was computed using the MuMIn R package (Barton 2016). Additionally, we computed simple GLMs with the number of fleas per host (count data) as response variable and each of the 7 explanatory variables. To evaluate error distribution (Poisson or negative binomial), as well as overdispersion on the count or zero data, we visualized flea counts in Cleveland dot plots and histograms. Whenever graphical analyses were not conclusive, we obtained models using different error distribution (Poisson and negative binomial GLM) and accounting for zero inflated data (zero-inflated GLM) and later compared them using the likelihood ratio test (LRT) (Zuur et al. 2009). The measurement of 'false zeros' (absence of fleas recorded when fleas are present) might be random or related to sampling protocols, but there are no biological reasons to believe that the explanatory variables have an effect. A constant was therefore used in all zero-inflated GLM to model false-zero probability on the count data. All GLM and statistics analyses were performed with the lmtest 0.

RESULTS
A total of 218 pichis were captured, with no recapture (104 in 2015 and 114 in 2016). The animals comprised 104 females and 114 males, 20 of which were juveniles and 198 adults. Additionally, of the total specimens captured, 113 individuals (20 juveniles and 93 adults) were captured in summer and 105 (all adults) during winter.
A total of 2,479 fleas, all of them identified as Phthiropsylla agenoris, were collected from 212 pichis, yielding a prevalence of 97% (MA = 11.3; MI = 11.7; see Table I). Only 6 pichis, all of them adult females captured in 2016 (4 in summer and 2 in winter), were not infested with ectoparasites.
The contingency tables showed that multidimensional analyses were significant for the following parameters: mutual independence tests for the prevalence of fleas plus 5 variables (Table II) Table  SI). Consequently, the independence of the 5 variables was tested using bidimensional contingency tables, which were significant for year and sex (i.e., dependence of the prevalence on these variables) ( Table II).
The GLM analyses using a negative binomial error distribution resulted in 64 models significantly different from the null model and having a high pseudoR 2 (Table SII). Among these 64 models, the full model explained the higher proportion of total deviance (~ 44.15%) and the one with only "capture year" (a single model) explained the lower proportion of total deviance (~ 38.87%). However, both had a very low probability of being the best model (2.13% and < 1% respectively; (Table SII). The 2 best models of the 64 selected had a low probability of being the best model (11.3% and 10.3%) and include a combination of 5 and 4 explanatory variables, respectively. In addition, when the 64 models were sorted by increasing AIC, the difference between consecutive pairs than 2 (Table SII). An averaged model with the 7 explanatory variables was therefore obtained from these 64 GLM's (Table III). The averaged model indicated that the variables, except for year and sex, were weakly related to the intensity of fleas (Table III).
Simple GLM results are summarized in Table  IV. The 7 simple GLMs required using different error distributions. All were significantly different from the null model, but the p-value for season was nearly non-significant. However, only the model with year as an explanatory variable, in which intensities of fleas were higher in 2015 than in 2016, explained about 40% of deviance (Fig. 1). The remaining models explained less than 4% of deviance. Intensities of fleas were higher in juveniles than in adults, and higher in males than in females (Fig. 2). The physical condition, weight, and body GLMs needed accounting for zero inflated data and all had significant (p < 0.0001) false-zero probabilities (z-values between -8.55 and -8.61). Intensities of fleas were lower in individuals with greater body weight and length, and higher in animals with poor physical condition (Fig. 2).  (Superina et al. 2009b). Only P. agenoris was found in the present study, possibly due to an influence of the environment. This species is more prevalent than M. grossiventris in Z. pichiy from the Patagonian steppe (M.C. Ezquiaga, personal communication). Phthiropsylla agenoris could have preferences for a more xeric habitat than M. grossiventris, a species more prevalent in mesic environments, such as the Pampa. In fact, the biology of Malacopsyllidae has been poorly studied, and its life cycle is unknown (Ezquiaga & Lareschi 2012). Juveniles were more parasitized than adults, males showed higher intensities than females, and individuals in poor condition were more parasitized than those in good body condition. This was not unexpected, as parasite load varies among hosts due to their different levels of immunocompetence (Krasnov 2008 and references therein). The latter varies among individuals of a host population depending on age, sex, and nutritional and/or reproductive status. For instance, the mammalian immune system develops continuously and reaches its maximum in early adulthood, because of which juvenile individuals are more prone to contracting infections or parasite infestations due to a reduced immunocompetence (Schmid-Hempel 2011). Additionally, previous studies have reported that males are less immunocompetent than females, probably because the higher levels of androgens in males suppress immune function (Poulin 1996, Khokhlova et al. 2004). In addition, different behavioral patterns of males and females (related, among others, to territoriality, social interactions or diet) can influence exposure to infective stages of  parasites (Poulin 1996). Finally, previous studies in armadillos have shown that parasites are more abundant in hosts having poor body condition (Ezquiaga et al. 2017) and that the intensity of parasitism can be both cause and effect of its interaction with the host condition (Beldomenico & Begon 2009).
Sampling year explained about 40% of deviance, and intensities of fleas were higher in 2015 than in 2016. We hypothesize that this factor is directly associated with precipitation, as we observed a change in the environment caused by increased precipitation (e.g., flooded areas, increase of water level of lagoons, more humid soils) since our second field season. The historical mean precipitation in Malargüe is 310 mm (data from Servicio Meteorológico Nacional Argentino, 1956-2017. Previous to the first field season (February 2015), the accumulated precipitation was normal to low. Due to the El Niño Southern Oscillation (ENSO) event in 2015-2016, one of the three strongest ever recorded (Vera & Osman 2018), precipitation levels in 2015 and 2016 were considerably higher than the historical mean (518 mm and 594 mm, respectively). In Mendoza Province, the effects of El Niño are more pronounced in the summer months, which is reflected in strong summer rains and hailstorms that generate elevated ambient humidity levels (Dirección de Agricultura y Contingencias Climáticas 2016). In this context, it is reasonable to assume that the extreme environmental conditions during El Niño led to environmental changes that severely affected the life cycle of fleas. Environmental factors, such as relative humidity, may strongly affect the survival, birth and death rates of fleas and result in changes of life-history parameters of fleas, such as their abundance, reproductive rates, or patterns of parasitism (Krasnov 2008, Eads et al. 2016. Considering that the annual cycle of a particular flea species in a particular locality corresponds with seasonal climatic fluctuations (Krasnov 2008), a minimal change in precipitation levels could alter the survival rate of these fleas. In the same way, ecological studies in Argentina and other countries of South America have shown a higher prevalence of fleas in arid zones than in humid areas. For example, the prevalence of fleas was ≥ 40% in arid regions such as Ñacuñán, Mendoza Province, Argentina  Linardi et al. 1991), suggesting that humidity is the factor that most influences flea abundance. The notably higher precipitation levels at our field site during our second field season may thus have led to lower survival and/or reproductive rates of P. agenoris. The sensitivity to higher humidity levels could also be the reason why we did not find any M. grossiventris, although this flea species had been found on pichis in the same area during field studies in 2001-2006(Superina et al. 2009b. A similar pattern related to humidity was observed in prairie dogs from New Mexico (Eads et al. 2016). These authors concluded that flea densities increase during droughts, which supports our findings.
We conclude that the projected increase in frequency of extreme El Niño events (Cai et al. 2014) may have profound effects on the lifehistory parameters of fleas and even lead to their (local) extinction.