Fear or food – abundance of red fox in relation to occurrence of lynx and wolf

Apex predators may affect mesopredators through intraguild predation and/or supply of carrion from their prey, causing a trade-off between avoidance and attractiveness. We used wildlife triangle snow-tracking data to investigate the abundance of red fox (Vulpes vulpes) in relation to lynx (Lynx lynx) and wolf (Canis lupus) occurrence as well as land composition and vole (Microtus spp.) density. Data from the Swedish wolf-monitoring system and VHF/GPS-collared wolves were used to study the effect of wolf pack size and time since wolf territory establishment on fox abundance. Bottom-up processes were more influential than top-down effects as the proportion of arable land was the key indicator of fox abundance at the landscape level. At this spatial scale, there was no effect of wolf abundance on fox abundance, whereas lynx abundance had a positive effect. In contrast, at the wolf territory level there was a negative effect of wolves on fox abundance when including detailed information of pack size and time since territory establishment, whereas there was no effect of lynx abundance. This study shows that different apex predator species may affect mesopredator abundance in different ways and that the results may be dependent on the spatiotemporal scale and resolution of the data.


Results
Landscape Level. Fox abundance increased with increasing proportion of arable land, and this variable was the most informative parameter explaining the variation in fox abundance at the landscape level (RVI = 1.0, Fig. 1a and Tables 1 and 2). Fox abundance at the landscape level was also positively related to lynx abundance (RVI = 0.78, Fig. 1b and Tables 1 and 2). There was no effect of wolf occurrence (independent of method used to classify wolf occurrence) on fox abundance ( Fig. 2a and Table 1) at the landscape level, nor was there any support for models including the vole density index (Table 1).

Figure 1.
Abundance of red fox (crossings 10 km −1 24 h −1 ) in relation to (a) arable land (%) and (b) lynx abundance (crossings 10 km −1 24 h −1 ). Model averaged predictions derived from the highest ranked models (Table 1) are shown (solid lines) with associated 95% confidence intervals (grey shading). Predictions for each explanatory variable were made while holding the other variables constant at mean values (arable land 8% and lynx abundance 1.5 crossings 10 km −1 24 h −1 ). Red fox abundance predictions were back-transformed to their normal scale for the figure. For model averaged parameter estimates see Table 2. Red fox and lynx abundance were surveyed by snow-tracking in wildlife triangles in south-central Sweden from [2001][2002][2003] Table 1. Highest ranked candidate models relating red fox abundance (crossings 10 km −1 24 h −1 ) to arable land (%), vole density index (low or high), lynx abundance (crossings 10 km −1 24 h −1 ), and wolf occurrence. We used four separate measurements of wolf occurrence: 1a) wolf abundance (crossings 10 km −1 24 h −1 ); 1b) a four-level wolf category (I: inside an observed wolf territory, II: inside an average wolf territory, III: inside a maximum wolf territory, or IV: outside a wolf territory); 2) wolf pack size (range 2-8 individuals), using only triangles inside wolf category I; and 3) time since territory establishment (year of territory establishment, 1-2 years after establishment, or ≥3 years after establishment), using only triangles inside wolf category I. Wolf Territory Level. At the wolf territory level there was a weak negative relationship between fox abundance and wolf pack size (RVI = 0.29) in combination with the proportion of arable land, although arable land and vole density index were the best predictors of fox abundance (RVI = 0.56 and 0.53, respectively, Fig. 3 and Tables 1 and 2). Including time since territory establishment in the analyses showed that there was a  Table 2. Model-averaged parameter estimates with standard error (SE) for each variable retained in the best models (ΔAIC c ≤ 2) in Table 1. We used four separate measurement of wolf occurrence, as defined in Table 1. Vole density index is a two-category variable where the parameter estimate is the difference in red fox abundance for low compared to high (intercept) vole density index. Time since territory establishment class is categorical, where the parameter estimate is the difference in red fox abundance at 1-2 years or ≥3 years after territory establishment compared to the year of establishment (intercept).

Figure 2.
Upper panel (a) shows mean (±SE) red fox abundance (crossings 10 km −1 24 h −1 ) in relation to a 4-level wolf category: (I) inside an observed wolf territory, (II) inside an average wolf territory, (III) inside a maximum wolf territory, and (IV) outside a wolf territory. Lower panel (b) shows mean (±SE) red fox abundance from wildlife triangles located inside estimated wolf territory borders (i.e., wolf category I) in relation to time since territory establishment (years) where year 0 is the year of territory establishment by a wolf pair. Red fox abundance was surveyed by snow-tracking in wildlife triangles in south-central Sweden from 2001-2003. negative temporal effect of wolves on fox abundance (RVI: time since territory establishment = 0.57, vole density index = 0.49, arable land = 0.47). Fox abundance decreased 1-2 years after wolf territory establishment and after these first two years fox abundance rebounded again, although to a slightly lower level than at the year of wolf territory establishment (Fig. 2b). In contrast to the landscape level, there was no relationship between fox abundance and lynx abundance at the wolf territory level (Tables 1 and 2).

Discussion
Although both wolves and lynx had a top-down effect on fox abundance in our study area during winter, bottom-up processes were more influential, as the proportion of arable land was the key indicator of fox abundance. The almost complete loss of apex predators during the 19 th century, in combination with geographical variation in ecosystem productivity, has been shown to be a key factor in the fox population increase in Scandinavia 8 . The agricultural landscape holds the highest densities of foxes 28 . Also, lynx in Scandinavia are attracted to agricultural fields because their main prey, roe deer, are clumped around this habitat [39][40][41] . Wolf territories, on the other hand, are composed by a small amount of open areas, such as agricultural fields, suggesting an adaptation to avoid humans 42 . Thus, in our study area, fox and lynx likely spatially overlap more than fox and wolves. The avoidance of open areas by wolves is further supported by the fact that the proportion of arable land in wildlife triangles inside wolf territories averaged 2% (range 0-20%), compared to 8% (range 0-76%) in the whole study area. The importance of vole density on fox abundance may be more pronounced in less favorable habitats 29 . Thus, the lower proportion of arable land within wolf territories may explain why the vole density index affected fox abundance at the wolf territory level, but not at the landscape level. Consequently, the lowest fox abundance was inside wolf territories that had a low proportion of arable land and low vole density index. This supports the conclusion that bottom-up processes were an important factor affecting fox abundance. At the landscape level, there was a positive relationship between fox and lynx abundance. Simulations where we added observation uncertainty to the lynx abundance measurements (i.e., coefficient of variation ranging from 0.2 to 3) revealed that our results were robust; it was a high probability (i.e., ˃97%) of a positive relationship between fox abundance and lynx abundance (see Supplementary Methods for further details). This result was surprising since several studies have shown a negative relationship between these two species 8,9,28,43 . For example, Elmhagen & Rushton 8 used historic data to demonstrate a mesopredator release of fox after lynx and wolves were exterminated from Sweden, indicating that foxes had been top-down limited before the extermination of these two apex predators, and Elmhagen et al. 9 found that this process was reversed when lynx recolonized parts of Finland in recent time. Our data did not allow us to determine whether the positive effect of lynx abundance on fox abundance in our study was a demographic or a behavioral response. Although previous studies have suggested a negative demographic effect of lynx on fox populations 8,9,13 , we suggest that the effect we observed was a behavioral response by fox to lynx abundance. Foxes are attracted to sites with olfactory cues from lynx, a behavior that should increase the possibility to access carrion from lynx killed prey 44 . However, from lynx scent marks, a fox might also gain information about the existence of lynx in its home range, a knowledge that may be used to adjust behavior in order to reduce the risk of encountering lynx 44 . Studies of multiple predator sympatry in African and Asian ecosystems have shown fine-scale behavioral tactics and trade-offs associated with resource availability [45][46][47] . Hence, we suggest that the positive relation between fox and lynx abundance indicates that foxes are attracted to lynx to increase their chances to find carrion, possibly facilitated by using a fine-scale spatiotemporal segregation strategy to minimize the risk of being killed. Two alternative hypotheses are that 1) both lynx and fox are attracted to a third, currently unknown, factor, or 2) species within the same mammalian family are likely to compete more than species in different families 48 , which causes fox to avoid wolves and not lynx. In any case, our findings contradict the hypothesis that the avoidance of apex predators could suppress mesopredator abundance more than predator-induced mortality 14 , at least with respect to fox and lynx in our study area.
Whereas we found no effect of wolves on fox abundance at the landscape level, we show a negative temporal effect of wolf occurrence at the wolf territory level. This effect was most pronounced the first two years immediately after wolf territory establishment. We suggest that this negative effect is explained by a combination of the novelty of wolf occurrence and their local abundance, i.e., pack size. Zimmermann et al. 38 showed that large wolf packs (6-9 wolves) in Scandinavia had a kill rate that provided just enough, or even less, prey biomass than required, while small packs (2-3 wolves) had a kill rate that provided more prey biomass than required, thereby leaving a larger surplus of carrion. Furthermore, the probability that foxes encountered wolves increased in large pack territories due to the tendency of large wolf packs to split up in sub-groups during winter 49 . The establishment of a wolf territory usually starts with a breeding pair, followed by an increased number of wolves the consecutive winter after the pair's first reproduction. Fox abundance one and two years after wolf establishment was reduced by 50% compared to the establishment year and with triangles outside wolf territories (wolf category III and IV in Fig. 2a). The reason that we did not see any negative or positive effect during the establishment year could be either that the pack consisted only of the establishing pair, i.e., the smallest possible pack size, and/ or that the pair was formed late in the winter. Consequently, we suggest that the pronounced negative effect on fox abundance the years immediately following wolf establishment was due to an avoidance reaction by foxes to the novelty of wolf presence which faded over time, or the increased mortality of naïve foxes. We suggest that a demographic effect is less likely than a behavioral effect as wolf scat analyses indicate that intra-guild predation by wolves on foxes is low in Scandinavia (<0.2% of wolf scats (n = 2091) contained hair from fox) 50,51 . However, if wolves kill foxes primarily to reduce competition (i.e., foxes are not consumed), wolf scat analysis might underestimate intra-guild predation. But studies of wolf predation in Scandinavia show that few wolf-killed foxes are found nearby the remains of wolf-killed ungulates 52, 53 which also indicates low intra-guild predation.
We argue that food competition is not a plausible cause for the negative effect of wolves on fox abundance. Wolf occurrence results in increased carrion availability, a food source frequently utilized by foxes in our study area 34 . Although wolves occasionally kill smaller prey species (hares (Lepus timidus, Lepus europeus) and forest grouse (Tetrao tetrix, Tetrao urogallus)) 53 , that are also prey for foxes, wolf predation on these species is likely not high enough to affect their abundance. After three or more years since the wolf territory was established, fox abundance rebounded. Unfortunately, our data set did not allow us to assess whether or when fox abundance rebounded starting at pre-wolf numbers. We were unable to conduct a similar analysis of possible time effects on the response of foxes to lynx, since lynx had already established in the area several years before our study started 13 . It is possible that the fox population underwent a learning process, comparable to what appeared to have occurred with wolves, which enabled them to co-exist with lynx. This process is likely more advanced with the lynx population when compared to the wolf populations, due to the longer time of co-existence between fox and lynx in our study area.
In our study, foxes were positively affected by lynx abundance, while wolf occurrence had a negative temporal effect on fox abundance, although support for the latter effect was weak. The observed response of foxes to these two apex predator species contradicts the results from Pasanen-Mortensen et al. 43 and Pasanen-Mortensen & Elmhagen 28 , who found a strong negative effect of lynx, and no effect of wolves, on fox abundance at the Eurasian scale. This is even more puzzling considering wolves, which prey mainly on moose, likely leave more carrion than lynx, which mainly rely on smaller roe deer 34,41 . Also, because the risk of intra-guild predation is greater when two species are of similar size, there should be greater risk for foxes in lynx areas than in wolf areas 48 . However, the extent of top-down control may depend on the density of the apex predators 11,54 . Unfortunately, a comparison between apex predator density in our small-scaled, short-term study (stable or slightly decreasing lynx population and a recolonizing wolf population) and previous studies conducted in Scandinavia is not straightforward due to different measurements of apex predators occurrence (presence/absence 28,43 , harvest statistics 8 ), and variation in spatial scales (country 8,9 , continent 28, 43 ), time span (long-term 9 ), and population phase (recolonized lynx population 9 ). Our study suggests that the relationship between fox, lynx and wolves might be more complex than observed at a larger scale, or by only using presence/absence data. Actually, Guillaumet et al. 55 found that by reducing the spatial scale of analysis there was a positive synchrony between Canadian lynx and coyote (Canis latrans) numbers, contrary to the negative pattern observed at a larger spatial scale. Furthermore, temporal partitioning of activity may facilitate the coexistence of foxes and apex predators 56 . The importance of spatiotemporal scale is further illustrated by our results that at a larger spatial scale (the landscape level) fox abundance and wolf occurrence were not correlated (c.f., Pasanen-Mortensen et al. 43 ), whereas at a smaller spatial scale (the wolf territory level) we found a negative temporal effect of wolf occurrence on fox abundance. Our results emphasize the importance of considering multiple spatiotemporal scales when investigating species interactions. Further studies at an even smaller scale than used in our study (e.g., locations of tracks in space and time within wildlife triangles or using GPS-collared individuals) are warranted to gain further insight into the impact of lynx and wolf on fox, and determine how this is influenced by habitat.

Materials and Methods
Study System. The study was conducted from 2001-2003 in the south-central part of Sweden within the counties of Värmland and Örebro (referred to as the landscape level, 58°50′-60°N, 12°-17°E, approximately 26 000 km 2 ). The area is dominated by boreal forest interspersed with numerous lakes, rivers and mires that contribute to a more variable landscape than is found in other areas within the boreal forest zone. The proportion of arable land increases from north-northeast to the south-southwest. About 74% of the total land area of Värmland, and 67% of Örebro, is forested and the corresponding area of arable land is 6% and 16% 57 . The number of days with snow ranges from approximately 75/year in the south to approximately 175/year in the north 58 . For a detailed description of the study area see Wallgren et al. 59 .
In the 1980s, the density of the fox population decreased dramatically due to an epizootic outbreak of sarcoptic mange 60 , and subsequently recovered in the 1990s 61 . Lynx reoccupied south-central Sweden in the beginning of the 1950s 26 and, by 2002-2003, the total Swedish lynx population was estimated to 1300-1600 individuals 62 . During the study period the number of lynx family groups (i.e., female with kittens, the monitoring unit within the national lynx monitoring program) 62 in the study area was 69 (2001) and 49-55 (2003). Roe deer are the main prey for lynx within the study area 41 . Wolves returned to the study area in the early 1980s through natural re-colonization and, in 2002-2003, the total Swedish wolf population was estimated at 84-100 individuals 63 Fig. 4). For each triangle visit, we calculated the relative abundance of each species as crossings 10 km −1 24 h −1 , following Kurki et al. 37 . Kurki et al. 37 demonstrated that abundance based on snow-tracking was significantly correlated with hunting bags of two mesopredators (fox and pine marten (Martes martes)) although other factors also contribute to variation in bag. The use of snow-tracking as an index of abundance is further exemplified by Lyly et al. 67 , but may not work for all species or areas 68 .
As available data from snow-tracking and locations from VHF-collared wolves (two locations per week on average) did not give a full representation of wolf territory borders, we measured distances between the centre of each wildlife triangle to the centre of the nearest wolf territory (using ArcView 3.0, ESRI Corporation, Redlands, CA) and classified wolf occurrence for each wildlife triangle according to the following four categories: (I) inside an observed wolf territory if the triangle centre was located within the confirmed borders of a wolf territory, (II) inside an average wolf territory if the triangle centre was not located inside an observed wolf territory but within a distance corresponding to the average wolf territory radius from the nearest territory centre (18.0 km, average wolf territory size = 1017 km 2 ) 70 , (III) inside a maximum wolf territory if the triangle centre was not located inside an observed or an average wolf territory but within a distance corresponding to the maximum wolf territory radius from the nearest territory centre (23.0 km, maximum wolf territory size = 1676 km 2 ) 70 , or (IV) outside a wolf territory if the distance from nearest territory centre was longer than the radius used in category III. Thus, a triangle belongs to only one category.
Habitat Composition and Vole Density Index. As arable and pasture land are known to have a positive impact on fox density due to food abundance 8, 37 , we calculated the percentage of these habitats (pooled and referred to as arable land), using the Swedish Land Cover Data 71 (similar to the European Union classification system Corine Land Cover (CLC)). This was done within a 1 km buffer zone around each wildlife triangle (mean = 8%, range: 0-76%) in order to account for the effect of land composition on fox abundance (Fig. 4). For a detailed description of the methodology and classification see Wallgren et al. 59 . Furthermore, the fox population tends to fluctuate following the variation in vole density with a time lag of one year 35,72

Statistical Analyses. We conducted all statistical analyses in R version 3.2.2 74 using General Linear Mixed
Models in the lme4 package 75 with fox abundance (crossings 10 km −1 24 h −1 ) as a response variable. Since fox abundance is a continuous variable, and its variance increased with increasing mean (i.e., similar to a log-normal distribution where σ = μ 2 ), we used log-transformed fox abundance and a Gaussian error distribution in all models to fulfil the assumptions of linearity. We used wildlife triangle identity and year as random effects in all models to account for repeated measurements and year effects. To test the effect of apex predators on fox abundance at the landscape level, we included four explanatory variables: wolves (classified according to wolf abundance or wolf category, as described in (a) and (b) above, exclusive of each other), lynx abundance (crossings 10 km −1 24 h −1 ), proportion of arable land (%), and vole density index (2-level categorical: low or high). We also included pairwise interactions between habitat and all other explanatory variables as well as an interaction between wolf and lynx abundance.
To assess the effect of wolf pack size and time since territory establishment on fox abundance, we restricted our analyses to data from wildlife triangles that were exclusively inside the observed borders of a confirmed wolf territory (wolf category I, n = 77, referred to as the wolf territory level). We used wolf pack size and time since territory establishment (3-level categorical: year of establishment (referred to as 0), 1-2 years after establishment, or ≥3 years after establishment) together with lynx abundance, proportion of arable land, and vole density index as explanatory variables in the models.
For each analysis, we compared candidate models using the sample-size corrected Akaike Information Criterion (AIC c ) and AIC weights (w i ) from the 'MuMIn' package 76 in R. Models with ΔAIC c <2 were used to generate model-averaged parameter estimates 77 . We used a bootstrap method implemented in R using the 'ez' package 78 to calculate predictions and 95% confidence intervals for mixed models. Additionally, we used AIC c weights on the full candidate model set to generate Relative Variable Importance weights (RVI) for each explanatory variable. Data Availability. All data analysed during this study are included in the Supplementary Data file.