Population dynamics models can be used to evaluate the effectiveness of predator control. On a restricted area, one key process is the rate at which removed individuals are replaced by immigration. Since this rate is difficult and costly to estimate by field study, we develop an analytical method to approximate immigration rate that makes use of data obtained through the removal process itself. In Britain, red fox Vulpes vulpes control is undertaken by gamekeepers on privately-owned shooting estates. The fox cull on each estate derives from both local reproduction and immigration. The proportional contribution of immigration to the cull can be expected to be greater on smaller estates. We describe a mechanism by which the average annual cull per unit estate area on a very small estate approximates the annual rate of immigration. We used fox culling records from 534 estates across seven different landscape types and a Bayesian hierarchical model to relate the density of foxes culled to estate area, with immigration rate assumed to be equal to the model intercept. The posterior predictive distribution of annual immigration rate was lognormal with a median of 2.41 fox km-2 year-1 and a CV of 0.84. Posterior median estimates of immigration rate varied between landscapes, ranging from 0.86 to 4.13 fox km-2 year1. Immigration rate was higher in arable and pastural landscapes compared to upland landscapes. Variation in immigration rate broadly matched differences in fox density characteristic of the regional landscape type. This study presents a widely applicable method for quantifying immigration rate in populations that are subject to depletion, e.g. through culling. The use of the fox immigration rate estimate as an informative prior distribution in population dynamics models could help in evaluating effects of control on local fox populations and lead to improved control strategies.
Red fox Vulpes vulpes populations are culled throughout much of their range to reduce their impact on prey species and the incidence of zoonotic disease (Macdonald and Reynolds 2004). Across a variety of landscapes in rural Britain, intensive fox control is typically associated with shooting estates aiming to generate a harvestable surplus from small game populations (Tapper 1992, Reynolds and Tapper 1996, Heydon and Reynolds 2000a). Potentially, the local impact of culling on fox density can be determined from data on culling effort and success using population dynamics models. However, as they are rarely recorded systematically, cull data are often characterised by short or sparse time series. This makes it challenging to fit models as the data may contain little information about key model parameters. In such data-poor situations, use of a Bayesian modelling framework can help improve convergence and reduce uncertainty by constraining parameter estimates within reasonable biological limits through explicit incorporation of prior knowledge on model parameters via informative prior probability distributions (McCarthy 2007).
Populations are subject to births, deaths and movement in and out of the population. Culling mortality may either be ‘additive’ to natural mortality or ‘compensatory’, where a removal by culling causes increased survival, reproduction or immigration into the population. At a local scale, replacement of culled animals by compensatory immigration may be a key demographic process in culled populations, with strong implications for control strategy (Smith 1985, Hone 1994, Lieury et al. 2015). Replacement of culled animals through neighbouring foxes expanding their territories into the culled area, or by movement of foxes originating further from the culled area, will both be perceived as immigration on the culled area. If the rate of immigration from outside a culled area is high then culled animals may be rapidly replaced (Reynolds et al. 1993, Robinson et al. 2008). Culling has both economic and animal welfare costs and these become difficult to defend if the impact is insufficient to achieve the aims of culling. Where replacement is rapid, more culling effort will be required to achieve the same reduction in average density than where replacement is slow. Being able to make reliable inferences about immigration rate can therefore help to determine the culling effort required for effective control of a local population.
Despite its importance, immigration rate is notoriously difficult to estimate for a variety of reasons. For example, immigrants must be distinguished from resident individuals and there exist only labour intensive field techniques for this (Abadi et al. 2010). As for most other carnivore populations, little knowledge exists about immigration rates in fox populations, or about how factors such as habitat and surrounding fox density may affect immigration rates. Methods most often used to estimate immigration in terrestrial bird and mammal populations rely on marking individuals. Immigration rates may be estimated in long-term intensive studies where all individuals in the population are marked each year, allowing a direct count of unmarked individuals in the next year (Slough and Mowat 1996); or from the collection of capture—recapture data on marked individuals and, depending on how the data were sampled and which age classes were marked, using different models to infer immigration (Nichols and Pollock 1990, Peery et al. 2006, O'Hara et al. 2009). Capture—recapture data may be combined with population counts or indices to estimate immigration using integrated population models (Besbeas et al. 2002, Abadi et al. 2010). It is also possible to use molecular genetic techniques to detect immigration by determining genetic differentiation between multi-locus genotypes of individuals in different areas (Wandeler et al. 2003).
There are several problems with using these methods for a cryptic species such as the fox. Live capture is difficult, so intensive effort is required to mark or take genetic samples from individuals. In Britain, the approved capture methods for foxes are free-running neck snares and cage traps, but both methods have capture rates below 1 capture/100 trapnights, and usually an order of magnitude lower (Baker et al. 2001, Short et al. 2012). Consequently, trapping for capturerecapture studies must take place over such an extended time period that model assumptions are violated (Sadlier et al. 2004). If capture causes trap shyness, then recapture rates may be even lower. Finally, estimation of immigration rate using marked individuals is scarcely feasible in intensively culled populations because of the rapid turnover of individuals. An alternative way of estimating immigration rate using data readily available from culled populations is needed.
In Britain, many shooting estates volunteer data to the National Gamebag Census (NGC) on the annual numbers of game and predator species killed, including red fox (Tapper 1992). Fox culling is typically carried out by gamekeepers employed by an estate and is confined within the estate boundaries, hence this is ‘restricted-area culling’, sensu McCullough (1996). The majority of estates cover an area of <10 km2 (unpubl. NGC data, GWCT), which we describe as a local rather than regional scale. While the NGC is a valuable indicator of temporal trends in regional fox numbers (Battersby 2005), cull data may also contain information about the rate of replacement of culled foxes by immigration at local scales.
The fox cull on an estate can include both net local reproduction (defined as births minus non-culling deaths) and net immigration (defined as immigrant minus emigrant foxes). The contribution of local reproduction and immigration to the cull will depend on the size of the estate and the type of landscape in which it is located. Cull size will be positively related to estate area. The relationship is expected to be nonlinear asymptotic because 1) the number of foxes available to be culled is finite; 2) the ratio of boundary to area increases as estate area decreases; and 3) the cull size within different landscapes reflects both the number of foxes available to cull and the culling effort applied, such that as estate area increases, the average effort of a single gamekeeper is spread more thinly until there is a clear need to employ an additional gamekeeper. On smaller estates, the increased boundary relative to area will mean proportionally more immigration than on larger estates. Smaller estates have fewer breeding dens (Hewson 1986) and are subject to higher culling intensity, meaning that compared to larger estates in the same landscape, local reproduction will be both smaller and more likely to be culled. As estate area decreases, immigration can therefore be expected to contribute more to cull size than local reproduction, depending on the landscape type and the carrying capacity of the estate.
The immigration rate onto an estate is assumed to be primarily a function of fox density and productivity in the surrounding region, a high immigration rate requires a plentiful source population. In the source population, both fox density and productivity are ultimately determined by food resources, but culling intensity also has a complex influence. At relatively local scales, culling may have an impact beyond the estate boundaries when individual foxes are active on both sides; and culling on neighbouring estates will directly reduce immigration rate by reducing the density of potential source populations nearby. The intensity of culling across larger regions might also be expected to reduce immigration rates. In Britain, fox culling intensity is influenced by regional variation in culling practices (Tapper 1992, Heydon and Reynolds 2000a), with culling found to be a key determinant of fox density at this scale (Heydon and Reynolds 2000b, Heydon et al. 2000). In contrast, variation in fox density between seven landscape types across Britain was found to be closely related to habitat variables in six of them, while relationships of density with culling indices in these landscapes were ambiguous (Webbon et al. 2004). We discuss the apparent contradiction between these studies later, but here we note that there is evidence that habitat-related variables (e.g. food resources) and culling intensity are correlated, with both being important determinants of regional fox density and therefore of immigration rate.
In this study, we developed a conceptual model describing how immigration rate may be estimated from cull data, and used annual fox bag records from the NGC to quantify immigration rates of foxes onto British estates. Although some of the variation in annual fox bags from the NGC was expected to be associated with variation in local culling intensity among estates, we hypothesised that some variation would be explained by landscape-related factors (e.g. food resources or culling intensity) acting regionally through immigration rate. Hence, we aimed to evaluate differences between landscape types. The purpose of this approach was to construct an informative prior for immigration rate, for subsequent use when fitting population dynamics models to assess the impact of culling on populations of managed species.
Material and methods
NGC data and landscape categorisation
We queried the National Gamebag Census database (Game and Wildlife Conservation Trust, Fordingbridge, UK; extracted on 15 April 2016) for shooting estates which submitted annual fox cull records during 1996–2000. This spanned the survey period of Webbon et al. (2004), allowing comparison of our immigration rate estimates with their estimates of fox density in different landscapes from faecal density. A total of 534 estates submitted data for at least one of these years (mean 3.39 years of data per estate). Data were summarised for each estate by calculating the mean area and mean number of foxes killed per year. This dealt with those estates which did not contribute data consistently (327 estates did not contribute for all five years) and those estates for which the area changed over time as new parcels of land are bought, sold or otherwise incorporated into the managed area (48 out of 1166 consecutive annual records showed a change of estate area from the previous year). Use of mean data also smoothed the unknown differences in fox culling effort during this period.
Each NGC estate was classified into one of the seven landscape categories used by Webbon et al. (2004): arable a, arable b, arable c, pastural a, pastural b, marginal upland, and upland. These landscapes were determined by groupings of land class strata which summarise variation in ecological, physio-geographical, and human geographical attributes across Britain (Bunce et al. 1996a, b, Walsh and Harris 1996). Estate records included area and a grid reference; estate boundaries were unknown in most cases. For all estates a landscape categorisation was determined within a circular buffer equal to the estate area and centred on the grid reference provided, using a geographical information system (MapInfo Professional 9.5, Pitney Bowes Software Ltd. 2008). Some larger estates included more than one landscape type: for those estates the category with the largest proportion of total area was chosen. For 27 estates having digitally mapped boundary data, the landscape categorisation determined by the boundary was compared to the categorisation determined by the circular buffer using compositional analysis (Aitchison 1986, Aebischer et al. 1993). This analysis was performed using the ‘compana’ function from the adehabitat package (Calenge 2006) in the R statistical software (< www.r-project.org>).
Conceptual model
We represent the spatial distribution of fox abundance on estates by a 2-D grid of values at discrete points, where the abundance at each point is assumed to be representative of a cell around the point. If Nij is the abundance in the row i, column j cell of the grid, and k represents the index combinations of the four neighbouring cells (i — 1, j; i + 1, j; i, j — 1; ; i, j + 1), the derivative of population abundance with respect to time in cell ij is
Here, R is the per capita birth rate, M is the instantaneous non-culling mortality rate, vk,ij is the rate of immigration from cells k into cell ij, hij,k is the emigration rate from cell ij across cells k, Eij is the culling effort in cell ij and p is a parameter determining susceptibility of foxes to culling effort. Similar differential equation systems are widely used in the fisheries literature to represent spatial mixing of harvested fish (Walters et al. 1999, Walters and Martell 2004 p. 279). Thus, the first two terms amount to net local reproduction rate, the third and fourth terms to net immigration rate (ignoring density-dependent migration between cells), and the last term to the cull rate.
Small estates are expected to have limited net local reproduction as there is less productive area, i.e. breeding dens and feeding habitat. The mean density of foxes across rural Britain is just over 1 fox km-2 (Webbon et al. 2004). Given this relatively low value, as the estate (cell) area approaches zero, the value for Nij can be expected to become very small and decrease monotonically to zero. Immigration can occur over very short timescales, with intensive radio-tracking data suggesting that foxes exploit newly undefended areas adjacent to their own territory after about one week (unpubl. data, GWCT). On this time-step the mean value for R, the averaged weekly per capita birth rate, is 0.054 cub fox-1 week-1 (Porteus 2015). Over the same time-step, the mean value for M, the weekly instantaneous non-culling mortality rate, is 0.007 week-1 (Porteus et al. 2018). Given Nij on a small estate is very small, the terms in Eq. 1 relating to emigration, 
local births, RNij, and non-culling mortality, MNij, can therefore also be expected to be negligible, so it is reasonable to assume that any fox cull on small estates derives primarily from immigration from over the boundary. The rate of population change, dNij/dt, on a very small estate can thus be approximated by removing emigration and net local reproduction from Eq. 1, leaving
Here, 
is the average immigration from neighbouring cells k. Therefore, as estate area approaches zero, and the rate of population change is close to zero, the immigration rate on a very small estate could be expected to be approximated by the cull rate.
Cull rate on a very small estate can be estimated from NGC fox culling data. To remove any effect of estate area on annual cull size, we standardised the cull size by estate area and examined the relationship between cull density (i.e. foxes killed km-2) and estate area. Cull size and estate area are positively related. If this relationship was linear, the relationship between cull density and area would be horizontal; whereas if as expected it was non-linear asymptotic, the relationship between cull density and area would be negative. On a hypothetical estate that has an area approaching zero, the annual cull density could therefore be expected to approximate annual immigration rate per km2. This suggests that we can estimate the rate of replacement by immigration as the intercept of a regression model of cull density on estate area.
Estimation model
Estimation of immigration rate from the NGC data used a regression model relating the fox cull density (D = number of foxes killed annually / estate area) to estate area A. The estimated intercept was assumed to be equal to the annual immigration rate (as detailed above). The regression model also assumes population equilibrium over time, i.e. that the number of foxes killed annually on each estate is sustained by net local reproduction and immigration. We undertook initial data exploration using linear and non-linear regression models. This suggested that a better fit to the data was achieved using a log-linear model with ln(D)) as the dependent variable, which also enabled constant variance assumptions to be satisfied. A desirable feature of using a log-linear model is that the estimated intercept is in natural log-space, meaning that once transformed into the same real-space as the cull density, immigration rate will take only positive values. A log-log model was not considered because the estimated intercept is undefined in real-space for zero area.
A Bayesian hierarchical framework was used to account for effects of landscape category on immigration rate. In a hierarchical model, the regression parameters for each category are assumed to come from a common cross-category probability distribution specified by hyper-parameters that describe the cross-category variability in the regression parameters (Gelman et al. 2004, McCarthy 2007). The observation errors, ε, in the ln(D)) estimates were assumed to be normally distributed, giving the model:
where Di is the cull density on the ith estate and acat and bcat are the intercept and slope parameters for each landscape category. The observation errors are assumed to be constant across landscapes. Immigration rate, v, as fox km-2 year-1, in each landscape was then computed by transforming the intercept values within the model:
Vague priors were used for all parameters to minimise the influence of the priors on the parameter estimates. The standard deviation in the normal observation errors was assumed to be uniformly distributed between 0 and 10. The priors for a and b were assumed to be normally distributed with mean (µa, µb) and standard deviation (σ,a, σb) hyperparameters. The mean hyperparameters were assumed to be normally distributed with mean 0 and standard deviation of 1000, the standard deviation hyperparameters were assumed to be uniformly distributed between 0 and 10. Sensitivity to these prior specifications was tested by varying the standard deviations by one order of magnitude in either direction. The hyperparameters for the intercept were used to describe the posterior predictive distribution for the intercept. Following transformation into real-space, this predictive distribution is suitable for use as an informative prior for immigration rate where landscape type is unknown.
Bayesian analysis was performed using WinBUGS 1.4 (Spiegelhalter et al. 2007), implemented from within R using the R2WinBUGS package (Sturtz et al. 2005). BUGS code can be found in the Supplementary material Appendix 3. Samples from the joint posterior distribution were obtained using Markov chain Monte Carlo (MCMC) simulation. The posterior was estimated from two chains of 25 000, following an initial burn-in of 50 000 samples. Convergence of the Markov chains to the posterior distribution was diagnosed using the R coda package (Plummer et al. 2006). Gelman Rubin convergence statistics (Gelman et al. 2004) were <1.1 for all parameters and Geweke's Z-scores (Geweke 1992) did not fall within the extreme tails of a standard normal distribution, suggesting that the chains had fully converged.
The estimation model assumes a positive relationship between the cull size and area variables. Standardising data on cull size for area effects to account for larger cull sizes on smaller estates can result in spurious correlation, because the cull/area variable (cull density) is inevitably correlated with area (Jackson and Somers 1991). Compared to the situation where a relationship between cull and area exists, the fitting of a regression model to standardised data when there is no relationship between cull and area will therefore result in larger intercepts and more negative slopes due to larger cull density values being found at small areas. We examined the possibility that our analysis using cull density resulted in spurious correlations that influenced the results but we found no significant effects (Supplementary material Appendix 2).
Relationship of immigration rate with fox density
Fox density estimates for each landscape were taken from Webbon et al. (2004). The influence of fox density on immigration rate across landscapes was explored by linear regression modelling. The model intercept was fixed at zero because immigration onto estates within a landscape cannot occur without foxes being present in the surrounding region.
Results
NGC data and land categorisation
Compositional analysis found no difference between the landscape compositions of estates calculated using either the known estate boundary or a circular buffer centred on the estate location (Wilks' A = 0.705, 
= 9.43, p = 0.15, or p = 0.11 by randomisation). In the absence of known boundaries for all estates we therefore used circular buffers to determine the landscape composition. There were unequal numbers of estates within each of the seven landscape categories, and although estates were located throughout England and Scotland, few estates were in Wales (Fig. 1). For 72% of estates, the geographically nearest estate lay within the same landscape type, implying that the containing region was to a large extent confounded with landscape. Estate size varied among landscape types, with mean areas of 8.2 km2 (pastural a), 8.6 km2 (arable a), 9.7 km2 (arable b), 12.8 km2 (arable c), 19.6 km2 (pastural b), 28.9 km2 (marginal upland) and 54.1 km2 (upland). Relative to the total area of each landscape category, there were more estates in upland and arable a categories, and fewer estates in the pastural categories.
Relationships of the fox cull to estate area
In each landscape, mean annual cull was positively correlated with mean estate area, though the correlation coefficient was only weakly positive in arable c, marginal upland and upland (Fig. 2). The relationships between mean cull density and mean estate area were all negative but the correlations were relatively weak (Fig. 2). Correlations between ln(D)) and area were generally stronger, providing an additional reason for using a log-linear model. The slope of the relationship between ln(D)) and area in each landscape was negative, with the 95% credible interval overlapping zero only in the arable a landscape (Table 1). The standard assumptions of linear regression all appeared to have been satisfied. A log-linear model fitted using MLE to the observed data pooled across landscapes was significant (p < 0.001) and had an explanatory power of 27.1% (adjusted R2).
Immigration rate estimates
Posterior probability distributions for immigration rate in different landscapes were very different, with those from arable a and arable b landscapes encompassing significantly greater values than marginal upland and upland landscapes (Table 1, Supplementary material Appendix 1 Fig. A1.1). Posterior median estimates of annual immigration rate for different landscape categories ranged from 0.86 to 4.13 fox km-2 year-1 (Table 1). Coefficients of variation (CV) ranged from 0.11 to 0.22. Posterior distributions were more precise from landscapes with data from more estates. The results were not sensitive to alternative prior specifications. The posterior predictive distribution was lognormal with a median of 2.41 fox km-2 year-1 and a CV of 0.84. This median value was close to that obtained from the model fitted to data pooled across landscapes, which was 2.38 fox km-2 year-1.
Figure 2.
Relationships between (a) mean annual cull and estate area, and (b) mean annual cull density and estate area within each landscape category (with associated Pearson correlation coefficients). Numerals on the right of each row corresponds to the landscape category (I = arable a, II = arable b, III = arable c, IV = pastural a, V = pastural b, VI = marginal upland, VII = upland).
Relationship of immigration rate with fox density and landscape type
The pattern in estimated immigration rates onto estates in different landscape types broadly matched the differences in fox density at the landscape scale, as arable and pastural landscapes had both higher fox density and immigration rates when compared to upland landscapes (Fig. 3). The significance of this positive relationship (p < 0.005) was dependent on the intercept being set to zero (the adjusted R2 of the model is not reported for this reason).
Table 1.
Summary statistics from the posterior probability distributions of lognormally-distributed immigration rate (v) and normally-distributed slope (b) of the relationship between In(cull density) and estate area in different landscapes and the posterior predictive distribution as determined by the hyper-parameters of the model.
Discussion
This study describes a novel method for quantifying immigration onto restricted areas during a culling operation, using the cull data themselves. The method assumes that the contribution of immigration to the cull increases with decreasing size of the culled area, in this case a shooting estate. This is reasonable because on smaller estates there is less area to support resident animals and more boundary (relative to area) across which immigrant animals can move. On small estates, intensive culling will ensure no within-estate reproduction, hence the average annual cull over a period of years must reflect immigration. On an extremely small estate, the annual cull rate was thus considered a suitable approximation of the annual rate of immigration.
Figure 3.
Relationship between posterior estimates of immigration rate in each landscape category with fox density estimated from faecal surveys along linear features (Webbon et al. 2004). Dashed line shows the fit of a linear model where the intercept is fixed at zero.
Estimated immigration rates of foxes were clearly related to the landscape type in which the estates were located. Immigration rates onto estates in arable a and b were four times greater than those in upland landscapes. The landscape-related pattern in immigration rates broadly matches landscape-related differences in fox density (estimated by faecal density counts; Webbon et al. 2004). Upland landscapes had lower estimates of fox density and immigration rate, while arable and pastural landscapes had higher estimates of both. However, within the arable and pastural categories there was no apparent relationship between immigration rate and density, as immigration rates onto estates in arable a, arable b and pastural b were higher than expected from fox density alone, and immigration rates onto arable c and pastural a were lower than expected. The reasons for this are not immediately clear, but are most likely to be due to an interaction between the factors determining fox density — prey availability and culling — and how each may influence immigration rate, as discussed below. We also cannot rule out some confounding between landscape type and estate area, as immigration rates were lowest in landscapes with the largest estates. Landscape types have a distinctly regional distribution (Supplementary material Appendix 1 Fig. A1.2) and to a large extent region and landscape are confounded, as shown by geographically closest estates tending to be in the same landscape (Fig. 1). Food resources are expected to vary among landscape types, and the observed variation in fox density is assumed to reflect this (Webbon et al. 2004); but culling practices also vary between regions and may limit fox density at a regional scale in parts of Britain (Tapper 1992, Heydon and Reynolds 2000a, b, Heydon et al. 2000).
The effect of gamekeeper effort also contributes to the larger cull per km2 on small estates, and is not independent of the increased immigration on such estates. Generally, the effort of a single gamekeeper is more concentrated on smaller estates thereby leading to more intensive culling. Some large estates employ several gamekeepers to ensure that a single gamekeeper's effort is not spread too thinly. In the subset of NGC data used in this analysis the number of gamekeepers employed does increase with estate area (Supplementary material Appendix 1 Fig. A1.3), but landscape type may be an influence, with arable estates employing more gamekeepers on an estate of a given area compared to upland estates. While information on the number of gamekeepers per estate is submitted to the NGC annually by contributing estates (Tapper 1992), it is unreliable as an indicator of predator control effort. On estates where game birds are hand-reared for release, gamekeepers generally have considerably less time for predator control compared with estates aiming to produce wild game, where predator control is a key ingredient (Tapper 1992, Reynolds and Tapper 1996). Game-bird releasing occurs predominantly in lowland arable landscapes and this probably explains the smaller area per gamekeeper in these landscapes (Supplementary material Appendix Fig. A1.3). Finally, motivation and skills are highly variable among individual gamekeepers, and may vary from year to year with personnel changes and personal life events; we used multi-year mean data in part to account for such variations in effort over time. In other respects, variation in fox control effort had to be regarded as unattributable noise around the relationship of interest.
Immigration rate is typically an unknown demographic parameter whose estimation requires expensive specialist study. In a context where tagging-based field methods are impractical, e.g. populations subject to intensive culling, our alternative approach offers an accessible starting point for further exploration. Our use of a Bayesian hierarchical modelling framework for analysis yielded plausible estimates of fox immigration rate for different landscape types. These estimates could be used for broad guidance on culling effort, as the rate of fox removal must exceed the immigration rate for control to have any effect on local density. However, our approach also provides informative results to be carried forward as a prior for population dynamics models in which immigration is a key process. The value of using suitable informative priors such as this has been repeatedly shown, whether they are derived from expert knowledge (Martin et al. 2005, Kuhnert et al. 2010), published data (McCarthy and Masters 2005, Martin et al. 2013), or other analytical methods (McAllister et al. 2001, 2010). In our own research (in prep.) we are using informative priors for immigration rate within population dynamics models to understand and improve the effectiveness of fox control at a local level (i.e. individual estates). Where suitable data on culling (or harvest) are available, this approach to quantifying immigration rate has potential to be useful in the management of a range of species.
Acknowledgements
The National Gamebag Census database is owned by the Game and Wildlife Conservation Trust (GWCT) and we are grateful for access to the data. We thank Nicholas Aebischer for advice on the analysis.
Funding — The study was funded by the GWCT.
References
Notes
[1] Supplementary material (available as Appendix wlb00416 at www.wildlifebiology.org/appendix/wlb-00416). Appendix 1–3.
