Top–down influences of upper trophic levels on lower trophic levels are a widely studied phenomenon in ecology (Terborgh and Estes 2010). These cascading effects often are among the most important factors in structuring many ecological communities (Ripple et al. 2014). Typically, these effects are identified in systems that contain at least three distinct trophic levels, with species occupying the highest level providing positive, indirect effects to those species occupying the lowest trophic level through limitations on the abundance of organisms in the intermediate trophic level. Among the most widely recognized examples are the orca–sea otter–kelp system of the Pacific coast of North America (Estes et al. 1998) and the wolf–elk–aspen system of Yellowstone National Park (Ripple et al. 2001). A variation of this concept is the mesopredator release hypothesis, which posits that a similar three-level hierarchy can occur within a single trophic level of carnivores provided that there is still size structuring of the species. Under mesopredator release, the abundance of the smallest member of the carnivore guild is positively associated with the largest member of the guild through mediation of the intermediate member's population size (Prugh et al. 2009). Thus, cascading effects of large carnivores are not limited strictly to trophically structured systems. This effect has been seen in numerous mesocarnivore systems across the globe (Elmhagen et al. 2010).

Levi and Wilmers (2012) provided evidence of cascading effects of recovering wolf populations within a single trophic level. They found evidence that recovering wolf Canis lupus populations in Minnesota suppressed coyotes Canis latrans, an intermediate predator, which subsequently benefitted smaller canids in the form of red foxes Vulpes vulpes. Interference competition is a known behavioral mechanism occurring between wolves and coyotes, particularly in areas with recovering wolf populations (Berger and Gese 2007, Merkle et al. 2009). Similarly, direct mortality due to wolves is commonly observed in coyote population where the two species are sympatric (Ballard et al. 2003). Thus, a positive association between wolf and fox populations may be expected through competition-mediated coyote populations (Newsome and Ripple 2015). The effects of wolves on coyote populations can also benefit other non-canid species that are negatively affected by coyotes through predation (Berger and Conner 2008). At local scales, wolves have been shown to alter the behavior of coyotes in the Great Lakes region (Flagel et al. 2017). However, it is not known if these localized effects translate into population-level impacts among these species, or if these effects impact populations of other canids.

In Wisconsin, wolves were extirpated by 1960 due to a long history of persecution that included government sanctioned bounties (Wydeven et al. 2009). Wolves began recolonizing Wisconsin in the mid-1970s, likely dispersing from Minnesota, and exhibited an initial period of limited growth characterized by an Allee effect (Stenglein and Van Deelen 2016) followed by an extended period of rapid population growth (Fig. 1). Harvest seasons in 2012–2014 limited wolf population growth until re-listing of the Great Lakes wolf population occurred in December 2014, and the population has continued to expand since that time (Wiedenhoeft et al. 2018). Wolves in Minnesota were never extirpated and currently occur at much higher densities than in Wisconsin. The differing recovery and management histories of wolves in Minnesota and Wisconsin may limit the transferability of the results of Levi and Wilmers (2012), potentially limiting our understanding of the ecological impacts of wolf recovery in Wisconsin. Thus, we sought to determine the relative impacts of recovering wolf populations on coyote and fox population growth in Wisconsin and how these impacts may influence mesocarnivore communities.

Methods

Snow track surveys

We used data from the Wisconsin Department of Natural Resources (WDNR) winter furbearer survey program. This is a long-term monitoring program initiated in 1977 to monitor distributions and relative abundances of furbearers (species harvested for their pelts) and furbearer prey in 20 counties in northern Wisconsin. Although initially targeting furbearer populations in the northern third of Wisconsin only, the program expanded to include a southern zone comprising nine counties and to include wolves beginning in winter 1998–1999. Within each county, two 16 km surveys were conducted annually along two pre-designated routes separated by ≥16 km (Fig. 2). Logistical constraints often limited the number of surveys conducted, with the total number of surveys conducted within a year ranging from 8–37 in the northern region and 3–10 in the southern region. These routes occurred along remote roads with minimal traffic and snow plowing. Generally, surveys were completed between November and December when roads were still passable and snow had not accumulated to levels that restrict furbearer movement (>30 cm). Some surveys were conducted following 1 January if survey conditions remained poor during the standard survey months. Surveys were conducted during daylight hours after a snowfall which ended prior to 18: 00 the previous day to ensure fresh tracking substrate. Surveys were not conducted at temperatures ≤–17°C. Additionally, surveys were not conducted during the state's firearm harvest season for white-tailed deer Odocoileus virginianus or during the ten days following this harvest season to avoid consequences of increased human activity (opening day of Wisconsin's firearm deer season is the Saturday before the Thanksgiving holiday in mid- to late-November). Surveys were conducted by at least two individuals traveling in a vehicle at speeds slow enough to allow observation of all tracks (12–16 km h–¹). After a track was located, surveyors followed the track on foot for accurate identification and to determine if the tracks belonged to one or several animals. To the best of the observer's ability, each observation denoted the tracks of an individual animal; multiple observations of the same animal's tracks were not recorded. Generally, tracks were assumed to be of one animal if the species-specific tracks were observed <0.5 km apart. Observers record the number of individual animals detected per 0.8-km segment along the survey transect, resulting in 20 segments per transect. All surveys were conducted by trained WDNR personnel.

Figure 1.

Annual mid-winter minimum wolf population in Wisconsin, 1980–2017.

Figure 2.

Study regions and distribution of track-survey transects within Wisconsin.

Currently WDNR uses a naïve occupancy estimate derived from the snow-track surveys to monitor relative abundance of these species. Relative abundance for each species was quantified as the average maximum number of tracks observed per transect (Rees 2015). Because wolves were not included in the track survey protocol until winter 1998–1999 (Rees 2015), we chose to quantify wolf abundance using the Wisconsin DNR mid-winter minimum population estimate rather than the track survey index (Wiedenhoeft et al. 2018). These data come from a combination of extensive radio-telemetry and targeted snow-tracking surveys and thus likely provide improved estimates of wolf abundance compared to the track surveys, which we felt was critical to understanding wolf impacts. These counts extend from winter 1979–1980 to winter 2017–2018, providing a 38-year time-series of wolf, coyote and fox abundance for the northern region. Although track surveys in the southern zone were initiated in winter 1998–1999, continuous surveys did not begin until winter 2003–2004, thus providing us with a 13-year time series in the southern zone.

Statistical analysis

The abundance indices derived from the track surveys appeared capable of capturing long-term trends in population size, but may be insufficient for accurately characterizing year-to-year variation in population growth rates due to high levels of sampling variation and uncertain functional relationships to true abundance (MacFarland and Van Deelen 2011). Because of this, we sought to reduce sampling variation in our coyote and fox abundance estimates by applying a state-space population growth model to our coyote and fox time-series in each region (Humbert et al. 2009). This model uses a stochastic version of the exponential growth model that partitions deviations from exponential growth in year-to-year abundance estimates into both process and observation error, yielding annual estimates of population size that are free of sampling induced variation and instead reflect only biological variation in population growth. We used the fitted values from these models as derived indices of annual abundance for foxes and coyotes. We assumed this approach would help remove sampling variation inherent in the track-survey estimates, for example due to differing sampling effort and sampled transects among years, and provide a more accurate representation of annual fluctuations in abundance (Fig. 3). Indirect estimates such as these frequently assume a linear relationship between the index and true abundance despite the fact that the functional relationship may be more complex (MacFarland and Van Deelen 2011). However, for our purposes we assumed that our derived indices were sufficiently representative of relative changes in true abundance for each species.

Levi and Wilmers (2012) provided an elegant explanation of how linear models of annual population growth rates can be decomposed into density-dependent growth and discrete Lotka–Volterra competition components. Briefly, we consider the annual log difference in our derived population indices (ln [n_t+1 – n_t]) as an additive linear function of population abundance of the focal species (density-dependence) and the competing species (Lotka–Volterra) in the previous year t.

Figure 3.

Annual trends in fox and coyote relative abundance in northern and southern study regions. Dashed black lines represent Wisconsin DNR estimates (Rees Lohr 2017), while gray lines represent predicted values derived from state-space population growth model.

where r_n(t) represents the growth rate of the population in time t, β₀ and β₁ represent the carrying capacity and intrinsic growth rate of the population transformed into regression coefficients, β₂ and β₃ represent the strength of the effects of species 2 and 3 on the population growth rate, and n_2t and n_3t represent the abundance of species 2 and 3 at time t. We refer readers to Levi and Wilmers (2012) for a detailed description and derivation of this model. In this framework, the abundance of a given species in the previous year would be expected to negatively affect its population growth rate in the presence of negative density dependence, while the effects of abundance of others species would serve to measure prey abundance (foxes for coyotes) or competitive effects (coyotes for foxes, wolves for coyotes). To evaluate the relative influence of density-dependent and competitive effects on fox and coyote population growth rates we used an information theoretic approach in which we built a suite of competing models for each species in each region that represented all possible additive combinations of our predictors and evaluated their support based on AICc values and model weights (Burnham and Anderson 2002) while accounting for uninformative parameters in cases of model uncertainty when the null model was within 2 AIC units of the best model (Arnold 2010). We also used model averaging, in which we weighted parameter estimates based on model weights, to determine the final structure of each model (Lukacs et al. 2010). All analyses were conducted in the R programming language (<  http://www.r-project.org www.r-project.org>) using the MuMIn package (Barton 2018). For each covariate we calculated relative variable importance (RVI) as the sum of the model weights in which a particular covariate occurred.

Results

Factors affecting annual population growth rates of canids varied between species and between regions. For fox population growth rates in the northern region, the top three models suggested that by both the abundance of foxes (i.e. density-dependence) and the abundance of coyotes (Table 1) were influential. Collectively, the three models including these terms together and univariately contained nearly 60% of the total model weight for this suite (Table 1). These two covariates each had a negative effect on fox population growth (Table 2). However, the overall strength of these relationships appeared minimal when examining model coefficients and their associated p-values (Table 2). Models for fox population growth in the southern region appeared to be strongly influenced by coyote and wolf abundance based on model weights (Table 1). As with models from the northern region, the strength of this relationship appeared to be minimal when examining model coefficients and p-values (Table 2). In both regions the null model was within 2 AICc units of the top model, suggesting that in each case these may be uninformative parameters. The top individual models of fox population growth in the north and south regions yielded r² values of 0.158 and 0.754 respectively.

Coyote populations in the northern region appeared to be strongly influenced by density-dependence, with the model containing only the density-dependent term yielding 40.5% of the total model weight (Table 1). The effect of coyote abundance in this region on coyote population growth was negative and had the highest relative importance of any parameter (Table 2). In the southern region we found no support for any parameter influence coyote population growth, as the null model was the most strongly supported (Table 1). Additionally, relative importance values for all three parameters were low when compared to other model suites (Table 2). The top model for coyote population growth in the north region yielded an r² of 0.115.

Table 1.

Model selection results for eight competing models of fox and coyote population growth rates in two regions of Wisconsin. Model parameters are based on the abundance of wolves (N_w), coyotes (N_c) and foxes (N_f) in the previous year. Dot notation refers to the null (intercept only) model.

Table 2.

Model averaged parameter estimates, p values and relative variable importance (RVI, sum of w_i values from models in which parameter occurs) from of coyote and fox population growth rates in two regions of Wisconsin.

Discussion

We found limited evidence of a three species interaction chain among canids in Wisconsin. These results are counter to those found by Levi and Wilmers (2012) in neighboring Minnesota. There are substantial similarities between Wisconsin and Minnesota with regards to land use and canid management, thus the exact cause for this discrepancy is unknown. It is possible that the data used in our analysis were of lower quality than those used by Levi and Wilmers (2012), suggesting that our failure to detect a three species interaction chain in Wisconsin could be the result of sampling error. However, that we were able to detect a small effect of coyotes on fox populations, at least in the northern region, suggests that our data were not wholly limiting. Instead, it is possible that the ecological context of our system, relative to Minnesota, may facilitate the coexistence of wolves and coyotes. For example, the wolf population in Minnesota is more than twice the size of the Wisconsin wolf population (Erb et al. 2015, Wiedenhoeft et al. 2018), which may suggest that wolf limitation only occurs or may only become detectable once certain population densities are reached, and that the population in Wisconsin has yet to reach such a threshold. Similarly, a smaller wolf population may only result in localized suppression of coyotes (Berger and Gese 2007) rather than at a broader spatial scale, and there is already compelling evidence to suggest that wolves can impact coyotes at local scales in the Great Lakes region (Flagel et al. 2017). Because the coyote and fox abundance estimates we used were spatially aggregated, we had no means of assessing localized impacts of wolf populations. Additional examination of spatial patterns in canid abundance, particularly as they change over time, may elucidate these potential effects. The most compelling evidence for wolf impacts on coyotes comes from Minnesota (Levi and Wilmers 2012) and Yellowstone (Merkle et al. 2009), areas with substantially higher wolf densities than Wisconsin. Thus, it is possible that our failure to detect large-scale evidence of wolves suppressing coyotes is a function of density-related mechanisms. If wolf populations in Wisconsin continue to increase, such effects may become evident in the future. It is also possible that localized habitat and landscape conditions surrounding our sampling transects mediated canid interactions (Flagel et al. 2017).

In the western United States, wolf limitation of coyote populations is well documented (Berger and Conner 2008, Merkle et al. 2009). One potential explanation for this is that coyotes are somewhat dependent on scavenging of wolf kills in systems where large herbivores, that are too large to be regular prey of coyotes, are the dominant source of food for wolves (Merkle et al. 2009). This leads to direct interactions between wolves and coyotes that can serve to limit coyote populations. In Wisconsin, the dominant prey species for wolves is white-tailed deer Odocoileus virginianus, a species that is also readily killed by coyotes (Ballard et al. 1999). In much of northern Minnesota, where Levi and Wilmers (2012) found evidence of wolves suppressing coyote populations, moose Alces alces can comprise a major component of wolf diets (DelGiudice et al. 2009). Thus, it is possible that in Wisconsin the abundance of a shared prey species may serve to limit direct interactions between wolves and coyotes by reducing the amount of scavenging of wolf kills done by coyotes. A similar phenomenon may also explain why we failed to detect any competitive influences in our southern study region, where our top supported models of both fox and coyote population growth rates were our null models. High primary productivity in Wisconsin's southern region is indicative of high prey biomass in the form of small mammals and other shared prey between coyotes and foxes. Thus, our failure to detect strong interspecific relationships in either study region could be attributable to high prey abundance. However, it should be noted that, in the southern study region, this may be the result of the limited length of our time-series in this region (n = 12 years).

Because there have been no studies conducted in Wisconsin that directly assessed interactions among canids, we urge caution in the interpretation of our strictly correlative results. Although we were not able to find strong evidence to support the plausible biological hypothesis of a three species interaction chain among wolves, coyotes and foxes, there is still a distinct possibility that wolves may be impacting coyote populations through behavioral modifications (Merkle et al. 2009) rather than through a direct influence on population growth. We suggest that targeted studies of these species where they are sympatric in this region could help to identify the mechanisms that support what appears to be a state of coexistence between these two species. Conversely, the mechanism behind the apparent limitation of foxes by coyotes, at least in our northern study region, is unknown. It is possible that this is the result of direct mortality of foxes due to coyote predation (Gosselink et al. 2007), but it is also possible that interference competition between the two species exists in this region (Gosselink et al. 2003). Again, without a targeted study to directly assess interactions between the two species, any inferences regarding their interactions in Wisconsin remain purely speculative.