Introduction

Invasive alien species (IAS) are a worldwide issue and in many cases demand some form of population management. The muskrat (Ondatra zibethicus L.) is invasive in Europe and is for important reasons considered a serious threat to human safety, economy and biodiversity in several countries (EU Deputy Direction of Nature 2015). Its habit of excavating extensive burrow systems in earthen dams (levees), banks and other water retaining structures is considered dangerous, especially in regions below sea level that would be flooded if these structures were damaged (BCM 2006; Bayoumi and Meguid 2011). To minimise these and other negative effects, muskrat populations are actively controlled in some of these countries. The efficacy of this population management has been challenged (Pelz 1996; Zandberg et al. 2011), and in this paper, we report a study designed to compare the growth of muskrat populations with and without a control programme.

The muskrat is a medium-sized rodent, native to North America, with a semiaquatic lifestyle. Basic reviews of its biology are given in Boutin and Birkenholz (1987), Heidecke and Seide (1990) and Perry (Perry Jr. 1982). Natality, mortality and dispersal are all affected by population density and show strong seasonal variation (Errington 1963; Clark and Kroeker 1993). Muskrat populations are present in large parts of NW Europe, following introduction in 1905 (Artimo 1960; Long 2003; Le Louarn and Quéré 2011). Recent European legislation (EU regulation no 1143/2014 on the management of invasive alien species) has identified the muskrat as posing a high risk, and the species is now on the list of IAS of Union concern. Member states are required to take appropriate action if muskrat populations are found on their territories. Derogations from this obligation are permissible only if methods are unavailable or if it can be demonstrated with reasonable certainty that the long-term costs will be disproportionately relative to the benefits of control.

In several European countries, muskrat control is not an issue because population sizes have declined for unknown though apparently natural causes. This has happened in Sweden (pers. comm. P-A Åhlen), Finland (data Finland Natural Resources Institute LUKE, http://statdb.luke.fi), Poland (Brzezinski et al. 2010), Lithuania (Ulevicius et al. 1999; Butautytė-Skyrienė et al. 2011), Russia (Danilov 2016) and Slovakia (Rajsky and Rajsky 2003 in Lammertsma and Niewold 2005). In contrast, muskrat control remains a resource-demanding activity in Flanders (Belgium), the Netherlands and parts of northern Germany. This is because large parts of these areas are below sea level and the integrity of the water infrastructure is of vital importance. The number of animals killed (by mechanical means only; no poison is used) in these programmes is high, amounting until recently to several hundreds of thousands per year (van Loon et al. 2017a). Negative side effects of the control programmes include bycatch of fish, birds and mammals (Bos et al. 2017), disturbance of nature areas and blocking of ecological corridors by muskrat traps.

Several authors have doubted the efficacy of muskrat control programmes, arguing that annual catch rates by trappers were too low relative to recruitment to result in a significant population decrease (Pelz 1996; Barends 2002; Zandberg et al. 2011). In the UK, the effectiveness of muskrat control was demonstrated in five well-documented cases that resulted in complete eradication—but all were small, recently established and isolated populations (Gosling and Baker 1989; Robertson et al. 2017). In the well-established populations in mainland Europe, some evidence for an effect of muskrat control is provided by historical correlations (van Loon et al. 2017a), population dynamic modelling (Bos and Ydenberg 2011; van Loon et al. 2017b) and field experience in Flanders (Stuyck 2008; VMM 2010). But experimental evidence that would allow the effects of trapping to be isolated from influences such as disease, climatic factors, food availability, predation and possibly other factors is lacking.

A true experiment with the random allocation of suspended or zero trapping to replicated sites is impossible in the Netherlands due to the strong aversion of the Regional Water Authorities to possible negative consequences. However, in the context of a nationwide study with a number of parallel studies, we identified nine study sites, all in the Netherlands (2013–2017), in which muskrat control either (i) was ongoing (ongoing trapping), (ii) had been suspended for 3 years (suspended trapping) or (iii) had ceased more than 8 years previously (no trapping). To take advantage of this opportunity, we planned and carried out systematic capture-mark-recapture (CMR; Otis et al. 1978) procedures in each of these areas, in order to estimate indices of muskrat abundance that would enable comparison of demographic measures in these well-established populations. The predictions are that (1) the muskrat abundance index will be higher in sites without muskrat control than in sites with control and (2) the muskrat abundance index will increase after trapping is suspended, altering the index values resembling the ‘ongoing trapping’ situation to those resembling the ‘no trapping’ situation.

Material and methods

Study sites

The Netherlands has a mild maritime climate with average winter temperatures of 2 °C. Water levels are nearly everywhere carefully regulated and stable, with possibly different target levels for summer and winter. The diverse predator community includes terrestrial (fox Vulpes vulpes L.; grey heron Ardea cinerea L.), aquatic (otter Lutra lutra L.; pike Esox lucius L.) and aerial predators (white-tailed sea eagle Haliaeetus albicilla L.). The non-native predator American mink (Neovison vison S.) is also present, almost exclusively originating from escapes from fur farms (Dekker and Hofmeester 2014).

We compared a total of nine sites in the Netherlands on which (i) muskrat control was suspended for 3 years (suspended trapping, ST), (ii) ongoing control by kill trapping was conducted (ongoing trapping, OT) or (iii) control had been abandoned 8 years previously (no trapping, NT). These treatments overlapped during the period from spring 2014 to autumn 2015. In the remainder of the manuscript, we refer explicitly to ‘livetrapping’ to distinguish it from ‘kill trapping’, for which we (also) use the terms ‘trapping’ and ‘catch’.

Suspended trapping

Muskrat trapping was suspended from January 2014 until March 2017 at a site around the village of Dinteloord (51.6° N, 4.3° E) and at a site near and in the town of Lelystad (52.4° N, 5.4° E). The Dinteloord site (3936 ha) is dominated by agriculture, while the Lelystad site (2526 ha) is urban—although with many parks and green areas (see Fig. 1 and Table 1). In the years prior to the cessation of trapping, the mean annual trapping effort in these areas amounted to 480 and 800 h, respectively. During the experiment, a limited trapping effort (< 5% of previous annual totals) had to be maintained to ensure the integrity of strategic dikes and dams in a few specific situations at which no muskrat presence was allowed. Kill trapping was resumed March 20, 2017.

Fig. 1
figure 1

Location of the two experimental areas where muskrat control was temporarily suspended (ST; spring 2014 to March 2017), the six stratified randomly selected reference sites (OT) and the site without any muskrat control (NT). Also indicated is a pilot site at the German border that features in the ‘Discussion’. See Table 1 for details regarding the sites

Full size image
Table 1 Characteristics of the nine study sites. AS codes identify an ‘atlas square’. The column ‘stratum’ refers to categories of atlas squares in the Netherlands (low, medium and high) differing in (i) linear length of waterways, (ii) historical muskrat trapping effort and (iii) historical muskrat catch. The last record refers to a pilot site, visited early in the study to practise field methods, that features in the ‘Discussion’
Full size table

No trapping

The NT site was in the 3600-ha area of marshland in the ‘Oostvaardersplassen’ (see Fig. 1 and Table 1), a nature reserve also containing 4000 ha of grassland grazed by Heck cattle (Bos taurus L.), Konik horses (Equus caballus L.), red deer (Cervus elaphus L.) and geese (Anser anser L., Branta leucopsis B.; Cornelissen et al. 2014). This is the only area in the Netherlands currently without muskrat control. No trapping has been carried out in the marsh since 2006 (van Eerde 2006).

Ongoing trapping

Six sites (see Fig. 1, Table 1), each a 5 × 5 km ‘atlas square’ based on the national Dutch grid reference (Vogelbescherming Nederland 2007), were selected from a large-scale national experiment (Bos et al. in prep.). Pairs of sites located near the towns of Raalte (52.3 °N, 6.2 °E), Sneek (53.0° N, 5.6° E) and Boskoop (52.0°N, 4.6°E) represent, respectively, three strata (low, medium and high presence of muskrat habitat). The three strata (see Bos et al. 2016) differ in three measures that covary strongly and positively, namely, (i) linear length of waterways (see Table 1), (ii) historical muskrat trapping effort and (iii) historical muskrat catch. The length of waterways (km; defined in van Loon et al. 2017a) is an estimate of the amount of muskrat habitat. The historical muskrat trapping effort and historical muskrat catch refer to the mean annual muskrat catch and mean annual effort over 3 years (2009–2011) preceding the start of the large-scale national experiment.

The ‘low’ strata sites near Raalte (OT1 and OT2) have sandy soil, little muskrat habitat (90–91 km of waterway), little annual historical trapping effort (< 100 h y−1) and few muskrat captures. The ‘high’ strata sites near Boskoop (OT5 and OT6) have peat soil, much more muskrat habitat (680–728 km of waterway) and high trapping effort (1900–3300 h y−1). The sites near Sneek (OT3 and OT4) are also on a peat soil and represent the intermediate stratum (with 361–378 km of waterway and 400–760 h y−1 trapping effort).

Capture-mark-recapture

CMR is a standard ecological methodology to estimate survival rates and population densities (Otis et al. 1978), and it enabled us to follow population parameters during the experiment. We livetrapped muskrats on all nine study sites (see Fig. 1, Table 1), fitting each individual with a pair of earmarks carrying a unique identity number (Monel tag #1). At each site, 40 livetrapping cages (20 × 20 × 60 cm) were deployed. Livetraps were placed along water bodies on floating rafts, or directly on banks, using the same locations throughout the study, placed ~ 350 m apart. Professional trappers in and around the study sites subsequently reported marked retrapped muskrats, including those killed or found dead.

In each of the six OT sites, six CMR sessions were undertaken between April 2013 and October 2016. Livetrapping sessions lasted eight nights (Monday evening until Friday morning in two successive weeks, with no livetrapping during the intervening weekend) and took place in April (directly after the spring dispersal period) and October (just before the autumn dispersal period). The livetrapping cages were placed in the centre 9 km2 of each 5 × 5 km atlas square. At the OT sites, the regular trapping programme (using kill traps) was temporarily suspended during the eight-night CMR sessions.

In each of the ST sites, there were seven CMR sessions between spring 2014 and spring 2017. The timing and duration of the livetrapping sessions was identical to that in the OT sites (see above).

In the large NT site, livetrapping was concentrated along large ditches giving small boat access, as the vegetation in the main marsh area was impenetrable.

Here, livetrapping sessions were six nights in duration (Monday through Wednesday nights in successive weeks, with no livetrapping during the intervening weekend) and took place from spring 2014 until autumn 2015 in the eastern marsh (an estimated area of 6.4 ha and 2.2 km of waterway). The autumn 2015 session was shortened at the request of park management, to minimise potential disturbance of a recently established otter. Two additional sessions took place in the autumns of 2016 and 2017 in the western part of the marsh (an estimated area of 4.1 ha and 2.6 km of waterway). The data for this site were thus collected somewhat less consistently than those at the other sites.

Data analysis

Muskrat abundance was assessed from the live capture data in two ways by (1) mark-recapture modelling and (2) using an index (captures/trap night).

We analysed the CMR data with a robust design mark-recapture model (Kendall et al. 1997), using the package RMark (Laake 2013), which is an interface in R for the programme MARK (White and Burnham 1999). We tested robust design models using primary occasions (spring or autumn) and secondary occasions (trap nights within a season) to estimate survival probability between the primary occasions and resighting probability as well as population abundance per primary occasion. We did not take dead captures into account, as the percentage was low (< 1%).

Assumptions of the robust design models are the following: (1) the population is ‘closed’ during the secondary occasions; i.e. no births or mortalities and no immigration or emigration should occur within that time span, and (2) all individuals are correctly identified (Kendall et al. 1997). Because the trap nights spanned 10 days at maximum, we assumed that the population was effectively closed during the secondary occasions. We also assumed no tag loss. Models were fitted separately per site due to differences in the number of livetrapping sessions. If the data allowed, an estimate was made of the number of muskrats prone to capture per session (i.e. within the neighbourhood of the livetraps). On each site, we compared models in which survival and resighting probability could vary between primary occasions with models in which survival and resighting probability were constant over time, and models in which immigration and emigration rates could be random or fixed at zero. Population abundance was allowed to vary between sessions, but the low sample size could not support models in which emigration and immigration differed between sessions.

Model selection was assessed using Akaike’s information criterion for low sample sizes (AICc). We used the primary occasions to test for goodness of fit of the data with release.gof() in RMark for the three sites combined, and we found no overdispersion (χ2 = 7.15, df = 9, P > 0.5). For this analysis, the dataset for the NT site was limited to the four sessions from 2014 to 2015 in the eastern marsh.

Sufficient numbers of live captures were reached for the CMR analysis only at the ST and NT sites. In order to be able to compare muskrat abundance among all treatments, we also calculated for each session an index of muskrat abundance from the number of live captures relative to the number of trap nights, expressed as captures per 100 trap nights. To account for the fact that OT sites were not entirely independent but matched in pairs, values were first averaged per pair and session before modelling. We used a generalised linear model (GLM) with a quasi-Poisson error distribution to test the relationship of the muskrat abundance index to treatment (OT, ST and NT) and season (spring and autumn). We started with a full model including all interactions and used backward selection to remove insignificant interactions until no further improvement was possible. Models were assessed using F tests. Influential points were identified using Cook’s distance, and their effect on the model outcome was studied by replacing them with the 95 percentile value. All analyses were performed in R version 3.4.2 (R Core Team 2017). The same models were evaluated for live captures of unique individuals.

Dead captures and estimates of density

The Dutch muskrat control programme has extensive historical data on catch rates in each atlas square on the national reference grid. As background for the analyses of live captures, the catch rates (expressed as muskrats per kilometre of waterway per year) in atlas squares corresponding to, or overlapping with, each of the study sites are given in Table 2.

Table 2 Background information on catch rate (dead) and live capture rate during the study period for atlas squares corresponding to, or overlapping with, each of the study sites. The catch rate refers the sum of muskrats killed per kilometre of waterway per year (n/km/y). Live capture rates refer to the mean of the total or unique number of individuals per 100 trap nights over the study period (n/100 trap nights). The last record refers to a pilot site, visited early in the study to practise field methods, which features in the ‘Discussion’. Source: Dutch Water Authorities
Full size table

The percentage of marked individuals reported killed by trapping during the experiment was calculated per treatment. We estimated muskrat density in the ST sites at the end of the experiment (20 March 2017) by multiplying the total catch in the 10 weeks after kill trapping was resumed, with the ratio [the number of marked individuals presumed alive at the end of the experiment] to [the number of marked individuals reported dead in the 10 weeks after trapping was resumed].

Data availability

The datasets generated during and/or analysed during the current study are available in the University of Groningen repository, DataverseNL https://hdl.handle.net/10411/YQPYZM.

Results

Mark-recapture modelling, survival analyses and population estimates

Estimates of muskrat numbers could be made over several years at the ST sites in Dinteloord and Lelystad. At both sites, the estimate rose toward the end of the experiment (see Fig. 2). In ST sites, the top model had recapture rates varying per session, constant survival (0.24 ± 0.06 s.e. for both sites) and no parameter for emigration or immigration. Though it was possible to fit a mark-recapture model for the NT site, muskrat abundance could not be estimated. The top model for this site had constant survival (0.37 ± 0.08 s.e.) over the sessions and varying rates of emigration and immigration (0.21 ± 0.16 s.e.; see electronic appendix in Online Resource 1 for tables with model selection results and all parameter estimates).

Fig.2
figure 2

The number of muskrats prone to capture in the two ST sites, (A) Lelystad and (B) Dinteloord, estimated using capture-mark-recapture modelling. The black dots refer to capture sessions for which the model was able to produce an estimate of population size, with the 95% confidence interval. Open dots refer to the number of unique individuals caught during those capture sessions in which no population estimate could be produced

Full size image

Index of muskrat abundance

During the study, we made 935 live captures of muskrats, comprising 289 unique individuals. The capture rate (live captures per 100 trap nights) in the NT treatment was 25.8 ± 2.7 s.e. (n = 1 site, 6 sessions), almost 50 times higher than that in the OT treatment (0.54 ± 0.14 s.e.; n = 6 sites, 6 sessions; see Fig. 3). Neither treatment (NT nor OT) shows any hint of a trend over time. In contrast, in the ST sites, the index of muskrat abundance strongly increased over time, from a level matching that of the OT treatment in 2013 (after which control was suspended) to that matching the NT treatment in 2017 (Fig. 3). The final backward-iterated model contained significant effects of treatment, season, year and interaction treatment × year (Table 3). The model outcome was robust for the replacement of the most influential point. The pattern for live captures of unique individuals was essentially the same (Table 3) with highly significant effects of treatment, year and interaction treatment × year.

Fig. 3
figure 3

Indices of muskrat abundance (live captures/100 trap nights) in nine study sites in the Netherlands. The lines indicate the predicted values from the best model and the grey area refers to the 95% confidence interval over the prediction. The black circle with red background is an additional data point which is not in the model, but which features in the ‘Discussion’. See Table 1 for details regarding the sites

Full size image
Table 3 Parameters included in the backward-iterated GLM Poisson models (which agreed best with the data) to explain variation in live capture rates and their significance. (A) total live capture rate, (B) live captures of unique individuals. The models used a log link function and a quasi-Poisson error distribution. Parameters between brackets are not included in the final models
Full size table

Kill trapping, marked individuals captured dead and estimates of density

Kill trapping in the ST sites in the 10 weeks directly following the end of the experiment (21 March to 31 May 2017) caught 533 (Lelystad) and 706 (Dinteloord) muskrats. This increased threefold or more the 2017 catch rate recorded in the atlas squares including the ST sites (Lelystad 2.19 muskrat km−1 y−1; Dinteloord 1.11 muskrat km−1 y−1) compared with the previous 3 years (2014–2016) during which trapping had been suspended (Lelystad 0.23–0.69 muskrat km−1 y−1; Dinteloord 0.26–0.37 muskrat km−1 y−1; see Table 2).

Similar increases were not observed in the OT sites, where in most cases the 2017 catch rate was the same as during the experiment. Obviously, there are no kill captures in the NT site, but in the zone surrounding the nature reserve, about 0.64 muskrat km−1 yr.−1 (± 0.14 s.e.) were caught annually.

Out of all individuals marked in the OT sites, 40% were killed and reported during regular kill trapping. Only a single marked individual from the ST area (< 1%) was kill trapped outside the boundaries of the study site.

Given the survival estimate of 0.24 per half year in ST sites, the number of marked individuals present at the end of the experimental period was estimated at 21 and 56 for ST sites Lelystad and Dinteloord, respectively (Table 4). The majority of marked individuals was caught and reported after kill trapping was resumed in March 2017: 19 in Lelystad and 38 in Dinteloord. Together, this led to estimated population densities of 2.2–3.0 muskrats per kilometre of waterway at the end of the experiment (spring 2017, Table 4).

Table 4 The number of marked individuals reported dead in the first 3 months after regular kill trapping was resumed in the ST treatment, in relation to total catch in the same period. The estimated number of unique individuals presumed to be alive was calculated from the numbers marked per season multiplied by survival (0.24 per half year). Using the ratio of marked animals present and reported, a population estimate has been made
Full size table

Discussion

The important finding of this study is that muskrat numbers strongly increased after the experimental suspension of control by trapping. At the start of the suspension of trapping, the muskrat abundance index on the ST sites was low and matched that at the sites where kill trapping was ongoing (OT sites). Once trapping ceased, the index rose to match that of the nature reserve site where kill trapping had been abandoned 8 years ago (the NT site). This increase contrasts with the temporal pattern of kill trapping at the six OT sites (see Table 2) and with the declining trend in catches in the Netherlands as a whole (van Loon et al. 2017a). The CMR models estimated considerably higher muskrat numbers at the ST sites by the end of the experimental period, corroborating the above findings. The findings are further bolstered by the large numbers of muskrats trapped in the ST sites once the experiment ended. Parallel increases in catch rate were not observed in OT sites. Thus, in the current Dutch landscape, with the currently prevailing diseases, food availability, climatic conditions and composition of the predatory fauna, a threefold increase in muskrat abundance was measured over 3 years under field conditions when muskrat control was suspended (Fig. 2). These experimental results provide compelling evidence that the correlation between trends in muskrat captures and control effort found by van Loon et al. (2017a) is causal.

The robust design mark-recapture models allowed estimation of change in muskrat numbers prone to capture at the ST sites. Because muskrats are territorial and do not mix randomly, an estimate of the density of muskrats from these numbers depends on their home range and the area sampled. The density of traps was probably too low to cover effectively all home ranges of present muskrats in the entire area. We succeeded however in obtaining a reliable estimate of density using the information on numbers of muskrat trapped after the experiment stopped (Table 4), arriving at an estimate of density between 2.2 and 3.0 muskrat km−1. These values are low in comparison to the range of values published elsewhere (range 4–48 muskrat km−1; Doude van Troostwijk 1976; Leboulengé and Leboulengé-Nguyen 1981; Brooks and Dodge 1986).

That parameters for emigration and immigration could not be estimated in the robust design capture-recapture models for the ST sites is probably a consequence of low sample size. The survival estimates (per half year) are 0.24 ± 0.06 s.e. for the ST sites and 0.38 ± 0.08 s.e. for the NT site which is in the same order as values measured in North America (0.29–0.46; Clay and Clark 1985; Clark and Kroeker 1993). The limited sample size in OT sites is interpreted as a result of the low density of muskrats in areas with muskrat control, given that we invested the same effort as in Clark and Kroeker’s (1993) study in experimental marshes in Manitoba but caught considerably fewer muskrats. Given this situation, we believe the use of the index of abundance (captures/trap night) to make comparisons between the treatments was appropriate as an alternative to mark-recapture modelling.

Muskrat abundance in the Dutch landscape

Generally, all six OT sites had a low muskrat abundance index throughout the study period, indicating low population size. We believe these three pairs of sites are representative of sites on the Dutch landscape with high (OT5 and OT6 near Boskoop), medium (OT3 and OT4 near Sneek) and low (OT1 and OT2 near Raalte) amounts of muskrat habitat. The ST and NT sites were all located in areas with a medium amount of muskrat habitat, and the sites developed or had a high muskrat abundance index in the absence of control. All this suggests that muskrat population sizes are nowadays low at a national scale, in comparison with what would likely happen in the absence of effective control. This is consistent with the finding that current muskrat numbers in the country are low in comparison to the recent past, which was a conclusion in a previous population dynamic model study (van Loon et al. 2017b) and which has also been deduced from the historical catch data (Unie van Waterschappen 2017; van Loon et al. 2017a).

An interesting exception is found near the village of Leuth (51.8° N, 5.9° E), on the Dutch border with Germany. At this site, a field pilot with the CMR method resulted in a particularly high value of the index of muskrat abundance of 20.4 captures/100 trap nights (n = 1 site, 1 session), given that the area is under regular year-round muskrat control (Fig. 3). The location was chosen as a pilot site for practising our field methods, because it was anticipated that chances of live captures would be substantial there. Both field experience and historical catch data suggested that immigration of muskrat is high at this border; the pilot site is characterised by a relatively high catch rate (Table 2; the mean is 0.84 ± 0.08 s.e. muskrat km−1 y−1), considering its situation in the ‘low’ stratum. We presume this is related to higher muskrat densities on the German side of the border, where the emphasis on muskrat control is lower. Our live capture data confirm the high encounter rates with muskrats here and corroborate the conclusion made by several authors (Matis and Kiffe 1999; Bos et al. 2016; van Loon et al. 2017b) that the role of migration cannot be neglected in population dynamics of muskrats.

Consequences for management

Effectiveness of trapping is one of the basic premises of the control programme, and this study delivers experimental evidence that muskrat density is causally related to control. It is consistent with experience gathered in Flanders (Stuyck 2008; VMM 2010) and the UK (Gosling and Baker 1989) where muskrat populations have been reduced effectively using mechanical means. However, other criteria need to be considered as well when evaluating the need for and usefulness of muskrat control. These criteria include a relationship between numbers and serious damage (given in Bos et al. 2016), ethical considerations (Warren 2007; Zandberg et al. 2011) and a lack of feasible alternative methods to prevent this damage other than control. It is recommended to carefully study the feasibility of alternative methods to prevent damage and provide a cost-benefit analysis (CBA: Reyns et al. 2018). In the meantime, it is recommended that the current low densities are valued as the result of a large investment (sensu Clark 2010), which should not be discarded without proper information on the alternatives.

Conclusions

Under the prevailing conditions in the Dutch landscape and given the current investments in trapping, muskrat control limits muskrat population size.