Vole outbreaks may induce a tularemia disease pit that prevents Iberian hare population recovery in NW Spain

Study site

Our study site is in an intensive agricultural landscape in NW Spain known as “Tierra de Campos”, which occupies part of three out of nine provinces of Castilla-y-León region (Palencia, Valladolid, and Zamora). This area is considered the main “hot-spot” of tularemia in Spain and Southern Europe¹⁶ and is characterized by higher-than-average vole abundances during outbreaks¹⁷.

Iberian hare abundance index

Yearly occurrence of vole outbreaks in NW Spain between 1996 and 2020 (i.e., 1997, 1998, 2007, 2008, 2014, 2019) were identified based on reports in the news (historical reconstruction¹⁸) and more recently (from 2009 onward) using common vole abundance indices obtained from live-trapping monitoring (i.e.^4,19).

To study the Iberian hare population trends we used regional hunting statistics available from the regional government (Junta de Castilla-y-León, CAZDATA Project, https://medioambiente.jcyl.es/web/es/caza-pesca/cazdata-banco-datos-actividad.html [Cited 2022 Sep 23]), which included hunting records as well as the number of hunting licences from 1974 to 2020. We used the number of hunted hares divided by the number of hunting licences each year as an abundance index for hares in “Tierra de Campos” (compiling data from the provinces of Palencia, Zamora and Valladolid). CAZDATA Project is an initiative proposed by the Hunting Federation of Castilla y León, which has the support of the regional government and, more importantly, the commitment of almost 60% of the hunting societies in the community to implement a system for monitoring hunting activity. Since this information is gathered by hunters for the benefit of the hunting activity, we are confidence on its reliability to carry out the present study.

Francisella tularensis prevalence in Iberian hares

We compiled data on F. tularensis prevalence in Iberian hares from 2007 to 2016 using previously published information from a passive surveillance program carried out by the “Regional Network of Epidemiological Surveillance” (Red de Vigilancia Epidemiológica de la Dirección General de Salud Pública) of Castilla-y-León region²⁰. This provided us with information on hare tularemia prevalence (amount of positives/number of screened individual) each year within the three provinces from “Tierra de Campos”.

Statistical analyses

To study Iberian hare population trends, we calculated an index of yearly hare population instantaneous growth rate (PGR) using the hunting bag data (hare abundance index) from 1996 to 2020. Hare PGR was calculated as follows:

where ln stands for natural logarithm, AI_t is Iberian hare abundance index on year t. and AI_t+1 is the Iberian hare abundance index on year t + 1. PGRs were estimated yearly from 1996 to 2019. This dependent variable was fitted to a Generalized Linear Mixed Model using the glmmTMB function (GLMMTMB, package glmmTMB²¹) and a gaussian family distribution and identity link function. The categorical variable vole outbreak year (i.e., with two levels: years with (1) or without vole outbreak (0), hereafter “Vole”) and “Province” (i.e., with three levels: Palencia, Valladolid and Zamora), and their interaction were used as explanatory variables. “Year” of sampling was included as a random factor (i.e., 1996–2019). Significance of the fixed effects in the models was calculated with Type II tests using the function Anova in the car package²². We previously checked the model for overdispersion and distribution fitting using function simulateResiduals (package DHARMa²³, simulations = 999). The variable PGR expresses the change between year t and t + 1. We included AI at t as a covariate in the model, in order to take into account density-dependence in hare PGR (the extent to which the abundance changes in between year t and t + 1 depends on the abundance during year t). For this to make biological sense, we rescaled the covariable AI so that it has mean equal to zero. Thus, the effect of the other predictor variables in the model (i.e., “Vole” and “Province”) was interpreted as the effect that these variables have on PGR when the abundance value is at 0. Thus, the effect of “Vole” and “Province” on PGR will be obtained by the mean value of abundance.

We assessed the effect of vole outbreak years on the Iberian hare’s population PGR by running a multiple Pearson correlation (function ggscatter) between PGR and AI, considering both, PGR for all the years of the study period (i.e., 1996–2019) and only those years where vole outbreaks were detected (i.e., 1997, 1998, 2007, 2008, 2014, 2019).

Finally, we tested for difference in the prevalence of F. tularensis on Iberian hare’s during years with or without vole outbreaks using a GLMMTMB²¹ with a binomial family distribution and a logit link function, where prevalence of F. tularensis in hares was the dependent variable, and “Vole” outbreak years and “Province” (i.e. Palencia, Valladolid and Zamora) were the responses variables. In this case, the variable “Vole” outbreak years included three levels (i.e. 0 = no vole outbreak, 1 = vole outbreak year, 2 = one year after vole outbreak), to assess if F. tularensis prevalence in hare also persist one year after a vole outbreak. “Year” of sampling was included as a random factor (i.e., 2007–2016). Due to the limited sample size, we did not include the interaction between “Vole” and “Province” to not overfit the model. We also previously checked the model for overdispersion and distribution fitting using function simulateResiduals (package DHARMa²³, simulations = 999). All analysis were carried out using the R statistical computing environment²⁴.

Source: Ecology - nature.com