Population dynamics of two beaver species in Finland inferred from citizen-science census data J. E. BROMMER,1, R. ALAKOSKI,1 V. SELONEN,1 AND K. KAUHALA2 1Department of Biology, University of Turku, FI-20014 Turku, Finland 2Natural Resources Institute Finland Luke, It€ainen Pitk€akatu 3 A, FI-20520 Turku, Finland Citation: Brommer, J. E., R. Alakoski, V. Selonen, and K. Kauhala. 2017. Population dynamics of two beaver species in Finland inferred from citizen-science census data. Ecosphere 8(9):e01947. 10.1002/ecs2.1947 Abstract. A species’ distribution and abundance in both space and time play a pivotal role in ecology and wildlife management. Collection of such large-scale information typically requires engagement of vol- unteer citizens and tends to consist of non-repeated surveys made with a survey effort varying over space and time. We here used a hierarchical single-census open population N-mixture model, which was recently developed to handle such challenging census data, to describe the dynamics in the Finnish population sizes of the reintroduced native Eurasian beaver (Castor fiber) and the invasive North American beaver (Castor canadensis). The numbers of beaver winter lodges (i.e., family groups) were counted by volunteers in the municipalities of Finland every third year during 1995–2013. The dynamics of both species followed Gom- pertz logistic growth with immigration. Initial abundance of North American beavers increased with prox- imity to the introduction sites as well as with the amount of water in the municipality. The intensively hunted North American beaver population declined and the Eurasian beaver population increased during the study period. The model generated reasonable estimates of both total Finnish and local numbers of lodges, corrected for the incomplete detection. We conclude that the single-census N-mixture model approach has clear potential when using citizen-science data for understanding spatio-temporal dynamics of wild populations. Key words: beaver; citizen science; monitoring; population dynamics; population ecology. Received 27 January 2017; revised 3 April 2017; accepted 12 June 2017; final version received 17 August 2017. Corresponding Editor: Rebecca J. Rowe. Copyright: © 2017 Brommer et al. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.  E-mail: jon.brommer@utu.fi INTRODUCTION Understanding variation in distribution, abun- dance, and density in space and time is at the heart of ecology (Andrewartha and Birch 1954). A particular challenge in current ecology is to study abundances of organisms over spatial and tempo- ral scales, which are larger than what typically can be surveyed by researchers or research teams. Such large-scale spatio-temporal resolution is needed, for example, to understand how environ- mental change will alter the distribution patterns of species (Parmesan and Yohe 2003, Root et al. 2003) and what the community-level and ecosystem consequences of such changes could be (Sekercioglu et al. 2004). Further examples are the sustainable management of wildlife populations (Caughley and Sinclair 1994) and the control of invasive species (Williams et al. 2002), which require proper knowledge of distribution and abundance or density of species concerned. Addressing such large-scale ecological issues critically depends on data collected by large numbers of volunteers, who follow a simple pro- tocol, the so-called citizen-science data (Silver- town 2009). Citizen-science data present an unusual opportunity for ecologists to conduct research on hitherto challenging spatial and ❖ www.esajournals.org 1 September 2017 ❖ Volume 8(9) ❖ Article e01947 temporal scales (Dickinson et al. 2010, Bonney et al. 2014). Despite clear benefits, citizen-science data also present considerable analytical chal- lenges before proper inferences can be drawn. For example, spatio-temporal variation in survey effort partly drives changes in range margin of breeding birds in several countries when ranges are inferred directly from atlas data (Kujala et al. 2013). Because standard statistical models fall short, novel implementations or completely new procedures need to be developed for the analysis of citizen-science data (Dickinson et al. 2010, Hochachka et al. 2012). One avenue, for exam- ple, is to develop and explicitly take on board metrics of data quality (Isaac and Pocock 2015, Kelling et al. 2015a, b). In recent years, hierarchical models (sensu Ber- liner 1996) have become increasingly common in ecology and wildlife management. A powerful asset of ecological hierarchical models is that the state process is described separately from the observation process. The state process captures the spatio-temporal dynamics of interest, whereas the observation process links the changes in latent (i.e., unmeasured) population abundance to the census data (Royle and Dorazio 2008, Kery and Schaub 2012). By explicitly correcting observa- tions for the errors arising from the observation process, driven in particular by the incomplete detection inherent to count data, hierarchical models attain a clearer description of the ecologi- cal processes of interest, thereby avoiding com- mon biases occurring when using census data directly (e.g., Freckleton et al. 2006, Knape and de Valpine 2012). In the context of citizen-science data, variation in the propensity of volunteers to report their observations is likely to strongly influ- ence the observation process parameter “detection probability.” Nevertheless, provided the observa- tion process can be described, hierarchical models hence provide a fruitful platform for analysis of citizen-science data. From this perspective, a particularly exciting recent development is the extension of Royle’s (2004) N-mixture model to one-count-per-census- period data by Dail and Madsen (2011) and Hostetler and Chandler (2015). This “latest gener- ation” of open population N-mixture models presents a powerful approach for making good use of citizen-science data when modeling time- series data. This is because citizen-science data typically consist of a single count made during a survey period, whereas Royle’s (2004) original formulation of the model for open populations requires repeated counts to be made in each survey period (the so-called primary period) over two or more so-called secondary periods. The Dail and Madsen (2011) generalization of this approach requires time-series data, but a single census per survey period suffices. Statistical approaches for analysis of single-census data col- lected during one survey period have also been developed (e.g., Solymos et al. 2012). The single- census formulation of the open population N-mixture model presents an integration of large- scale population dynamics within a framework of explicit estimation of the latent population abun- dances. Such an integration has exciting applica- tion in the study of population ecology (Buckland et al. 2007), invasion ecology (Hooten et al. 2007), and range expansion (Pagel and Schurr 2011), as the approach opens up the research avenue of spatial population dynamics from census data, while correcting for bias caused by the observa- tion process. With this theoretical framework in place, it becomes an empirical question to evalu- ate the performance of these models to real-life data on wildlife populations. In this paper, we explore the use of the Hostetler and Chandler formulation for single-census open population N-mixture models for assessing spatial and temporal variation in abundance of the native Eurasian and invasive North American beaver populations in Finland. After being hunted to extinction from Finland in the latter half of the 19th century, the native Eurasian beaver (Castor fiber) was reintroduced in the country in the 1930s with individuals from Norway. However, also some North American beavers (Castor canadensis) were brought to Finland because their status as a distinct species was not known. Individuals of both species were released to several locations, but Eurasian beaver population survived and increased only in western Finland, whereas North American beavers flourished especially in eastern Finland. North American beaver population increased and spread much faster than the Eura- sian beaver population. Currently, the North American beaver occupies a larger range than the Eurasian beaver and has expanded its range close to the distribution area of the native beaver (Fig. 1). From a game management perspective, ❖ www.esajournals.org 2 September 2017 ❖ Volume 8(9) ❖ Article e01947 BROMMER ET AL. the target is to control North American beaver populations and to achieve a wider distribution and a sustainable level of native Eurasian beaver populations in Finland. To this end, information on the temporal dynamics and spatial variation in abundances of both species is an important aspect. The observations we use were obtained from a monitoring count of the number of active winter lodges of beavers in Finland. Unlike in summer, each beaver family uses only one lodge in winter. Hence, the number of active winter lodges pro- vides a reasonable estimate of the number of beaver families. This count is conducted in autumn every third year (1995–2013). The volun- teers carrying out the surveys are hunters who, in connection with hunting activity in autumn, report the number of active beaver winter lodges in the area covered by their hunting association. Hunters were used as volunteer surveyors because, firstly, Finland has approximately 300,000 hunters (6% of the population), and a sys- tem dispersing information and instructions to all hunters (Finnish Wildlife Agency). Second, most hunters belong to a hunting association. In total, there are about 4000 hunting associations in Fin- land, each covering an area ranging from 2000 to 10,000 ha. By collecting reports of beaver winter lodges from each hunting association, the scheme facilitates effective reporting and avoids multiple identifications of the same lodges. We here con- sider the administrative borders of municipalities as boundaries for sites, thereby grouping the reports of all hunting associations belonging to a municipality. The single-census open population N-mixture model (Dail and Madsen 2011, Hoste- tler and Chandler 2015) is an attractive approach for modeling counts of beaver lodges, because it is a single hierarchical model to assess population dynamics while accounting for the spatio- temporal uncertainty in the count data across municipalities. Clearly, not all winter lodges in a municipality are detected during a survey. Such uncertainty is inherent in the beaver census data as the level of reporting by volunteers differs both spatially (across municipalities) and across years. MATERIALS AND METHODS Study species The Eurasian beaver and the North American beaver are morphologically and ecologically simi- lar species. They are monogamous, herbivorous, crepuscular, and nocturnal animals that defend a territory by marking it with anal gland secretion and castoreum (Willson 1971, Nolet and Rosell 1994). Both species feed mainly on broadleaved trees, from which aspen (Populus tremula) and birches (Betula sp.) are the most favored, and in summer also aquatic plants are included in the diet (Lahti and Helminen 1974). In addition to deciduous trees nearby, beavers require an Fig. 1. Map of Finland with the municipalities occu- pied by Eurasian beavers (in blue) and by North American beavers (in red). The various shades of blue and red denote the game management districts within the national range of the respective beaver species. The part of Finland depicted in gray denotes the area where beavers are assumed to be absent as winter bea- ver lodges were never reported from these areas. ❖ www.esajournals.org 3 September 2017 ❖ Volume 8(9) ❖ Article e01947 BROMMER ET AL. aquatic habitat where they build a lodge (either free-standing lodge, bank burrow, or a combina- tion of these), and trees are also used for lodge and dam building (Willson 1971, Lahti and Helminen 1974, M€uller-Schwarze 2011). Beavers are known to build dams in order to regulate the water level at the site of the lodge, but lodges can also be in stream, river, and lake shores without dams. Beavers require a territory of approxi- mately 2 km of shoreline (Hartman 1994). Mostly two-year-old subadults disperse up to 3–13 km a year (in the Eurasian beaver) but often stay close to natal habitat if there is suitable habitat available (Hartman 1994, M€uller-Schwarze 2011). Also adult beavers may abandon the lodge and move to a new area when food resources become scarce (i.e., when they have consumed all suitable trees). Data on beavers Monitoring counts have been carried out every third year since 1995. Local hunters counted the number of active winter lodges in autumn (September–November), usually in connection with moose hunting, in the area of their hunting association. Food storage or other signs of activ- ity should be seen in the vicinity of an inhabited winter lodge. Hunters reported the number of winter lodges including the absence of winter lodges (i.e., count of zero). In addition, the type of winter lodge (bank burrow, bank lodge, or free-standing lodge) was recorded, as well as the estimated number of beavers, hunting bag of beavers, and the damage caused by beavers (flooded area or cut trees). Hunters returned the report to the Finnish Game and Fisheries Research Institute (nowadays part of the Natural Resources Institute Finland, Luke) via the Fin- nish Wildlife Agency. Winter lodges of the two beaver species cannot be distinguished. We considered winter lodges to belong either to Eurasian or North American beavers depending on the municipalities (Fig. 1). The rationale for this geographic division is based partly on historic information on the intro- duction and spread of the two beaver species (see Introduction), and partly on findings of the identities of beavers hunted in the municipalities. In particular, beaver skulls were collected from hunters for species determination, especially from the hunting districts where both species might occur (Fig. 1; eastern parts of Satakunta, western parts of Pohjois-H€ame, and southern parts of Pohjanmaa). DNA and morphometric analyses of the skulls confirmed species identity. In other areas, the identity of the species was presumed according to the history of beavers in Finland. Thus, for the period considered in this study, we consider the municipalities of Eurasian and North American beavers as distinct sets of local populations. Model We fitted the data to single-census open popu- lation N-mixture models, which were originally formulated by Dail and Madsen (2011) and extended by Hostetler and Chandler (2015). We here briefly summarize the conceptual aspects of this model. Details of the full model are provided in Appendix S1. The open population N-mixture models assumed that several local populations (sites) were censused at some time interval. The initial abundance N at time t = 1 at site i was assumed to be Ni1 f ðh1Þ (1) where h1 is a vector of one or multiple parame- ters describing a distribution, either a negative binomial (h1 includes both mean initial abun- dance Λ and dispersion parameter a), a Poisson (h1 includes Λ), or a zero-inflated Poisson distri- bution (h1 includes both mean initial abundance Λ and the proportion of local populations with abundance of zero). Covariates x on initial abun- dance Λ for site i were included using a logarith- mic link, where Ki ¼ exp K0 þ X c Kcxci ! . (2) The dynamics in the consecutive n time steps are then given by the recursive Poisson stochastic NitPðNiðt1Þf ðhtÞ þ iÞ (3) where ht is a vector of one or more (possibly time-varying) parameters specified by a certain population-dynamical equation and ι is the expected number of immigrants in a local popu- lation. In the Dail and Madsen (2011) model, f(ht) includes separate estimators of survival and reproduction. In the extension of this model by Hostetler and Chandler (2015), the generalization to include Ricker logistic (ht includes both ❖ www.esajournals.org 4 September 2017 ❖ Volume 8(9) ❖ Article e01947 BROMMER ET AL. maximal potential growth rate r and carrying capacity K; Ricker 1954) NitP Ni;t1 exp r 1Ni;t1K    þ i   (4) and Gompertz logistic dynamics (ht includes both maximal potential growth rate r and carry- ing capacity K; Gompertz 1825), NitP  Ni;t1  exp  r  1 logðNi;t1Þ þ 1 logðK þ 1Þ  þ i  (5) as well as immigration ι. Covariates x for growth rate r at site i at time t were modeled using a logarithmic link as rit ¼ exp r0 þ X c rcxitc ! . (6) Under the assumption of a specific mixture model for initial abundance and consecutive dynamics, this state process describes the popu- lation dynamics of local population i over time t as latent abundance states Nit. These latent states are then compared to the observed data yit, assuming that yitB Nit; pitð Þ (7) where pit is the probability to detect a beaver lodge at time t in municipality i. The detection probability was hence inferred directly from the comparison of the time series in latent-state abundances N (which are assumed to be Marko- vian) to observed data y (Dail and Madsen 2011). The effect of the number of reports on detection probability for site i at time t was modeled assuming a logistic link as pit ¼ expit p0 þ preportsReportsit þ eit   ¼ exp p0 þ preportsReportsit þ eit   1þ exp p0 þ preportsReportsit þ eit   (8) where “expit” denotes the inverse logit, which is further written out explicitly. An overview of all parameters and covariates in the model is provided in Table 1. Below, we detail the covariates used. Municipalities and their landscape composition Each site i in the model outlined above was a municipality in Finland where beavers were counted (Fig. 1). Municipalities were defined according to borders in 2013 (National Land Sur- vey of Finland 2/2015), and some older monitor- ing counts in originally separate municipalities were thus combined according to what the municipalities’ borders were in 2013. Landscape variables were used as covariate to explain abun- dance at the first survey year (Eq. 1, i.e., as covariate for beaver abundance in 1995). Land- scape variables were computed in ArcMap 10.2. (ESRI 2011). Given the biology of beavers, we assumed that important metrics are those describing the length of watercourses, water areas, and mixed and deciduous forest (cf. Kau- hala and Turkia 2013). The total length of water- courses (streams with width ≤20 m, in units of 10 km) and the total area of water areas (lakes Table 1. Synopsis of the model terms, following the nota- tion of Hostetler and Chandler (2015), and the covari- ates used in the information-theoretical modeling. Parameter Explanation Λ Mean initial abundance w Proportion of sites potentially occupied by the species a Dispersion parameter of negative binomial distribution r Population growth rate in Ricker or Gompertz model K Carrying capacity in Ricker or Gompertz model p Detection probability ι Immigration, modeled as a number added each time step to each site c Recruitment probability in the Dail and Madsen (2011) formulation x Survival probability in the Dail and Madsen (2011) formulation rep Number of reports on winter beaver lodges reported in the municipality prox Proximity (Eq. 9) of the municipality to the species’ introduction site(s) str Area covered by watercourses in the municipality lak Area covered by lakes and larger water areas in the municipality for Area covered by deciduous and mixed forest in the municipality hunting Hunting pressure in the district to which the municipality belongs temp Geometric mean temperature in the three winters between censuses in the municipality snow Geometric mean snow cover in the three winters between censuses in the municipality ❖ www.esajournals.org 5 September 2017 ❖ Volume 8(9) ❖ Article e01947 BROMMER ET AL. and large rivers, in units of ha) in a municipality were computed using the large-scale topographi- cal map (1:100,000) in the open digital data of the National Land Survey of Finland (2/2015). In general, the distribution area of the North Ameri- can beaver covers the areas in Finland where large lakes are abundant, whereas in the distribu- tion area of the Eurasian beaver large lakes are scarce and aquatic habitats consist mostly of riv- ers, streams, and even small ditches between fields. We assumed that, from the perspective of beavers, the amount of the above-outlined aqua- tic habitats represents the size of each municipal- ity and therefore we did not include the total area size of a municipality as a covariate. The pooled percentage (%) of deciduous and mixed forests covering the landscape in a munic- ipality was computed using the Corine land cover maps of Finland in 2000 released by Fin- nish Environment Institute (SYKE 2/2015). In general, coniferous or mixed forests cover most of Finland. Deciduous forests occur abundantly only in the southwestern tip of Finland. Half (50%) of forest tree species in Finland consist of scots pine (Pinus sylvestris), 30% of Norway spruce (Picea abies), 17% of birch (Betula pubescens and Betula pendula), and 3% of other broadleaved species. Other broadleaved tree species include Eurasian aspen (P. tremula), alder (Alnus incana and Alnus glutinosa), European mountain ash or rowan (Sorbus aucuparia), and goat willow (Salix caprea; Ylitalo 2013). Hunting of beavers North American beavers can be hunted freely during the hunting season from 20 August to 30 April, although before 2001 a special license was demanded. Catches of North American beavers are reported voluntarily and their numbers are hence indicative. A license is demanded to hunt the Eurasian beaver and because there is a quota for this species, the reporting of the number of individuals hunted is accurate. In general, illegal poaching may cause substantial differences between reported and actual numbers killed. However, in the case of the North American bea- ver, hunting was not restricted during most of the years studied. For the Eurasian beaver in Fin- land, we believe that the official hunting statistics present a reasonable estimate. Officials received each year fewer applications for permits than annual hunting quota for Eurasian beaver, and obtaining permits was hence not highly restric- tive for this species. The official hunting statistics are gathered per hunting district (see Fig. 1) and are reported by the Natural Resources Institute Finland Luke (http://stat.luke.fi/metsastys). We hypothesized that the intensity of hunting would reduce population growth rate. We there- fore calculated a measure of hunting pressure (number of hunted beavers divided by the total area covered by all municipalities where beavers were counted in the hunting district). The hunt- ing pressure for all municipalities belonging to a hunting district was assumed to be equal. Other covariates We hypothesized that the initial abundance of beaver lodges in municipality i could be affected by the proximity of this municipality to the site(s) where the Eurasian and North American beavers were introduced. We hence calculated a weighted mean proximity based on the distance d of municipality i to introduction site j out of the k introduction sites where we exponentially dis- counted longer distances as Proxi ¼ Xk j¼1 expðdijÞ (9) where dij was the Euclidean distance (in units of 100 km). For the Eurasian beaver, there was only one introduction site (Noormarkku, Satakunta) where beavers survived and proximity of munici- palities in which the Eurasian beaver occurred was hence with respect to this site only. For the North American beaver, there were four introduc- tion sites and the proximity considers the spatial distances of each municipality to all these sites. We hypothesized that climatic variation could affect the growth rate (r) of the population. We extracted weather information from the Finnish Meteorological Institute (FMI). We considered all weather stations situated up to 200 km from the outer municipalities where North American and Eurasian beavers occurred, respectively. For each of these weather stations, the arithmetic average temperature and snow cover recorded in the winter months (December–February) were com- puted. For each three-year census period, we then calculated the geometric mean of three win- ters relevant to the time period. For example, to ❖ www.esajournals.org 6 September 2017 ❖ Volume 8(9) ❖ Article e01947 BROMMER ET AL. explain the dynamics between 1995 and 1998, the winters of 1995/1996, 1996/1997, and 1997/1998 were included. In order to arrive at a measure of winter weather for each municipality, we used ordinary kriging (Cressie 1993) over the geomet- ric mean winter weather of all weather stations with data for a particular census period. Ordi- nary kriging assumes a spatial autocorrelation in the data, where the correlation between two sites reduces following an exponential function of the distance between two sites. Based on the kriging parameters, geometric mean winter weather was predicted for the center coordinate of each municipality. The ordinary kriging and predic- tions were performed using the R package “gstat” (Pebesma 2004) in R (R Core Team 2015). Prior to analysis, we mean standardized prox- imity, the amount of watercourses and water areas, the amount of deciduous and mixed forest, temperature, and snow depth in order to con- sider deviation (both negative and positive) from the zero-mean intercept due to these covariates. Model selection and Bayesian model construction We followed the approach by Hostetler and Chandler (2015) and first selected the most parsi- monious model using Akaike’s information crite- rion (AIC)-based model selection where models were solved using a likelihood-based approach. These models were implemented in a likelihood framework in the package “unmarked” (Fiske and Chandler 2011) in R (R Core Team 2015). Hostetler and Chandler (2015) used simulations to show that AIC-based model selection was effi- cient in identifying the correct candidate model. As in Hostetler and Chandler (2015), after model selection, we constructed a Bayesian version of the most parsimonious model. The Bayesian ver- sion of the model was implemented in JAGS (Plummer 2003). Transferring to the Bayesian framework is done here because (1) it allows explicit calculations to be made with the latent- state variables and hence estimation of total population sizes for either the whole or part of Finland, (2) it allows a flexible addition of ran- dom effects to the model in an effort to improve model fit. A general approach for improving fit of a Bayesian model is to add random effects in order to explicitly model variance in model parameters, which would otherwise be assumed to be determined completely by the intercept and possible fixed effects included (e.g., Kery and Schaub 2012). Hostetler and Chandler (2015) added random effects to describe inter-annual variation in population growth rate and detec- tion probability, but did not specifically assess the fit of their models. Fit of the model to the data was assessed using Bayesian posterior predictive checks (PPC; Gel- man et al. 2004). For each iteration of posteriors, we calculated for each municipality i in which beaver lodges were counted at time t the Free- man–Tukey discrepancy as a statistic of fit com- paring the observed number of beaver lodges with the expected number of beaver lodges as Dobst ¼ Xk i ffiffiffiffiffi yit p  ffiffiffiffiffiffiffiffiffiffiffiNitpitp 2. (10) In addition, for each iteration of posteriors, the Freeman–Tukey discrepancy measure Drept was calculated following the same logic as Eq. 6 except using a replicate model-predicted observed number of beaver lodges and the prediction. Again, Drept was calculated for each municipality i where beaver lodges were counted at time t. That is, discrepancy was only calculated when there were observed data. The posterior draws for the replicate model-predicted census data Drept pro- vided the reference distribution under the null hypothesis that the model fits the data. A so- called Bayesian P-value describing model fit can then be computed as the proportion of times that the posterior of discrepancy in the observed data Dobst has a more extreme value than D rep t . The Bayesian P-value should be close to 50% under the null hypothesis that the model parameters predict the number of beaver lodges. While extreme Bayesian P-values (close to 0 or 1) sug- gest a poor model fit, it should nevertheless be noted that the procedure does not provide an a priori-defined rigorous threshold value for signifi- cant deterioration in model fit and is hence a descriptive technique. We performed a PPC for the AIC-based most parsimonious model, and, following Hostetler and Chandler (2015), for a model with (1) between-year residual variation in population growth rate (r), (2) between-year residual varia- tion in population growth rate (r) and detection probability (p). We then (3) further extended residual variation in detection probability to ❖ www.esajournals.org 7 September 2017 ❖ Volume 8(9) ❖ Article e01947 BROMMER ET AL. include both variation over municipalities and years. For comparison, we also ran models with random effects on other parameters, but this did not qualitatively change the results of model fit and these results are not reported in this paper. As a derived parameter, we calculated the sum of the abundance of lodges over all municipali- ties in each hunting district and for the whole of Finland for each census year, which hence esti- mated the total population of the respective bea- ver species in hunting districts and nation-wide, respectively. The Bayesian model (Text S1) was imple- mented in JAGS (Plummer 2003), based on three chains. For the data on the Eurasian beaver, we ran 50,000 iterations after a burn-in of 5000 itera- tions, sampling every 50th value. For the data on the North American beaver, after a burn-in of 10,000 iterations, an additional 300,000 model iterations were run where every 300th iteration was saved as posteriors. Mixing of the chains was assessed by eye and Gelman–Rubin statistic (Gelman et al. 2004). The details of the model construction, prior settings, and evaluation are presented in Appendix S1. RESULTS Description of the data used There were 168 and 29 municipalities where the North American and Eurasian beaver occurred, respectively (Fig. 1, Table 1). For the Eurasian beaver, the reporting intensity was higher than for the North American beaver (Table 1; 80% vs 65%, v2 = 15.6, df = 1, P < 0.001). When one or more reports were filed for a census occasion, reporting was reasonably active with on average approximately eight and four reports for North American and Eurasian beavers, respectively. Typically, 5–6 beaver winter lodges were reported per census occasion, although the variation in the number of lodges was considerable (Table 1). In 18% (213/1176) and 14% (30/203) of all census occasions, North American and Eurasian beavers, respectively, were reported to be absent. In 24% (40/168) and 10.3% (3/29) of the municipalities, no North American and Eurasian beavers, respectively, were recorded during the whole study period. In terms of the covariates used for the various model parameters, differences between beaver species reflected differences in landscape struc- ture and climate in the parts of Finland occupied by the species (Fig. 1, Table 2). The most note- worthy difference is that the hunting pressure on North American beavers tended to be larger than on the Eurasian beaver (Table 2). Model selection In general, model comparison based on AIC showed clear difference in AIC values of the dif- ferent candidate models for the dynamics of North American and Eurasian beavers (Appendix S1: Table S1). The negative binomial distribution best explained initial abundance in 1995 in both beaver species (section A in Appendix S1: Table S1). The number of reports filled in a municipality was an important covari- ate for detection in both species (section B in Appendix S1: Table S1). For both species, models including immigration had lower AIC and thereby fitted the data better than models with- out immigration. The Gompertz logistic model with immigration provided the most parsimo- nious description of the census data of both North American and Eurasian beavers (section C in Appendix S1: Table S1). The amount of water- courses and water areas in a municipality and its proximity to the sites where the species was introduced positively affected the initial abun- dance of North American beavers, and, in addi- tion to these variables, the amount of deciduous and mixed forest positively affected initial abun- dance in Eurasian beavers (section D in Appendix S1: Table S1). Warmer winter tempera- ture positively affected population growth rate in the North American beaver, but did not have consequences for the population growth in Eura- sian beavers (section E in Appendix S1: Table S1). The dynamics of North American and Eurasian beavers in Finland We constructed a Bayesian version of the most parsimonious models for North American and Eurasian beavers (details in Appendix S1). Poste- rior predictive check of the most parsimonious model for the North American and Eurasian bea- vers revealed a poor fit to the data (pPPC = 0.01). Following Hostetler and Chandler (2015), we included additional random effects for between- year variance in population growth rate (r, pPPC = 0.01), followed by addition of between- ❖ www.esajournals.org 8 September 2017 ❖ Volume 8(9) ❖ Article e01947 BROMMER ET AL. year variation in detection (p, pPPC = 0.01). The model with between-year variation in population growth rate and between-year and between-site variance in detection probability showed reason- able fit to the data for the Eurasian beaver (pPPC = 0.12). For the North American beaver, however, this model still fitted the data poorly (pPPC = 0.01) and only when residual variance in detection probability was allowed to be heteroge- neous (i.e., year specific) did the model show rea- sonable fit to the data (pPPC = 0.22). Adding uncertainty (random effects) to other parameters did not increase the model fit for either species. Parameter estimates of the final Bayesian mod- els showed that population growth rate in North American beavers was very low (transformed r0 was close to zero, Table 3). The observed decline of North American beavers was largely caused by a low carrying capacity (approximately six winter beaver lodges per municipality, Table 3), which was surprisingly low especially when compared with mean initial abundance (Λ was approximately 36, Table 3). The effect of covari- ates when based on the Bayesian model showed that initial (i.e., 1995) abundance of North Ameri- can beavers was higher in municipalities with more bodies of water (streams and lakes), which were closer to the sites of introduction (Table 3). However, the positive effect of temperature on the population growth rate, which the AIC-based model selection suggested to be important (Appendix S1: Table S1), was in this final model not statistically significant (Table 3). Again, an important aspect of the Bayesian model for fitting reasonably well to the data was that heterogeneity in detection probability p across municipalities was larger in the last three survey years compared to the first five survey years (Table 3). After 2001, hunters did not need a hunting license for North American beavers any- more. Most likely, this increased heterogeneity (and reduced average level of reporting) was associated with this change in legislation. The parameters of the Bayesian model for the dynamics of the Eurasian beaver revealed a strik- ingly different pattern compared to the North American beaver. Despite having about equal initial mean abundance as the North American beaver (approximately 32 winter lodges per municipality, Table 3), Eurasian beavers had a clearly positive population growth rate and high carrying capacity (about 214 lodges per munici- pality, Table 3). In the Bayesian implementation of the model, AIC-selected covariates for initial Table 2. Descriptive statistics of the information used in the model. Characteristic North American Eurasian Municipalities 168 29 Total no. census occasions 1176 203 No. census occasions with ≥1 report filed (%) 768 (65) 162 (80) No. reports/census occasion with ≥1 report filed 8.1 (1, 79) 4.1 (1, 14) No. lodges/census occasion with ≥1 report filed 4.6 (0, 400) 5.6 (0, 164) Municipality-specific covariates (95% CRI) Proximity to introduction sites 3.21 (2.73, 3.45) 0.95 (0.92, 0.99) Streams (km) 2950 (243, 11046) 1540 (367, 4446) Lakes (ha) 15515 (246, 64076) 2650 (115, 10315) Deciduous and mixed forest (%) 22.8 (13.4, 35.5) 16.3 (7.2, 24.9) Municipality- and year-specific covariates (95% CRI) Hunting pressure (beavers per 1000 ha per 3 yr) 0.37 (0, 1.80) 0.020 (0, 0.055) Snow depth (cm) 27.3 (17.7, 33.9) 15.4 (8.1, 21.9) Temperature (°C) –7.7 (–10.0, –5.0) –4.7 (–7.4, –2.5) Notes: Presented are the mean with, when relevant, the minimal and maximal value, and 95% credible interval (CRI) between brackets. Each municipality has a census occasion every third year during the study period consisting of seven censuses in total (1995–2013). The total number of (no.) census occasions was hence seven times the number of municipalities censused. We list in how many of these census occasions at least one report was made, and statistics for the number of reports, as well as how many beaver lodges were reported, per census occasion with at least one report made. Municipality-specific covariates were used to model initial abundance and are specific to the start of the study period (1995). Municipality and year- specific covariates varied both over space and time. Hunting pressure is the number of beaver killed per 1000 ha in the game management district in the three-year time interval between the censuses. Temperature and snow depth are averages for December–February in the three-year period between the censuses in the municipalities studied. ❖ www.esajournals.org 9 September 2017 ❖ Volume 8(9) ❖ Article e01947 BROMMER ET AL. abundance were not significantly different from zero (Table 3). Detection probability of beaver lodges was low (approximately 20% and 35% for North American and Eurasian beaver, respectively, at the average level of reporting as indicated by dashed vertical lines in Fig. 2), although it rapidly increased when additional reports were made in the municipality, especially for the Eura- sian beaver (Fig. 2, Table 3). The total abundance of beaver lodges in Finland varied during the study period (Fig. 3). Qualitatively, the pattern of estimated number of lodges corresponded rea- sonably well with the total numbers counted (black line in Fig. 3) and a simple index of abun- dance, the number of lodges divided by the Table 3. Parameters estimated by the Bayesian version of the single-census open population N-mixture model by Hostetler and Chandler (2015) for the number of lodges of North American and Eurasian beaver in 169 munici- palities in Finland 1995–2013 censused every third year. Parameter Scale Transformed scale Data scale Mean 2.5% 97.5% Mean 2.5% 97.5% North American beaver a 0.69 0.50 0.91 Λ exp 3.6 3.2 3.9 35.9 25.5 50.6 K exp 1.8 1.4 2.9 6.2 0.24 17.3 r0 exp 3.9 18.2 2.2 0.02 108 0.11 ι 0.16 0.02 0.39 p0 logit 2.2 2.5 1.8 0.10 0.08 0.14 Covariate (parameter) Prox. to intro (Λ) exp 3.6 2.3 5.0 Streams (Λ) exp 0.02 0.01 0.03 Lakes (Λ) exp 0.02 0.01 0.03 Temperature (r) exp 0.61 1.0 2.0 Reports (p) logit 0.10 0.08 0.12 Variances Var in r over years 0.07 0.01 0.40 Var in p 1995–2004 0.37 0.21 0.57 Var in p 2007–2013 0.69 0.23 1.1 Eurasian beaver a 0.43 0.20 0.83 Λ exp 3.4 2.6 4.3 32.2 12.9 71.8 K 214 129 360 r0 exp 0.72 0.08 1.63 2.1 1.1 5.1 ι 0.33 0.08 0.75 p0 logit 1.7 2.4 1.0 0.15 0.09 0.27 Covariate (parameter) Prox. to intro (Λ) exp 15.7 4.2 34.2 Streams (Λ) exp 0.04 0.04 0.14 Lakes (Λ) exp 0.14 0.45 0.25 Forest (Λ) exp 0.14 0.06 0.35 Reports (p) logit 0.25 0.14 0.39 Variances Var in r over years 1.2 0.04 3.5 Var in p 0.78 0.39 1.4 Notes: Akaike’s information criterion-based model selection (Appendix S1: Table S1) led to the modeling of initial abundance following a negative binomial distribution with dispersion parameter a and mean initial abundance Λ. Mean initial abundance Λ, growth rate r, and detection probability p were modeled as being dependent on a number of covariates (see Appendix S1 for further explanation). Consecutive dynamics assumed Gompertz logistic growth with growth rate r0 in the absence of the effect of covariates, and carrying capacity K, as well as immigration (ι). Estimated abundances were assumed to be observed with binomial probability p0 for which the covariate was the number of observation reports filled by local volunteers in each munici- pality. Parameter values are noted on the scale in which they were estimated in the model, and parameters for the intercept values are also back-transformed to the data scale. Parameter estimates for covariates whose 95% credible interval did not include 0 are considered as statistically significant and appear in boldface. ❖ www.esajournals.org 10 September 2017 ❖ Volume 8(9) ❖ Article e01947 BROMMER ET AL. number of reports (red line in Fig. 3). Numeri- cally, nevertheless, the model allowed compar- ison of abundance changes within species across years and also between species by incorporating the difference in reporting (Fig. 3c, d). The num- ber of winter lodges of North American beavers declined during the study period (geometric mean population growth rate was 0.95, 95% credible interval [CRI]: 0.93, 0.97), but increased for Eurasian beavers (1.11, 95% CRI: 1.06, 1.17). DISCUSSION We here describe the dynamics of the number of family groups (winter lodges) of Eurasian and North American beavers in Finland using a recently developed single-census open population N-mixture model (Hostetler and Chandler 2015). An exciting feature of this model approach is that it estimates abundance by correcting for the detec- tion probability of the species censused. Thus, the model produces abundance estimates over both time and space, which can be used for numerical comparisons. The modeling results underline the similarity of the two beaver species studied, with respect to their population dynamics. Dynamics in both species is density dependent following Gompertz (1825) formulation, and immigration is an important aspect. Another striking feature is that in the North American beaver the initial abundance (i.e., abundance in 1995) is higher in municipalities, which are closer to the original introduction site(s). Hence, a spatial signature of the introduction is still observable six decades after the beaver species were introduced. This might indicate limits in dispersal capacity or that there are still suitable habitats near introduction sites, and therefore, distribution ranges have not expanded more (Nolet and Rosell 1994). Our findings demonstrate that Eurasian bea- vers are numerically increasing in Finland, while North American beavers, although currently more widespread, are declining in numbers. The primary reason behind this difference in dynam- ics appears to be the surprisingly low carrying capacity (K) in North American beaver lodges in combination with a low intrinsic rate of popula- tion increase (r) in the Gompertz dynamics for this species. It is currently not clear what the rea- sons are for the low values of these estimated parameters. In particular, North American beavers have a larger litter size than Eurasian beavers. However, it is likely that these parame- ters are strongly affected by the hunting on this species. Although our models suggest that inclu- sion of hunting pressure is not very important, it should be noted that our estimates of hunted North American beavers are coarse, both spa- tially (hunting statistics on the level of munici- pality are not available), and in terms of accuracy because hunting statistics are built by sending an inquiry to only a fraction (about 2%) of hunters. It is difficult to evaluate the effect of habitat parameters on North American beaver numbers, because hunting likely affects key population parameters of the North American beaver Fig. 2. The probability to detect a beaver winter lodge during a census of a municipality as a function of how many reports were received from volunteers. The posterior mean (line) and its 95% credible interval are plotted over the range of reports filed. The dashed vertical line indicates the mean number of reports filed per municipality. ❖ www.esajournals.org 11 September 2017 ❖ Volume 8(9) ❖ Article e01947 BROMMER ET AL. population and thus their local abundance in dif- ferent habitats in a manner which we cannot fully address within our model. Still, our find- ings in this respect are intuitive; North American beavers show a high mean initial abundance in municipalities with more water areas and water- courses. A comprehensive water system also enables beavers to move more broadly in search of food. Winter temperature, which AIC-selected modeling suggested to be a covariate affecting intrinsic population growth rate in North Ameri- can beavers, was not significant in the Bayesian implementation of the model. Likewise, all AIC- selected covariates for initial abundance for the Eurasian beaver were not significant in the final Bayesian implementation of the model. In gen- eral, the addition of random effects increases parameter uncertainty in population ecological models (Kery and Schaub 2012). Hostetler and Chandler (2015) also noted that estimates change when random effects were added in the Bayesian implementation of their model. Importantly, however, the addition of random effects did pro- duce models with a reasonable fit to the data as assessed using PPC, whereas the AIC-selected models fitted the data very poorly. Thus, we find that evaluation of fit of the single-census open population N-mixture model followed by— when needed—inclusion of random effect(s) is a crucial aspect of model building before drawing model inferences. Nevertheless, our list of covari- ates clearly does not cover all aspects relevant to beavers. For example, aquatic vegetation (Law et al. 2014) and agricultural plants may affect the Eurasian beaver population sizes also at large scale (instead of the covariates used here). The inclusion of spatio-temporal heterogeneity in detection probability is an essential aspect of the use of count data, because drawing inferences on “raw” count data implicitly makes the strong assumption that numbers counted are equal or at least proportional to true abundance (Royle 2004). Fig. 3. Population dynamics of the number of lodges of (a) North American and (b) Eurasian beavers in Fin- land during 1995–2013 based on monitoring every third year. For each species, the model prediction (blue line with 95% credible interval) and sum of all observed lodges (black line) are shown as well as the mean number of lodges counted per report (red line; scale on right-hand side). The lower panels show the average number of reports per municipality of (c) North American and (d) Eurasian beaver lodges filed during the study period. ❖ www.esajournals.org 12 September 2017 ❖ Volume 8(9) ❖ Article e01947 BROMMER ET AL. For this reason, a dimensionless “abundance index” is often constructed and temporal dynam- ics of species is inferred on the basis of this index (Sutherland 1996). However, a particular chal- lenge then arises when the interest is in compar- ing count data across species, as, for example, in our case, where one would want to get a handle on the actual abundance of the species. Doing so requires information on the detection probability (Kery and Schaub 2012). Our use of the single- census open population N-mixture model shows that even closely related and ecologically similar species may showmarked differences in detection probability and especially in how the probability of detection depends on reporting intensity. To achieve an 80% detection probability for Eurasian beaver winter lodges, approximately 13 reports are needed per municipality, whereas for North American beavers the required reporting intensity was approximately 40 reports per municipality. With greater reporting intensity, the probability to detect Eurasian beaver winter lodges increases faster than the probability to detect North Ameri- can beaver winter lodges. This difference in how detection changes with reporting intensity likely stems from the relatively higher densities in the municipality of Eurasian compared to North American beavers. For both species, nevertheless, the level of reporting required to reach 80% prob- ability of detection is markedly above the typical reporting intensity. For the North American bea- ver, it is furthermore apparent that the willingness to report beaver lodges decreased when a hunting license for this species was not demanded any more. Reasonable model fit for this species is only achieved by allowing considerable heterogeneity in detection probability, both across sites and across years. Thus, in order to fit the data, the model requires that each census in a given munic- ipality at a specific point in time is viewed as an event with an own “unique” probability of detec- tion. Such heterogeneity likely partly reflects the fact that there are unknown covariates for the observation process. Indeed, the “number of reports” is only a coarse description of factors affecting detection probability, primarily because it does not contain information on the survey effort behind a report, which conceivably varies both in time and space. Notably, however, hetero- geneity in detection probability may simply be an inherent aspect of census counts, especially when conducted by citizens. It is in general considered challenging to obtain reasonable model fit for open population N-binomial mixture models (Kery and Schaub 2012). Despite the challenges in fitting the single-census open population N-binomial mixture model, it is also evident that this approach facilitates comparison of large-scale population dynamics to be made within species over time and between species in terms of biologi- cally relevant properties. Data collection by volunteers presents an excit- ing opportunity for ecologists to work on large spatial and temporal scales (Bonney et al. 2014). At the same time, such data are inherently chal- lenging to analyze, and require harnessing new statistical modeling approaches (Hochachka et al. 2012, Isaac and Pocock 2015). Voluntary reporting of single species, as used for the beaver winter lodge data analyzed in this paper, is a census tech- nique likely to result in underreporting (Suther- land 1996). There is no non-arbitrary solution to the issue of high heterogeneity in voluntary report- ing other than to convince volunteers to increase their reporting intensity in all regions, also when beavers are assumed to be rare or absent. Technical advances, which lower the threshold of reporting, have been heralded as an important develop- ment in citizen science (Bonney et al. 2014). Furthermore, citizen-science monitoring schemes designed to census multiple species at the same time, based on checklists, automatically include a more systematic reporting of absences and low counts for each species compared to voluntary censuses of specific species such as the Finnish beaver monitoring scheme. Nevertheless, our findings show that the single-census open popula- tion N-mixture model by Hostetler and Chandler (2015) can be a good tool to obtain population ecological insight using large-scale citizen-science single-species monitoring data. Although the single-census-per-survey scheme in our case produced reasonable inferences on abundance and dynamics, we caution against interpreting this survey approach as a good over- all solution. Large-scale monitoring of any organism using volunteer efforts is clearly cost- effective, but comes at the price of collecting data with high heterogeneity in detection probability. Any approach that collects information to explic- itly estimate detection probability from the data (e.g., repeated censuses, distance sampling) will ❖ www.esajournals.org 13 September 2017 ❖ Volume 8(9) ❖ Article e01947 BROMMER ET AL. allow stronger inferences (Royle and Dorazio 2008), as illustrated in our findings by the need to include high uncertainty in detection. Neverthe- less, the Hostetler and Chandler (2015) modeling approach opens up exciting research avenues in terms of, for example, understanding invasive species dynamics or carrying out comparative studies across species of key population dynami- cal parameters, on the basis of relatively “cheap but messy” survey data. In general, further spatial aspects could include conditioning of the growth rate or the immigration rate on the spatial location of the site, for example, through auto-regressive arguments (Banerjee et al. 2004) where the popu- lation abundance in neighboring sites plays a role in determining the dynamics of a specific site. Furthermore, the distribution of the two beaver species in Finland will likely become increasingly overlapping, and the use of spatial information to infer species identity will likely become a valuable tool (cf. Conn et al. 2015), as well as inclusion of the details of the landscape in terms of water sys- tems (dispersal corridors) and watershed areas (dispersal barriers). ACKNOWLEDGMENTS We thank all the volunteers who have over the years participated in counting beaver lodges. We are grateful to A. Ermala for organizing the monitoring counts in 1995–2007 and R. Koivunen for the help in filing the data. This study was supported by Maj and Tor Nes- sling Foundation (to RA) and the Academy of Finland (to VS, to JEB). We thank reviewers for insightful com- ments, which improved the message of this paper. LITERATURE CITED Andrewartha, H. G., and L. C. Birch. 1954. Distribu- tion and abundance of animals. University of Chicago Press, Chicago, Illinois, USA. Banerjee, S., B. P. Carlin, and A. E. Gelfand. 2004. Hier- archical modeling and analysis for spatial data. Chapman and Hall and CRC Press, Boca Raton, Florida, USA. Berliner, L. M. 1996. Hierarchical Bayesian time series models. Pages 15–22 in K. M. Hanson and R. N. Sil- ver, editors. Maximum entropy and Bayesian methods. Springer, New York, New York, USA. Bonney, R., J. L. Shirk, T. B. Phillips, A. Wiggins, H. L. Ballard, A. J. Miller-Rushing, and J. K. Parrish. 2014. Next steps for citizen science. Science 343: 1436–1437. Buckland, S. T., K. B. Newman, C. Fernandez, L. Tho- mas, and J. Harwood. 2007. Embedding population dynamics models in inference. Statistical Science 22:44–58. Caughley, G., and A. R. E. Sinclair. 1994. Wildlife ecol- ogy and management. Volume 334. Blackwell Science, Cambridge, UK. Conn, P. B., D. S. Johnson, J. M. Ver Hoef, M. B. Hoo- ten, J. M. London, and P. L. Boveng. 2015. Using spatiotemporal statistical models to estimate animal abundance and infer ecological dynamics from survey counts. Ecological Applications 85 :235–252. Cressie, N. A. C. 1993. Statistics for spatial data. Wiley Blackwell, London, UK. Dail, D., and L. Madsen. 2011. Models for estimating abundance from repeated counts of an open metapopulation. Biometrics 67:577–587. Dickinson, J. L., B. Zuckerberg, and D. N. Bonter. 2010. Citizen science as an ecological research tool: chal- lenges and benefits. Annual Review of Ecology, Evolution and Systematics 41:149–172. ESRI. 2011. ArcGIS desktop. Release 10. Environmen- tal Systems Research Institute, Redlands, Califor- nia, USA. Fiske, I., and R. Chandler. 2011. Unmarked: an R pack- age for fitting hierarchical models of wildlife occur- rence and abundance. Journal of Statistical Software 43:1–23. Freckleton, R. P., A. R. Watkinson, R. E. Green, and B. J. Sutherland. 2006. Census error and the detec- tion of density dependence. Journal of Animal Ecology 75:837–851. Gelman, A., J. B. Carlin, H. S. Stern, and D. B. Rubin. 2004. Bayesian data analysis. Second edition. Chapman and Hall, Boca Raton, Florida, USA. Gompertz, B. 1825. On the nature of the function expressive of the law of human mortality, and on a new mode of determining the value of life contin- gencies. Philosophical Transactions of the Royal Society of London 115:513–583. Hartman, G. 1994. Ecological studies of a reintroduced beaver (Castor fiber) population. Dissertation. Swedish University of Agricultural Sciences, Department of Wildlife Ecology, Umea, Sweden. Hochachka, W. M., D. Fink, R. A. Hutchinson, D. Shel- don, W. K. Wong, and S. Kelling. 2012. Data- intensive science applied to broad-scale citizen science. Trends in Ecology and Evolution 27:130–137. Hooten, M. B., C. K. Wikle, R. M. Dorazio, and J. A. Royle. 2007. Hierarchical spatiotemporal matrix models for characterizing invasions. Biometrics 63:558–567. Hostetler, J. A., and R. B. Chandler. 2015. Improved state-space models for inference about spatial and ❖ www.esajournals.org 14 September 2017 ❖ Volume 8(9) ❖ Article e01947 BROMMER ET AL. temporal variation in abundance from count data. Ecology 96:1713–1723. Isaac, N. J. B., and M. J. O. Pocock. 2015. Bias and information in biological records. Biological Jour- nal of the Linnean Society 115:522–531. Kauhala, K., and T. Turkia. 2013. Habitat use of bea- vers: preliminary comparison between a native and alien species. Suomen Riista 59:20–33. Kelling, S., D. Fink, F. A. La Sorte, A. Johnston, N. E. Bruns, and W. M. Hochachka. 2015a. Taking a ‘Big Data’ approach to data quality in a citizen science project. Ambio 44:601–611. Kelling, S., et al. 2015b. Can observation skills of citi- zen scientists be estimated using species accumula- tion curves? PLoS ONE 10:e0139600. Kery, M., andM. Schaub. 2012. Bayesian population anal- ysis using WinBUGS: a hierarchical perspective. Aca- demic Press, Elsevier, Amsterdam, The Netherlands. Knape, J., and P. de Valpine. 2012. Are patterns of density dependence in the Global Population Dynamics Database driven by uncertainty about population abundance? Ecology Letters 15:17–23. Kujala, H., V. Veps€al€ainen, B. Zuckerberg, and J. E. Brommer. 2013. Range margin shifts of birds revis- ited–the role of spatiotemporally varying survey effort. Global Change Biology 19:420–430. Lahti, S., and M. Helminen. 1974. The beaver Castor fiber (L.) and Castor canadensis (Kuhl) in Finland. Acta Theriologica 13:177–189. Law, A., K. C. Jones, and N. J. Willby. 2014. Medium vs. short-term effects of herbivory by Eurasian beaver on aquatic vegetation. Aquatic Botany 116:27–34. M€uller-Schwarze, D. 2011. The beaver: its life and impact. Cornell University, Ithaca, New York, USA. Nolet, B. A., and F. Rosell. 1994. Territoriality and time budgets in beavers during sequential settlement. Canadian Journal of Zoology 72:1227–1237. Pagel, J., and F. M. Schurr. 2011. Forecasting species ranges by statistical estimation of ecological niches and spatial population dynamics. Global Ecology and Biogeography 21:291–304. Parmesan, C., and G. Yohe. 2003. A globally coherent fingerprint of climate change impacts across natu- ral systems. Nature 421:37–42. Pebesma, E. J. 2004. Multivariable geostatistics in S: the gstat package. Computers and Geosciences 30: 683–691. Plummer, M. 2003. JAGS: a program for analysis of Bayesian graphical models using Gibbs sampling. Pages 125–134 in K. Hornik, F. Leisch, and A. Zei- leis, editors. Proceedings of the Third International Workshop on Distributed Statistical Computing (DSC 2003), March 20–22, Vienna, Austria. R Core Team. 2015. R: a language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. Ricker, W. E. 1954. Stock and recruitment. Journal of the Fisheries Board of Canada 597:559–623. Root, T. L., J. T. Price, K. R. Hall, S. H. Schneider, C. Rosenzweig, and J. A. Pounds. 2003. Fingerprints of global warming on wild animals and plants. Nature 421:57–60. Royle, J. A. 2004. N-mixture models for estimating population size from spatially repeated counts. Biometrics 60:108–115. Royle, J. A., and R. M. Dorazio. 2008. Hierarchical modeling and inference in ecology. Academic Press, Elsevier, Amsterdam, The Netherlands. Sekercioglu, C . H., G. C. Daily, and P. R. Ehrlich. 2004. Ecosystem consequences of bird declines. Proceed- ings of the National Academy of Sciences of the USA 101:18042–18047. Silvertown, J. 2009. A new dawn for citizen science. Trends in Ecology and Evolution 24:467–471. Solymos, P., S. Leleand, and E. Bayne. 2012. Conditional likelihood approach for analyzing single visit abun- dance survey data in the presence of zero inflation and detection error. Environmetrics 23:197–205. Sutherland, W. J., editor. 1996. Ecological census tech- niques: a handbook. Cambridge University Press, Cambridge UK. SYKE. 2/2015. National datasets (25 m, 1 ha) SYKE (partly METLA, MMM, MML, VRK) EU datasets (25 ha, 5 ha) EEA, SYKE. http://www.syke.fi/avoin tieto Williams, B. K., J. D. Nichols, and M. J. Conroy. 2002. Analysis and management of animal populations. Academic Press, San Diego, California, USA. Willson, L. 1971. Observations and experiments on the ethology of the European beaver (Castor fiber L.). Viltrevy 8:115–306. Ylitalo, E., editor. 2013. Finnish Statistical yearbook of forestry 2013. Finnish Forest Research Institute, Vammalan Kirjapaino, Sastamala, Finland. SUPPORTING INFORMATION Additional Supporting Information may be found online at: http://onlinelibrary.wiley.com/doi/10.1002/ecs2. 1947/full ❖ www.esajournals.org 15 September 2017 ❖ Volume 8(9) ❖ Article e01947 BROMMER ET AL.