Spatial relationships at a coarse (monthly) temporal scale

For the monthly scale relationships, we used generalized linear mixed models (GLMMs) with negative binomial errors⁴⁴ to evaluate the relationships between the spatial patterns of red fox, our focal species, and those of the wolf. We set our models in the R package ‘glmmTMB’⁴⁵ and ‘lme4’⁴⁶. We built a full model for each study area. The number of detections of red fox in each camera-trapping location for each deployment was fitted as the response variable. In all models, the log (number of camera operating days) was included as an offset variable to standardize the response variables according to the actual sampling effort. We included the following predictors: (1) the detection rates of wolf; (2) people detection rate; (3) habitat (Maremma RP: oakwood; pinewood; shrub wood; ecotone/meadows; Białowieża PF: broad-leaved forest; coniferous forest; mixed forest; transitional woodland/shrub); (4) Distance from the nearest human settlement. We initially assessed whether to include the lynx among predictors for Białowieża PF models, considering its potential role as a competitor for the fox^26,47. However, the lynx has a low density in Białowieża PF³⁷, which is supported by the extremely small number of detections collected throughout our study (only 43 detections in 4 years). Hence, we decided to use only the wolf even to a better comparison between the areas’ models. In this study we put as predictor the variable ‘habitat’ on the camera location to control for potential effects of habitat type on the detection rate of our focal species, that can change in relation to different habitat. The habitat in each camera-trapping site was assessed through local land use maps (Maremma RP: see⁴⁸; Białowieża PF: clc.gios.gov.pl). The distance from the nearest human settlement was estimated in QGIS, by using the layer of land use and cover 2018 Corine Land Cover of Italy and Poland in which we included as ‘permanent human settlements: (a) continuous and discontinuous urban structures (with at least ten houses)²⁷; (b) industrial and commercial units. Data on people detection rates were not available for 2018 in Maremma RP, because people detections were not reported during this period. Thus, for Maremma RP, we conducted two analyses: a first one, considering only the years when people data were available and including all the above-mentioned predictors, and a second one, considering all the four study years, and where we did not include people detection rate among predictors. Since the results pointed in the same direction with both two approaches, in the main text we show only the results of models built through the datasets without the year 2018 that includes the people detections (for models with full dataset, see Supplementary Materials 1S–2S). We used the ‘Year’ as the random effect to account for potential differences between years. Absence of collinearity among linear predictors was checked preliminarily: pairs of predictors included in models did not show correlation coefficients >|0.6|⁴⁹.

Subsequently, all possible models were calculated for each study area, including all possible combinations of considered predictors, since all of them reflected different a priori hypotheses, and were evaluated through a model selection procedure based on a comparison of AICc scores (Akaike Information Criterion). We used the model selection with the nesting rule to avoid retaining overly complex models^50,51. We identified the best model as the most parsimonious one, i.e., the one having the lowest AICc^50,52. Moreover, we selected for inference all models with AICc ≤ 2^51,52, and among these, those which were not more complex versions of any simpler model⁵². Model selection was conducted using the R package ‘MuMIn’⁵³. We estimated the parameters (95% confidence intervals and B coefficients, which is the degree of change in the response variable for every 1-unit of change in the predictor variable) of the best model by using the R package ‘lme4’⁴⁶. Then, the best model was validated by visual inspection of the distribution of residuals⁴⁴ through the ‘DHARMa’ package⁵⁴. Model weight was standardized within the subset of selected models.