Introduction

Cyclic fluctuations of small mammals are characteristic of northern ecosystems and have been comprehensively studied in the last century (Korpimäki and Krebs 1996; Krebs and Myers 1974; Norrdahl 1995). Since the classic 10-year population cycle of lynx (Lynx canadiensis) and snowshoe hares (Lepus americanus) was described in the boreal forest of North America (Elton and Nicholson 1942), numerous studies have reported cyclic population dynamics of a diverse group of animals around the world, such as voles (Microtus spp.) and lemmings (Lemmus spp.) (Korpimäki and Krebs 1996), ptarmigan (Lagopus spp.), black grouse (Tetrao tetrix), hazel grouse (Tetrastes bonasia) (Lindén 1988; Ranta et al. 1995; Watson et al. 1998), and forest Lepidoptera (Myers 1988). These population cycles have profound effects in the functioning of the ecosystems where they occur (Ims and Fuglei 2005; Stoessel et al. 2018), influencing the dynamics and behavior of predators (Klemola et al. 1999; Korpimäki et al. 1991) and of other sympatric herbivores (i.e., the Alternative Prey Hypothesis; Angelstam et al. 1984; Breisjøberget et al. 2018).

In these systems, different predators within a community can share the same fluctuating prey. For example, in Scandinavia, vole populations follow 3–4 year cycles, which makes them a variable food source for a number of predators (Hansson and Henttonen 1985). Changes in prey availability may then influence interspecific interactions (Henden et al. 2009a; Stoessel et al. 2018). These interactions are largely based on exploitation competition, which operates when a competitor limits the availability of some shared resources (Morin 2011), and on interference competition, which involves direct negative interactions between species such as intraguild predation (Palomares and Caro 1999; Polis and Holt 1992). This puts many predators in a landscape of fear similar to that experienced by their prey (Laundré et al. 2001) and forces them to choose less optimal habitats and to alter their activity patterns in order to avoid competition and predation risk (Bischof et al. 2014; Hunter and Caro 2008; Palomares et al. 1998).

Boreal ecosystems are expected to suffer an increased vulnerability to climate change and anthropogenic disturbances (Gauthier et al. 2015; Li et al. 2018), and carnivore communities are dramatically changing under these impacts (Elmhagen et al. 2015). While specialist predators are declining, generalist predators like the red fox (Vulpes vulpes) seem to be favored by climate warming (Hersteinsson and MacDonald 1992) and by the expansion of agricultural lands and anthropogenic infrastructures (Elmhagen et al. 2015; Kurki et al. 1998). Additionally, during the last decades, rodent cycles are fading out in some parts of Europe (Hörnfeldt et al. 2005) and returning in others (Brommer et al. 2010), which might affect the interactions of sympatric carnivores dependent on rodents and change their community structure. Thus, there is a need for an improved understanding of how interspecific interactions may change, so that we can manage these species under rapidly changing conditions (Henden et al. 2009b).

In order to investigate predator species interactions in connection with prey fluctuations, we studied a predator guild in the Scandinavian boreal forests: red fox and pine marten (Martes martes), considered generalist predators (Kurki et al. 1998), and stoats (Mustela erminea) and least weasels (Mustela nivalis), considered rodent specialists (King and Moors 1979; King and Powell 2007b; Storch et al. 1990). These species rely largely on rodents as their main prey, which implies a high potential for dietary overlap and competitive interactions (Jedrzejewski et al. 1989; King and Powell 2007b).

Other bottom-up factors influencing accessibility to both main and alternative prey may also play an important role in structuring predator populations. Snow depth can limit access to voles, reducing the predator-hunting success depending on their size and hunting strategy. For instance, both red fox and pine marten are more limited by snow depth (Willebrand et al. 2017) than small mustelids (King and Powell 2007a). During periods of reduced availability of rodents, both stoats and least weasels are still able to reach rodents in their burrows, possibly contributing to the crash phase of rodent population cycles in boreal ecosystems (the Specialist Predator Hypothesis; Hanski et al. 1993; Hansson and Henttonen 1985; King and Powell 2007b; Sundell et al. 2013). On the other hand, red fox and pine marten (generalist predators) increase the proportion of alternative prey species in their diets when small rodents are few or inaccessible (the Alternative Prey Hypothesis; Breisjøberget et al. 2018; Helldin 1999; Jedrzejewski and Jedrzejewska 1992; Pulliainen and Ollinmäki 1996), which might facilitate coexistence among these competitive species (Schoener 1974). Accordingly, a range of factors, including net primary productivity and human land use, may affect the accessibility to main and alternative prey, and thereby potentially influence interactions among coexisting predators (Jahren et al. 2020).

The aim of this research was to assess the impact of rodent cyclic fluctuations in the interactions of medium- and small-sized competing predators using 8 years (2007–2014) of snow tracking data. Specifically, we assessed the effect of rodent population growth on the strength of predator’s interactions while accounting for the influence of determinant environmental variables, which can moderate the strength of these interactions (Elmhagen and Rushton 2007).

We hypothesize that the interspecific interactions among the mesopredator guild will change depending on rodent availability, following the rodent cycle. Exploitative competition is particularly important when dietary overlap is high and food availability is low. Under these circumstances, smaller predators are predicted to be more efficient at foraging than larger predators (Bagchi and Ritchie 2012). Stoats and least weasels are very well adapted to hunting small rodents, and at the same time, they can avoid predation by hiding under ground (King 1989). Therefore, in years of low rodent abundance, we expected weasels to have an advantage over red fox and pine marten (i.e., larger competitors), relaxing spatial avoidance and selecting areas preferred by rodents. On the other hand, red fox and pine marten will switch to alternative prey, thereby reducing interspecific competition (Angelstam et al. 1984; Randa et al. 2009).

During years of high rodent availability, red fox and pine marten will switch back to their main prey and interference competition may gain importance. Under strong interference competition pressure, larger competitors are usually in advantage (King 1989; Palomares and Caro 1999), and small-sized predators are more likely to use spatial avoidance as a coexistence strategy (Balme et al. 2017; Palomares et al. 2016). We therefore predict that, when rodent availability is high, small predators will avoid areas occupied by larger competitors in order to avoid aggressive interactions.

Material and methods

Study area

This study was carried out in Hedmark County (27,400 km2, 61° N 11° E), southeastern Norway. The southern part of the county is characterized by agricultural fields intermixed with large forested areas, whereas the northern part is less productive with fragmented alpine areas (Fig. 1). Forests are dominated by conifers, mainly Scots pine (Pinus sylvestris) and Norway spruce (Picea abies), intermixed with deciduous species such as birch (Betula pubescens and B. pendula), rowan (Sorbus aucuparia), aspen (Populus tremula), gray alder (Alnus incana), and willow (Salix caprea). Elevation ranges from 140 m above sea level (asl) in the south to a maximum of 2180 m asl in the north. Mean annual temperatures decrease with latitude and altitude, and therefore, the duration of snow cover changes from south to north (Pedersen et al. 2017). Predators detected were red fox, pine marten, and weasels (which included stoat and least weasel). Lager predators are also present in the study area. However, brown bear (Ursus arctos) is hibernating in winter, and densities of lynx (Lynx lynx), wolverine (Gulo gulo), and wolf (Canis lupus) are low. Consequently, we found too few tracks from these predators to include them the study.

Fig. 1
figure 1

Location of the centroids of small rodent survey areas (violet triangles) and predator snow transects (black dots) within Hedmark County, southeastern Norway

Full size image

Census data

We used 585 snow tracking transects of 2.93 km (±0.01 SE) length in average surveyed in January from 2007 to 2014 (Fig. 1) to estimate predator abundance. These transect lines were part of a Norwegian national monitoring program for large carnivores and were based on voluntary work from members of the Hedmark Chapter of the Norwegian Association of Hunters and Anglers (Tovmo and Brøseth 2014). The transect line density was 3–4 lines per 100 km2 (Tovmo and Brøseth 2014). On average, 383.5 (±17.89 SE) transects were surveyed each year under favorable snow conditions (i.e., 3.4±0.02 SE days since last snowfall). The number of wildlife tracks crossing the transect line was recorded for all species, together with snow depth and days since last snowfall. Tracks of stoat and least weasel were grouped together in order to reduce identification error given the difficulties to differentiate the tracks of both species in the field. Least weasels and stoats are referred to as “weasels” hereafter. An annual abundance index of predators was then calculated as the number of crossing tracks divided by transect kilometers and days since last snowfall.

Visual observations of rodents were obtained from line transect surveys of forest grouse species conducted in early August from 2006 to 2014 (Fig. 1). For each transect line, volunteers recorded whether small rodents had been observed. An average of six annual counts (SD=2.6) were conducted in 48 different survey areas (violet triangles on Fig. 1). The size of the survey plot averaged 56.0 km2 (SD=61.1), and in each area, an average of 15.4 (SD=10.9) transect lines (x̄=3.2 km ± 1.1 SD) with a total length of 47.1 km (SD=34.5) per survey plot were monitored (Breisjøberget et al. 2018). We calculated an index of small rodent abundance for each rodent survey area (Fig. 1) as the proportion of surveyed transects where rodents were observed (Breisjøberget et al. 2018). Since predator and rodent abundance estimates were obtained from transects that did not coincide in space, we interpolated rodent abundance index by using inverse distance weighting (IDW) in ArcGIS 10.3 (ESRI 2014) (Online Resource 1). We then extracted the interpolated rodent abundance value from each carnivore transect centroid and calculated rodent annual population growth as the difference between abundance index in year t and abundance index in year t−1 (Nt−Nt−1), giving an index of population growth that ranges from−1 to 1.

Habitat data

Snow tracking transects were related to habitat variables measured at their centroid point. These variables included elevation, latitude, relative density of human settlements, and relative density of agricultural fields, and they were obtained from digitized topographic land data from the Norwegian Mapping Authority (N250). We used the estimates for relative human settlements density and relative agricultural density calculated by Jahren et al. (2020) as follows: houses were transformed to a point layer that was later used to predict a planar kernel density map of relative settlement density. Kernel bandwidth was estimated by Gaussian approximation (Silverman 1986). Regarding the relative density of agricultural land, the geometrical center of agriculture fields was calculated and predicted planar kernel density by using agricultural-field size as z-value. Both estimated values, relative settlement density and relative agricultural density land, were then extracted to the transect centroid points.

Statistical analysis

All statistical analyses were carried out in R version 3.5.2 (www.r-project.com). We used Pearson’s correlation tests to check for correlation among continuous environmental variables, with a limit of r ≥ 0.6. Altitude and latitude were highly correlated (Pearson’s r = 0.84), and agriculture and house density were slightly correlated (Pearson’s r = 0.6). We decided to retain altitude and agriculture density for further modeling since they are important determinant factors for home range sizes and abundance of red foxes in Hedmark County (Jahren 2017; Walton et al. 2017).

Since the census data from carnivores and rodents were taken 4–5 months apart, i.e., August for rodents and January for carnivores, we used the index of rodent annual population growth from year t to test the effect on the index of carnivore abundance on year t + 1. In this way, we take into account the delay in the numerical response of predators to the abundance of their main prey (O’Mahony et al. 1999; Sundell et al. 2013). Therefore, we did not use carnivore data from year 2007 as we lacked rodent estimates to calculate rodent growth for year 2006.

We used structural equation modeling (SEM) to test how the interactions between competing predators changed in relation to the population growth of a shared prey. SEM provides a multivariate framework to develop and evaluate hypothesized networks of causal relationships, estimating the relative strength of direct and indirect paths within the system (Grace 2006; Grace et al. 2012).

We considered two alternative models, based on our predictions and on documented predator interactions in boreal ecosystems (Fig. 2). The first model represents a system where red fox interacts with weasels through pine marten (i.e., through an indirect pathway, Fig. 2a); while the second model represents a system where red fox interacts with weasels also through a direct pathway (Fig. 2b).

Fig. 2
figure 2

Path diagrams representing two alternative structural equation models (SEM). Arrows represent the direction of hypothesized causal relationships between measured variables, where predators affect each other through different direct and indirect pathways. We included the interaction of rodent growth with all the paths in each model; however, they are not represented in these diagrams for visual clarity

Full size image

To construct the SEM, we used generalized linear mixed models with negative binomial distribution to model overdispersion in a Bayesian framework, using the brm function in the brms package (Bürkner 2017) in R. We used the number of predator tracks per transect line as the response variable. Explanatory variables included the following: (i) abundance index of competing predators and (ii) environmental variables, represented by agricultural land, snow depth, and altitude. Spatial and temporal correlations were controlled using municipality and year as random effects. In order to account for track length and track accumulation over time, we included days since last snowfall and track length as offsets, thus effectively modeling the number of tracks per day per kilometer of transect (Hilbe 2014). In order to test the effect of rodent population growth in the different casual paths, we considered the interaction between rodent population growth and each of the fixed explanatory variables included in the models.

We used normal distributions with mean = 0 and standard deviation = 10 as weakly informative priors for the fixed variables and a t-Student distribution with three degrees of freedom as a non-informative prior for the random effects. We fitted the models using three chains and 3000 iterations, of which 1000 were discarded as warmup. We used leave-one-out cross-validation values using the function loo from package loo (Vehtari et al. 2019) as indicators of model fit and for model selection. We checked convergence by looking at the density distribution plots and with the Gelman and Rubin’s convergence diagnostic Ȓ (Gelman and Rubin 1992). We also calculated a Bayesian R2 (Gelman et al. 2018) for each model using the bayes_R2 function in the brms package (Bürkner 2017) to assess the variance explained by the fixed factors.

Results

In total, 15,257 red fox, 2752 pine marten, and 3164 weasel tracks were observed along 1713.78 km of transects during the 8-year survey. Two rodent peaks occurred during the study period: in 2011 and 2014 (Fig. 3). However, we only had carnivore data until the winter of 2014, so we could not test the effect of this last peak. The lowest mean rodent densities were observed in 2009 and 2012 (Fig. 3).The IDW results showed large-scale spatial synchrony of rodent cycles in Hedmark County among years (Online Resource 1).

Fig. 3
figure 3

Mean population abundance index of small rodents (percentage of tracks with rodents per survey area) ± 2SE and carnivores (tracks of red fox, pine marten, and weasels divided by days since last snowfall and track length) ± 2 SE during the whole study period (2006–2014). We aligned carnivore year with rodent year data collection (i.e., 4 to 5 months after rodent data collection)

Full size image

Between the two models tested, the model with direct interactions between all predators (Fig. 2b) had the lowest LOO value (LOOIC = 26,642.80), followed by model a (ΔLOO =  + 25.98; Fig. 2a). Therefore, we selected model b as the best fit for our data (Fig. 2b), with R2 values of 0.28 for red fox, 0.35 for pine marten, and 0.49 for weasels. For interpretation, we only considered coefficients whose 80% credible intervals (CRI) did not overlap zero, and we described separately the coefficients for the direct associations (Fig. 4a) and the coefficients for the interaction with rodent growth (Fig. 4b). At the species level, pine marten and weasels abundances were positively associated with rodent population growth (mean posterior distribution of 0.20, 80% CRI [0.09, 0.31] and 0.39, 80% CRI [0.26, 0.53], respectively; Fig. 4a, Fig. 5). Red fox, on the other hand, did not show a clear association with rodent population growth (mean posterior distribution of − 0.02, 80% CRI [− 0.06, 0.03]) (Figs. 4a and 5).

Fig. 4
figure 4

Final structural equation model (SEM) evaluating the relationship between predators, rodent growth, and environmental variables. a The coefficients for the direct associations (i.e., without the interaction with rodent growth). b The coefficients for the interaction with rodent growth (i.e., how rodent growth affects the direct paths in a). Values along arrows represent the relative magnitudes of positive (black) and negative (red) standardized path coefficients. Path coefficients whose 80% CRI overlap 0 are represented by gray arrows without values

Full size image
Fig. 5
figure 5

Posterior parameter distributions for the final structural equation model (Fig. 4) explaining direct associations and interactions of carnivores with rodent growth. The thick blue lines represent the 80% credible intervals of the posterior parameter distributions, and the thin lines represent the 95% credible intervals. Circles represent the posterior means

Full size image

Our best model showed positive direct associations among the three predators (Figs. 4a and 5). The density of agricultural lands showed a positive direct association with the three predator species, and snow depth showed a positive direct association only with pine marten. Altitude, on the other hand, showed a negative direct association with red fox and a positive direct association with pine marten (Figs. 4a and 5). Interestingly, these direct associations changed when we took into account the interaction with rodent growth (Figs. 4b and 5). The positive relationship between red fox and the other smaller predators was reduced when rodent abundance increased. Likewise, the positive association of pine marten and weasels with agriculture lands became weaker when rodents where increasing. Altitude was also positively associated with red fox when rodents were increasing (Figs. 4b and 5). For the detailed estimates of the final SEM, see Online Resource 2.

Discussion

In this study, we assessed whether interactions between sympatric predators in a boreal ecosystem are affected by the fluctuations of their shared cyclic prey. Our study shows that changes in prey resources modify how competing predators interact with each other, forcing them to adjust their competitive strategy depending on their diet specialization and competitive abilities. Pulsed resources have been shown to change the relative importance of top-down and bottom-up drivers in various ecosystems and our results reinforce the hypothesis that cyclic prey resources affect interspecific interactions within the predator guild (Greenville et al. 2014; Henden et al. 2009a; Stoessel et al. 2018).

There was a strong positive association between weasels and small rodent increase. This is in line with other studies and it highlights their degree of rodent specialization (Hanski et al. 1993; Jedrzejewski et al. 1995; Korpimäki et al. 1991). Pine marten also showed a positive association with rodent growth, but red fox was not significantly affected. Since red fox, pine marten, and weasels all rely largely on rodents as their main prey, especially in Scandinavia, they are expected to show functional and numerical responses to rodent fluctuations (Hanski et al. 1991). The red fox is considered a classic generalist predator that consumes small rodents when they are at high densities and switches to alternative prey when rodents decline (Angelstam et al. 1984; Jedrzejewski and Jedrzejewska 1992). Therefore, red fox abundance is not expected to follow rodent fluctuations as closely as the specialist predators, and alternative resources such as woodland grouse (subfamily Tetraoninae), mountain hares (Lepus timidus), ungulate carcasses, and anthropogenic food play an important role in stabilizing red fox population at low rodent densities (Carricondo-Sanchez et al. 2016; Killengreen et al. 2011; Selås and Vik 2006; Willebrand et al. 2017). Weasels are strict rodent specialists that prey on small rodents regardless of their availability and, thus, their survival and reproduction is strongly influenced by rodent fluctuations (Jedrzejewski et al. 1995; Sundell et al. 2013), as confirmed in our study. Pine martens, on the other hand, can be considered intermediate, because they may specialize on a few prey items such as small rodents and squirrels, but they can also switch to alternative food resources if needed (e.g., eggs, small birds, and carcasses; Jedrzejewski et al. 1993; Pulliainen and Ollinmäki 1996). Therefore, their abundance may correlate with rodent cycles, but not to the same extent as weasels.

Regarding species interactions, the three predators were positively associated with each other, which is likely to be a result of similarities in habitat and diet requirements. Positive associations between negatively interacting species are probably facilitated by landscape heterogeneity and differences in activity patterns; i.e., species may coexist by using different microhabitats and temporal niches (Bischof et al. 2014; Lesmeister et al. 2015; Viota 2012). However, this positive association between predators was strongly influenced by rodent growth: pine marten and weasels were less likely to be associated with red fox when rodents were increasing. As a generalist predator, the red fox tends to prey mostly on rodents during the increasing phase of the cycle. In addition, Dell’Arte et al. (2007) found that red foxes increased their predation on small mustelids when vole densities were high. Consequently, subdominant specialist predators like weasels are expected to experience higher impacts of interference competition when resource densities are high (Dell’Arte et al. 2007; Henden et al. 2009a). Our results suggest that under such conditions, avoidance of larger and dominant predators may drive smaller competitors to use less optimal habitats. This would explain the weaker association that we observed between red fox and mustelids during the increasing phase of the rodent cycle.

Conversely, during the decrease phase of the rodent cycle, pine marten and weasels showed a stronger spatial association with red fox. Under a low prey resource scenario, smaller predators are usually superior in exploitative competition because of their enhanced ability to hunt specific prey (King 1989). Furthermore, pine marten and weasels are considered very efficient predators, even under an increased risk of competitive interference, because of their ability to access food in a diversity of habitat strata (Bischof et al. 2014; Brainerd and Rolstad 2002). Hence, they can afford to select areas preferred by rodents during low phases despite an increased spatial overlap and potential competition with red fox. This has been shown for other pairs of sympatric carnivores such us red fox and arctic fox (Vulpes lagopus), which were positively associated only during the low phase of the rodent cycle (Stoessel et al. 2018). Moreover, under a scenario of resource constraints and stronger spatial associations among competing predators, it is likely that habitat heterogeneity and temporal niche segregation play an even more important role to facilitate species coexistence (Bischof et al. 2014; de Satge et al. 2017; Linnell and Strand 2000; Sergio et al. 2003).

All predators were positively associated with agriculture land, which represents productive areas where small rodents are abundant and easily accessible (Panzacchi et al. 2010). These areas are also associated with human settlements and anthropogenic food, which is an important alternative food source for red foxes (Killengreen et al. 2011; Rosalino et al. 2010). Therefore, productive agricultural lands might offer a continuous source of food for red foxes throughout the year, with anthropogenic food being an alternative resource when rodent numbers decrease. This could explain why the positive association between red fox and agricultural lands is not affected by changes in rodent abundance. Pine marten and weasels, however, were more strongly associated with agricultural land during the decrease phase of the rodent cycle. As rodent specialists, weasels tend to concentrate their foraging activity in more productive areas such as agricultural lands when food availability is low (Klemola et al. 1999), thus increasing their spatial association with red fox during low rodent phases.

Altitude had a negative effect on red fox abundance, yet this negative association was influenced by rodent growth. During the increase phase of the rodent cycle, red fox showed a positive association with altitude. Andreassen et al. (2019) found that vole amplitudes were higher at higher elevations, with greater maximum vole densities at high elevations during peak years. This could explain why red foxes move to higher elevations during the increase phase of the rodent cycle, since rodents seem to be more abundant at higher altitudes during this phase of the cycle. Furthermore, red foxes are expected to switch to alternative food sources when rodent abundance decreases, and this may induce a shift in habitat use towards lower altitudes close to human activities (Šálek et al. 2015).

Pine marten was the only predator to show a positive association with snow depth, and this association was even stronger during the increase phase of the rodent cycle. Pine marten seems to be more efficient in killing prey under the snow than red fox (Willebrand et al. 2017); therefore, they might use deep snow areas as a way of reducing interference competition when the overlap with red fox is higher. Alternatively, snow depth may be correlating with some other variable that we did not measure, with positive effects on pine marten abundance. Although some studies indicate that snow depth reduces red fox hunting success and probability to survive (Jedrzejewski and Jedrzejewska 1992; Lindström and Hörnfeldt 1994; Selås and Vik 2006; Willebrand et al. 2017), we did not find any particular association between red fox and snow depth. However, we did not measure snow compaction, which may affect predators’ habitat use (Pozzanghera et al. 2016).

Predator interactions tend to be very complex, as they are influenced by several intrinsic and extrinsic factors. Moreover, complexity increases when we consider the effects of a cyclic shared prey and strong seasonality. Here, we show that different phases of the prey cycle drive changes in the interactions among sympatric carnivores.

The results of this study are restricted to the interactions of a specific guild of predators in a boreal ecosystem. However, we believe that extrapolation on the functioning of other communities in which several predators share a cyclic prey can be made with caution. Under a climate change scenario, small mammal cyclic fluctuations are expected to change, for example, due to changes in seasonality. Given our results, we could speculate that these changes might alter the dynamics of the predator community by changing top-down and bottom-up interactions. This study highlights the importance of taking into account prey dynamics when studying interspecific interactions among predators. It also shows the need for monitoring in order to anticipate undesirable outcomes such as enhanced generalist predator abundances to the detriment of specialists predators.