2.4 STATISTICAL ANALYSIS
The results of PCR screening were used to classify each specimen as positive or negative to at least one infection and calculate the parasite richness in each specimen (i.e. , the number of different parasite targets detected in each sample). Parasite richness and presence/absence of infection were included as dependent variables in different generalized linear mixed models to estimate their relationship with the investigated covariates. The responses of the two species were tested separately and covariates were chosen following our ecological expectation. Specifically, in all the models we included the same covariates namely the percentage of green habitat and its fragmentation (ENN) as well as the abundance and diversity of floral resources. Moreover, each model encompassed the distance from hives and the number of hives in the surroundings. The amount of impervious cover, initially included in the models with the other covariates as a descriptor of urbanization, was excluded by the models because of the high collinearity with the percentage of green habitat evaluated through variance inflation factor (VIF) (see also the correlation plot in Supporting information, Figure S2). Presence/absence of infections were used as dependent variables in a Generalized Linear Mixed Model (GLMM) (Magnusson et al., 2017) with binomial distribution (accounting for binary presence/absence data) to evaluate changes in the probability of being infected in response to the mentioned independent variables. Changes in parasite species richness per sample in response to the considered covariates have been evaluated through GLMM with Poisson distribution (accounting to count data), in both the models sampling sites were included as a random effect. Final models were obtained through a backward stepwise model selection approach based on AIC (Zuur et al., 2009). Data analysis was performed using R (version 3.6.1; R CoreTeam 2019).