The basic model

We evaluated the alternative harvesting strategies for muskrats using a metapopulation simulation model. We required our model to have a spatial as well as a temporal component, to include the essential demographic processes and to be flexible in comparing the different strategies of population management. We therefore used a stage-structured stochastic dynamic metapopulation model. It was formulated in an existing model environment for reasons of robustness and transparency. Based on a review in Lindenmayer et al. (1995) comparing ALEX, VORTEX and RAMAS, we judged that Ramas Metapop 5.0 would be the most suitable model environment for our purpose, given the requirements presented above. Ramas Metapop 5.0 is well described (Akcakaya 2008), simplifying communication about calculations and technical details of the model.

The model simulates population development in a series of subpopulations and compares the outcome of different harvesting strategies. In each simulation, relevant aspects of the model were changed, such as season of harvesting, location of harvesting and intensity of harvesting (the proportion of animals harvested). The details of the general model are given in the appendix and in Akcakaya (2008). Here, we present the essential choices, such as those related to temporal and spatial scale as well as vital rates. As we are interested in seasonal harvesting, we used time steps equal to one season of three months. We distinguished seven stages to account for juveniles in three seasons and adults in four seasons. Juveniles reach adult size and weight by winter (Vincent & Quéré 1972), so we decided not to model juveniles in spring. Although in muskrats the sex ratio in adults is often biased towards males (Vincent & Quéré 1972, Parker & Maxwell 1984), we assumed a 50% ratio in adults and sex ratio was not modelled explicitly.

Because of our interest in a strategy of harvesting that is differentiated in space, we modelled a metapopulation that started with 49 subpopulations in 49 cells. The cells were arranged in a square of 3.5 × 3.5 km. The centre of each cell was 500 m from the centre of a neighbouring cell (Fig. 1). This level of scale was chosen in relation to known values of home range and population density, as well as ranges of dispersal by muskrats, and scale and grain of the Dutch landscape. Each subpopulation was defined by the same parameters (Table 1), for reasons of simplicity. These parameters are: type of density dependence (see below), survival and recruitment per stage (the stage matrix), initial abundance and, under assumptions of density-dependence, carrying capacity (K) and maximum population growth rate (R_max). Within simulations, there can be variation between subpopulations due to demographic and environmental stochasticity (see Appendix I and Akcakaya 2008). The cells are connected by dispersal. Dispersal to neighbouring cells is defined as the proportion of individuals from a subpopulation that move to another population at a given distance. It is calculated from a distance-dispersal function defined by four parameters, i.e. three parameters determining the shape of the function (a, b & c) and maximum distance (D_max).

The input parameters were obtained from the literature as far as possible (see Table 1). In order to deal with uncertainty, we used a range of values for each parameter involved and grouped parameter values that have positive and negative effects on population viability. We labelled these groups ′pessimistic′ and ′optimistic′, respectively (see Table 1). This short-cut, suggested by Akcakaya (2008), eases the burden of making large numbers of simulations varying each uncertain parameter individually. We took demographic parameters from three different sources (Smith & Jordan 1976, Clay & Clark 1985, Clark & Kroeker 1993). These sources allow the calculation of seasonal recruitment and survival rates, based on mark-recapture estimates over multiple years. We estimated dispersal parameters from data presented in Caley (1987) and Mallach (1971). We fitted the distance-dispersal function to their data by eye. Dispersal between cells is symmetric, and relative dispersal is assumed to be highest in juveniles in autumn and winter (see Appendix I, following Caley 1987). We estimated the capacity (K) based on Leboulengé & Leboulengé-Nguyen (1981), Doude van Troostwijk (1976) and Brooks & Dodge (1986). The estimate of maximum growth rate (R_max) has no reference to a source. It is consistent with the preliminary results of a non-linear regression analysis on capture data for the Netherlands (Bos et al. 2009), to which its values were calibrated. In our basic model simulations, we assume that density dependence operates, affecting all stages and all vital rates. This is in line with published findings (Errington 1954, 1963, Clark & Kroeker 1993). The type of density dependence was assumed to be ′scramble competition′. Density dependence of this type refers to situations where population growth diminishes at high densities of animals due to worsening returns, i.e. all animals receive fewer resources. In order to test the model robustness for this choice, we also ran simulations assuming ′contest competition′ or ′no density dependence′.

Validation

Validation of our models has been done in a qualitative way, by comparing model output with published information (i.e. Parker & Maxwell 1984, Errington 1954) on patterns of population trajectories, comparisons of harvested vs unharvested populations, effects of seasonal harvesting and values of harvest proportions.

Analysis of the model and evaluation of harvesting strategies

Each simulation lasted 40 time steps and was run 100 times. Forty time steps is the equivalent of 10 years. Short to medium time horizons are appropriate when models are used to compare alternative management options (Akcakaya 2008). One hundred runs is considered suitable for our aims. Our model results were compared using the expected minimum abundance, a parameter used in Population Viability Analysis (PVA; Akcakaya 2008). The expected minimum abundance is an index of propensity to decline. It can be used to compare results of two models: the area between two interval extinction risk curves is equal to the difference between the two expected minimum abundance estimates (McCarthy & Thompson 2001). Note that normally ′risk′ is the probability of an unwanted event, but in this case, the event of extinction is not unwanted. Other relevant model outputs are population density and harvest results.

First, we demonstrate the effect of harvesting in different seasons, as it is one of the basic assumptions that such effects are relevant (van Vliet & Lengkeek 2007). Over the whole metapopulation, we simulated a harvest in a specific season at proportions of 10 and 25% of the numbers present, and we compared the results in terms of population viability. We considered proportions of 10 and 25% moderate and intensive, respectively. We then proceeded with the evaluation of alternative strategies defined in our introduction. These alternatives are:

We ran simulations for each strategy at different proportions of harvesting. The year-round sustained effort strategy, at a proportion of 0% is of course equivalent to a no-harvesting scenario (i.e. 4). If harvest is relevant for a given strategy, it takes place during a time-step, before dispersal. We illustrate the effect of harvesting proportion in each strategy on: population density, total harvest and final yearly harvest. The strategies are ranked according to their outcome in terms of resulting reduction in population viability and number of animals harvested. A reduction in population viability is the biological objective of the management. This biological objective should be weighed against other objectives, such as financial costs, social constraints and ethical or practical considerations (Akcakaya 2008). Ethical considerations relate to the fact that animals are killed per se, the way that they are killed, and the rationale behind it (Warren 2007) and to some extent also the sheer numbers. In our study, the number of animals harvested is taken as a proximate for both financial costs and ethical considerations.

General model output

Without harvesting, the metapopulation in our model quickly grew to high levels, with strong seasonal fluctuations and a seasonal maximum determined by K (carrying capacity; Fig. 2). Growth is generally positive in spring and summer but negative in autumn and winter. The range of uncertainty in the parameter values is such that the estimated average maximum levels differ very much between the optimistic and pessimistic parameter groups. The equilibrium situation is reached in approximately four years. When we mention population densities in the rest of the text, we refer to autumn population densities. They are the peak numbers and are relevant for our problem. They are also indicative for the average population over the year.

Harvesting in different seasons

Our model predicts a strong effect of harvesting in the different seasons on population viability. Harvesting in the autumn or summer (square and round symbols in Fig. 3) has less effect on the minimum expected abundance than harvesting the same percentage in winter or in spring (triangles in Fig. 3). In addition, the number of animals harvested is lower when implementing harvest in winter or in spring. So, for a given number of animals killed, the predicted reduction in population viability is greater when harvesting in winter or spring.

Year-round harvesting - effect of harvested proportion

With an increase in harvest proportion, the average autumn population density declines (Fig. 4). The harvest is highest at intermediate intensity of harvesting. Above 15% harvesting per time step, the population is overharvested from an exploitation point of view, given the parameters of growth and carrying capacity chosen in this model. For convenience, we classify proportions of harvesting of up to 15% as ′moderate′ and proportions > 15% as ′intensive′. Note that the average harvest over the entire run of 40 time steps differs from the harvest in the final year, especially at high proportions of harvest. This difference is caused by the difference in initial and equilibrium abundance of muskrats.

Alternative strategies - space or time differentiated

The resulting population density in the central population for the year-round and the alternative strategies as defined in our introduction is shown in Figure 5. The central population is a population of interest in each of the strategies, including the space-differentiated strategy. Each strategy has an impact on the population density, but the effect of harvest proportion is of overriding importance. Under time-differentiated harvesting, the autumn population density is higher than under year-round harvesting. This shows that it is not useless to harvest in summer and autumn, although the effectiveness per animal was shown to be low for these seasons (see Fig. 3). Under space-differentiated harvesting, the population in the central cell is also reduced, but not as much as with year-round harvesting at the same intensity (see Fig. 5). For the space-differentiated strategy, it is of course useful to examine muskrat density in the remainder of the metapopulation as well. This is shown in Figure 6. Where high harvest proportions under the year-round strategy result in very low population densities both in the central and the metapopulation. This is not the case in the space-differentiated strategy, where the metapopulation is predicted to have a high population density, close to what would be achieved under no harvesting.

Evaluation of strategies

In Figure 7, we ranked the major alternative strategies on the basis of 1) our biological objective, i.e. the predicted reduction in population viability, and 2) an approximation of the ′costs′ associated with each strategy, i.e. the final yearly harvest. The graph allows us to jointly consider the two criteria (Akcakaya 2008). Note that the space-differentiated scenario is not properly captured in this diagram; one also has to consider the density or viability of the population in the main area of interest, where harvesting takes place.

We need a selection criterion to be able to select the most appropriate alternative. These selection criteria can be simple, e.g. lowest viability or lowest cost. But one could also include constraints in the consideration, which can be drawn into the graph as horizontal or vertical lines. Finally, one can use the cost:benefit ratio as a selection criterion, while satisfying either the cost or the viability constraint. This is matter of public debate, but we will elaborate upon it to some extent in the discussion section (below).

As we deduced already from Figure 5, the strategy of year-round harvesting is predicted to result in the greatest reduction in average autumn population density and thus in population viability (see Fig. 7). Now that we have concluded that harvest proportion has such an overriding effect, we felt obliged to represent each strategy for intermediate proportions of harvesting and for high proportions of harvesting. Of course these are only two out of a range of possibilities. With increasing harvest proportions, the viability is further reduced in each strategy (see Fig. 7, in which the filled symbols are above their open counterparts). For the year-round strategy, the final yearly harvest also decreases with harvest proportion (see Figs. 4 and 7, and compare intermediate to high harvest intensity year-round). But in the space-differentiated strategy as well as in the time-differentiated strategy, the final yearly harvest increases when increasing the harvest proportion to 25% (see Fig. 7).

The year-round strategy at high harvest intensity (25%) leads to the greatest reduction in viability. If we agree that the number of animals harvested is a good approximation of ′costs′, the year-round strategy at high harvest intensity also has the lowest cost:benefit ratio. This is indicated by the line that connects this point to the origin (see Fig. 7; all other points lie below this line). Note that the year-round strategy at intermediate harvest intensity (10%), results in the highest cost:benefit ratio, in comparison to the other alternatives. To use financial costs as a criterion, one has to assume a relation between the financial costs of each strategy and the amount of animals harvested. This is a very critical relation that we will discuss below.

Robustness

The ranking of alternative strategies is robust to different assumptions of density dependence (Table 2), and the ranking does not change when pessimistic or optimistic input parameters are used in the model.

General

Spatially explicit models are important tools for investigating scale-related questions in population ecology (Dunning et al. 1995). Our paper is an example of the application of such a model in the context of pest management. In the case of the muskrat, sufficient information is available (e.g. Errington 1963, Perry 1982, Boutin & Birkenholz 1987, Heidecke & Seide 1990), and the model provides us with a framework into which this information can be incorporated in a transparent way. Our objective was not to precisely predict population trajectories or abundances, and we strongly emphasise that the model results, as with any model, need to be interpreted with substantial caution. Nonetheless, as a tool to think quantitatively, to specifically formulate hypotheses and research questions, and as a means to communicate effectively, quantitative models are indispensable. The model we present here conforms to these objectives.

Robustness of results

All assumptions rank the management options in a similar order, in terms of population viability. In that sense, our model is robust. But this does not mean that all other model output, for each of these runs with different assumptions than our best estimates, is credible. For example, under the assumption of no density dependence, the metapopulation size keeps increasing at an exponential rate and so does the yearly harvest (data not shown). This pattern of unchecked population growth has not yet been observed in any other vertebrate species than humans (Allee et al. 1949); the evidence that density-dependent mechanisms affect the population dynamics of muskrats is strong (Errington 1954, Simpson & Boutin 1989, Clark & Kroeker 1993, Erb et al. 2001).

Comparisons with field examples - validity of the model

The model provides qualitatively realistic results, though of course quantitative validation is required. Population trajectories in the literature show great seasonal fluctuations (Errington 1954, Clark & Kroeker 1993) and so does our model. In the absence of harvest, the densities may become very high, but they are bounded by a maximum (Errington 1954, Clark & Kroeker 1993, Danell 1996). This maximum may vary strongly between habitat and years, by orders of magnitude. Such great variation in carrying capacity was not implemented in our model, and in that sense, our results should be regarded as somewhat idealised. On the other hand, the Dutch habitats are also more stable than the prairie potholes and river systems in the U.S., where such great annual fluctuations have been observed, as they are carefully managed for agricultural production, and the climate is moderated by the presence of the sea. Generally speaking, population models predict that harvested populations have lower average densities than unharvested populations (Boyce et al. 1999). Our output is in line with that finding, but for muskrats, no rigorous field studies are at hand to support it, except the work of Errington (1963) on fur refuges in Iowa. Errington (1963) stated that within the refuges, densities are generally higher than outside. There is one National Park in the Netherlands where a core area is left unharvested, but no data on population numbers have been collected from this area so far. Parker & Maxwell (1984) report on an experiment with controlled harvesting in different seasons. Their results, in line with our predictions, show that combined harvesting in spring and autumn leads to stronger effects upon the population than harvesting in either autumn or spring alone.

Reported harvest proportions vary between 50 and 90% of the autumn population (Smith et al. 1981, Parker & Maxwell 1984, Clay & Clark 1985). Some of these proportions are suspect, because the underlying estimate of the population suffers from difficulties (Boutin & Birkenholz 1987). Parker & Maxwell (1984) harvested about 60% either in autumn or spring, without an obvious decline in population or harvest. Harvesting at those proportions in both seasons caused the population to decline. In autumn, 74% is suggested to be the optimal sustained yield according to a model by Smith et al. (1981). Most of these values can be compared directly to the proportions we used for the autumn, because the literature estimates often refer to animals that were harvested during the autumn trapping season. Our model suggests only a moderate decline in population viability relative to an unharvested situation, at an harvest proportion of 25% in autumn alone (see Fig. 3). A year-round harvest at 35% per time step results in a strong decline in our model. During autumn and winter alone, this corresponds to a compounded harvest of 56% of the autumn population. These two observations are not inconsistent with the above mentioned proportions in the literature. For the other seasons, the comparison is problematic for one reason, because in our study the proportions refer to one time step rather than a full year, and the animals harvested in spring and summer are to some extent born in the year following the autumn that we should refer to.

Current muskrat management

The current situation in Flanders (Belgium) is probably a good representative of year-round intensive management at high harvest proportions. The catch history indicates that muskrat densities in Flanders are probably much lower than in the Netherlands (i.e. in 2008). In Flanders, muskrat management has strongly intensified over the past 15 years. The only reliable data have been collected as of the year 2000, when > 39,000 animals were harvested by the Ministerial agency, the Vlaamse Milieu Maatschappij (VMM), along 2,866 kilometers of major waterways. Other agencies, private parties and the local councils were also campaigning against muskrats at that time, killing an unrecorded amount of animals. By the year 2008 the VMM was the sole responsible party harvesting muskrats. The catches had been reduced during 2000-2008, whereas the area searched had increased dramatically. Currently, the managed populations are presumably small in Flanders. All the effort, which is still year-round and intensive, yielded 5,300 animals in 2008, which corresponds to 1.9 animals.km⁻² or 0.4 animals.km⁻¹ of shorelength. There are no claims of damage and more or less independent checks of population density by controlling teams of trappers also point at low densities. Interestingly, a large proportion of the catches is made in 5-km buffer zones within the boundaries of Flanders. Current labour investment for muskrat management in Flanders (Belgium) appears not to be significantly higher per unit area than for the Netherlands (12 man hours.km⁻² in Flanders vs 14 ± 2.1 (S.E.) man hours.km⁻² (N = 12 provinces) in the Netherlands).

The current Dutch strategy is also a year-round sustained effort. There has been quite some variation in harvest and effort, between different regions within the country over the past 20 years. In general, the current situation in the Netherlands can be considered an intensive management at high harvest proportion and consequently low population sizes. But there are areas where harvest proportion actually appears to be intermediate rather than high. On average, 8 ± 21 animals are caught per km² or 1.1 ± 4.8 animals.km⁻¹ of shorelength (Bos et al. 2009). This is considerably more than in Flanders, but similarly to Flanders, the Dutch trend is declining.

We believe it is possible to increase harvest proportion everywhere in the Netherlands to such a level that populations decline to publicly agreed low levels, without an exceptional rise in costs at the medium or long term. This belief is supported by experience from Dutch and Flemish professional trappers (H. Post, pers. comm.), and by data series of effort and catches in the Netherlands for different provinces and habitat types (Bos et al. 2009).

Role of predation

Interestingly, muskrat numbers in Poland have strongly declined without an intensification of harvest effort. A recent analysis by Brzezinski et al. (2010) documents the decline in the hunting bags of muskrats in Poland after the invasion of American mink Neovison vison and provides an analysis that identifies mink predation as one of the most important factors affecting muskrat numbers. Note that this refers to a situation in which the muskrat is unadapted to mink predation (c.f. Soper & Payne 1997). In mainland North America, American mink are also primary predators of muskrat, but they are especially successful when muskrats are stressed for reasons of food-shortage or disease (Shier & Boyce 2009). In the Netherlands, no viable population of American mink has established yet, in spite of a more or less continual release of individuals from fur farms over the past decades (Muskens & Dekker 2010). The chances of future establishment of American mink are limited. The American mink is not protected and can legally be killed. Besides, mink farming in Germany and Belgium has come to a halt, and maybe this will happen in the Netherlands as well. The native counterpart of American mink, the European mink Mustela lutreola is extinct in the major part of its ancient range. Another potential predator of muskrat is the European otter Lutra lutra, formerly extinct in the Netherlands but now being reintroduced. At the moment, the role of predation in population regulation can only be small in comparison to the effects of trapping.

Management alternatives

In deciding upon the most promising alternative strategy, one has to consider harvest, effort, population viability and average population size (given in Figs. 5 and 7), but also practical feasibility of a strategy, the damage level that is publicly acceptable and the effects upon other natural values. These will be discussed below in relation to the different alternative strategies.

Conclusions

Our model points at differences in effectiveness of the alternative strategies and allows us to evaluate these differences quantitatively with respect to muskrat numbers. Now, in order to also quantify cost effectiveness of each alternative, increase precision and validate the model, clearly identified field information needs to be gathered. Our model thus clarifies the need and usefulness of field experiments dealing with harvesting strategies. It also allowed us to focus such experiments with respect to spatial scale, the variables that need to be manipulated and the parameters that need to be measured, which was our main objective.

We need insight into the costs of harvesting at different proportions for different seasons, landscapes and population densities. We also need to establish the relationship between population density on the one hand and economical damage or safety risk on the other hand. A final constraint that needs to be quantified is the publicly acceptable level of damage per region of interest. These aspects are identified as the major gaps in knowledge that hamper proper policy making at the moment. The benefits that can be derived from guiding expensive control programmes like these with information derived from well-designed field studies, are likely to outweigh the costs of such research.