Large carnivores have been increasing generally in Europe and North America during recent decades (Chapron et al. 2014, Ripple et al. 2014). Their expansion into human-dominated landscapes has entailed increased damages to livestock husbandry and concerns for human safety (Kaczensky 1999, Linnell et al. 2001, Penteriani et al. 2016). Large carnivore populations are regulated through legal hunting in some European countries. In Finland and Sweden, hunting is an action inseparable from management of the most numerous carnivores, such as the brown bear Ursus arctos (Swenson et al. 1995, Kojola and Heikkinen 2006).

The brown bear has a low reproductive rate and is sensitive to high harvest rates (Zedrosser et al. 2001, Bishof et al. 2018). Owing to late maturation and long reproductive lifespan, generation time is far longer than for other large carnivores in Europe (Tallmon et al. 2004). Another life history character, meaningful for the sustainability of management, is natal philopatry of females (Swenson et al. 1998, Støen et al. 2006).

Due to human-caused mortality in 1800s and early 1900s, the era when bears were principally considered as a threat to livestock, the bear population in Fennoscandia gradually decreased and split into Scandinavian and Karelian subpopulations (Curry-Lindahl 1972, Pulliainen 1990, Swenson et al. 1995). During the last 50 years, both Scandinavian and the western segment of the Karelian population, and Finland's bears have increased in size and distribution (Chapron et al. 2014). The northern part of Fennoscandia is a reindeer Rangifer tarandus husbandry region, where bears cause substantial damage by killing calf reindeer (Sivertsen et al. 2016). Therefore, the human-caused brown bear mortality rate in northernmost Finland has been higher than in south of reindeer husbandry region (Kojola and Heikkinen 2006).

Bear densities might be low in Finnish reindeer husbandry areas owing to higher harvest rates (Kojola and Heikkinen 2006), which might reduce gene flow between the Karelian and Scandinavian populations. Former analyses in Finland indicated that the population structure shifts to increasing male biased from east to west (Kojola and Laitala 2000) and from south to north (Kojola and Heikkinen 2006). Population densities are highest in eastern Finland (Kojola et al. 2006) and peripheral Finnish populations have been receiving immigrants from Russia during the last decades, which is indicated through increasing genetic diversity among Finland's bears (Hagen et al. 2015). However, immigrants might be mostly sub-adult males that are most prone to move far from their natal home ranges (Støen et al. 2006).

Immigration rates between the two bear populations in Fennoscandia (Karelian and Scandinavian) are very likely to be linked to the past human persecution. Historically the distribution range of brown bear covered almost the entire forest zone of Europe but human persecution resulted in population declines and disconnections (Curry-Lindahl 1972). In Fennoscandia brown bear populations in Finland and Scandinavia has shown substantial recovery during the last 50 years (Chapron et al. 2014, Swenson et al. 2017) but genetic analyses provide evidence that gene flow between Karelian and Scandinavian populations are weak (Kopatz et al. 2014, Schregel et al. 2015). While there has been increasing influx of brown bears from Scandinavia into Finland, migration from Finland into Scandinavia appeared to be very low (Kopatz et al. 2019).

The management plan for Finland's brown bear population (Ministry of Agriculture and Forestry 2007) proposed that regional harvest plans for reindeer husbandry region (the only district that is bordering Scandinavia) should enhance immigration from Finland to Scandinavian population. In this study we examined changes in local population structure in Finland's brown bear population to see whether the proportion of adult females that produce potential immigrants from Finland to Scandinavia was increasing. The second question was the expected higher proportion of males by geographic distance from suggested core areas in northwestern Russia. Furthermore, sex ratio along the distance gradient might be changed because of doubled population size (cf. Liukko et al. 2016) might accelerate natal dispersal by females (Støen et al. 2003). Because female dispersal in the brown bear may increase before the population reaches saturation density (Støen et al. 2006, Zedrosser et al. 2007), we might expect the male bias to become less pronounced near the edges of the population over time. On the other hand, bears in Finland live mostly in human-dominated landscapes, where intentional or accidental killing by humans (Kojola et al. 2006) may have a stabilizing effect on the potential changes in population structure (Jerina and Adamic 2008, Krofel et al. 2012). The time window in our study is 21 years, which comprises two brown bear generations (Tallmon et al. 2004).

Data

We analyzed data on age and sex of dead brown bears in Finland during 1996–2016 (Fig. 1). Most of these 1225 bears were killed by legal hunting (93.3%), management removals of problem individuals (2.6%), traffic collisions (1.5%) and for other reasons (2.6%). The brown bear hunting season was from 20 August to 31 October and harvest was limited by regional quotas. Hunters are encouraged to fill out a form about every bear they kill. The form includes information about the geographic coordinates of the kill site and the sex of the bear, and is sent to Natural Resources Institute Finland (Luke) along with the second premolar tooth and a piece of muscle from the harvested bear. Age determination was based on cementum annuli in the root of the premolar and was executed in Matson's Laboratory, Montana, USA. Finland is divided into four zones for the management of the brown bear population (bear management area BMA, Fig. 1), based on estimated differences in bear population density and structure (Ministry of Agriculture and Forestry 2007). The BMAs will be referred to as; stable population area (STAB), area of dispersing population (DISP), area of developing population (DEVE) and reindeer husbandry area (REIN). We combined the bears' age into subadults (1–4 years for males, 1–3 years for females) and adults (>4 years for males, >3 years for females, Swenson et al. 1995).

Statistical modeling

Dependent variable was Bernoulli distributed variable with the event ‘dead bear is female'. We constructed generalized linear mixed model (GLMM) where the tested independent variables (predictors) were 1) the age class of the bear, 2) bear management area BMA, 3) distance (km) from the Russian border and 4) year. We treated these variables and their two-way interactions as fixed effects. Geographic coordinates (EUREF standard, divided by 1000) were used as descriptor of the spatial correlation structure. Two-way interactions were tested by two different ways. We started from the model with all possible combinations and then, step by step dropped one, two, three etc. combinations from the model. Then we constructed models where each two-way interaction term was entered one by one into the model. We tried also models with three-way interactions but these models did not converge properly.

We tested two correlation structures in the logistic mixed effects model for the lowest level observations (dead bears nested within year). The correlation structures could not be tested simultaneously in the same model, because R package MASS and its function glmmPQL (Venables and Ripley 2002) did not allow the definition of two different correlation structures in the various levels in the same model. In the first stage a longitudinal autocorrelation (AR1) within years was tested. A feature in R package MASS (Venables and Ripley 2002) is that it does not compute the significances for the random effects or the autoregressive coefficients. However, it computes the 95%'s confidence intervals for the parameters if they were asked for. Furthermore, a generalized linear models using the function glmmPQL in package MASS uses a pseudo-likelihood estimation and does not compute AIC:s. However, the other mixed model packages in R have deficiencies in the definition of the different correlation structures (e.g. lme4, Bates et al. 2015).

In the first mixed model the year was treated as a random variable assuming autoregressive correlation structure. In the model phi value (estimate 0.045, 95% confidence interval –0.02 to 0.10) did not differ from zero. Therefore we used compound symmetry as the variance–covariance structure for year in the second model. The second model took account spatial autocorrelation and because the range estimate (0.09) and the 95% confidence interval (0.03–0.28) for the correlation structure indicated that it differed from zero, we considered the second model as the final one. Exponential spatial autocorrelation was used in this model (Eq. 1, Table 1).

Figure 1.

Brown bears dead in Finland during 1996–2016 in a map showing the Finnish Wildlife Agency (RFWA) regions and bear management areas (BMA; REIN=reindeer husbandry, DEVE=developing population, DISP=dispersal zone, STAB=stable population).

Considering the spatial autocorrelation in the model did not change the fixed effects estimates and their significances, but it might be theoretically reasonable to use it. Spatial autocorrelation (or autocovariate) models address spatial autocorrelation by estimating how much the response variable at any one site reflects response values at surrounding sites (Dorman et al. 2007). Theoretically, it could be assumed that the probability to find a dead female in a certain place increased, if there were another dead female in the neighborhood. This is because females are spatially organized into kin clusters (Støen et al. 2005).

We computed the model considering a possible spatial autocorrelation using R function glmmPQL in R package MASS (Venables and Ripley 2002). The logistic mixed model with spatial autocorrelation term can be described as (Eq. 1):

where y is the probability of the event, i.e. ‘dead bear is female'. Binomial (n,p) denotes the binomial distribution with parameters n describing binomial sample size, in our case the number of all the bears and p describing the proportion, or probability of the occurrence of dead female bear. ln is a logit-link function and f(.) describes the linear function with arguments X_ij (i.e. fixed predictors) and β (i.e. fixed parameters). Term µ_i represents a random year effect and p is the coefficient of the autocovariate A_ij (called as range, estimated as a result of variance–covariance matrix, assumed to account for spatial exponential structure at the level ij denoting bears nested within years). The dispersion (cf. Browne et al. 2005) in the logistic model was estimated and expressed in the term µ_ij.

Table 1.

Parameter estimates and type III analysis of deviances for Finland's brown bears, 1996–2017. χ² test value was used for the analysis of deviance, and t-value for the test of the estimated coefficients. SE denotes standard error of the estimate, df the degrees of freedom and p the significance (probability) of the test. 95% CI denotes 95% confidence intervals. The area under ROC of marginal model (fixed effects only) is 0.60 and for marginal model (fixed and random effects) is 0.61.

We computed a receiver operating characteristics (ROC) curve to determine the classification efficiency of the logistic model. If the value under the ROC curve approaches 1, it means that the logistic model can predict the events (coded as 1 in the binomial response variable) adequately. The ROC curve was computed using the R package ROCR (Sing et al. 2005).

The model tested fixed effects by using ANOVA type III tests. Only significant (p < 0.05) independent variables and their interactions were presented in Table 1. We report analysis of deviance table with joined χ² Wald type III tests for deviances and t tests for the model intercept and for the coefficient of the age group difference from zero. All the analyses were done in the R (< www.r-project.org>).

Results

The GLMM model we conducted did not provide evidence of a significant trend in sex ratio (the probability of which the nearest shot bear is a female) over the study period (χ² = 1.73, df = 1, p = 0.188, Fig. 2). Two-way interaction between year and bear management area (BMA) on sex ratio was not significant (χ² = 5.38, df = 3, p = 0.146).

Both in the first mixed model (the year as a random variable assuming autoregressive correlation structure) and the final mixed model taking account spatial autocorrelation, shown in Table 1, the significant fixed effects were age group and BMA, but distance from Russian border was not related to the probability that the nearest dead bear was a female (p = 0.163). Interaction between BMA and distance was not significant (χ² = 1.02, df = 3, p = 0.786). We tested all other potential two-way interaction terms between fixed effects (age group, BMA, distance, year) but any of them were not significant (p values >0.05).

Sub-adult males predominated in our data (Fig. 3) and the probability that the nearest dead bear was a female was highly significantly related to the age group (Table 1). The probability was lowest within the reindeer husbandry area (REIN, Fig. 4). The reindeer husbandry region differed highly significantly from two other bear management areas (stable population STAB and dispersal zone DISP) while not from the developing population (DEVE, Table 1, Fig. 4).

Figure 2.

Annual sex ratio in brown bears dead in Finland during 1996–2016.

Figure 3.

Age distribution of male and female brown bears dead in Finland during 1996–2016.

Figure 4.

Probabilities that the nearest dead brown bear was a female within different bear management areas (a) and by the bear age class (b) during 1996–2016 in Finland (REIN=reindeer husbandry, DEVE=developing population, DISP=dispersal zone, STAB=stable population).

Discussion

The proportion of sub-adult males usually increases and the proportion of females decreases from the core toward the periphery in brown bear populations (Jerina and Adamic 2008, Krofel et al. 2012). The probability that the nearest dead bear was a female was linked to age class which is probably related to the higher risk of sub-adult males than females to get harvested by hunters (Kojola et al. 2003, Krofel et al. 2012). One reason for the higher risk are broader movements by sub-adult males than sub-adult females Støen et al. 2003, 2006, Zedrosser et al. 2007). The low proportion of females among dead bears provides evidence that the reindeer husbandry region is a peripheral part of Finland's brown bear population. Owing to females' higher natal philopatry, peripheral areas in bear populations are male-biased (Swenson et al. 1995, 1998, Støen et al. 2006). The low proportion of females in northern Finland is likely to have a negative effect on metapopulation viability of brown bears in Fennoscandia because gene flows are dependent on migration. Migration from the natal range is affected by population density (Swenson et al. 1998, Støen et al. 2005, 2006, Zedrosser et al. 2007), and connections between two populations by the distance between natal home ranges of potential migrants (Linnell et al. 2005, Kojola et al. 2009).

Occasionally males can move far from their natal home range (Barton et al. 2019), but genetic differentiation indicates low overall gene flow between Finland and Scandinavia (Kopatz et al. 2014, Schregel et al. 2015). Recent genetic analyses (Kopatz et al. 2019) found low levels of gene flow from Finland into the Scandinavian population. Contrary to our expectation, the sex of the bear was not related to the distance from the Finnish–Russian border. An earlier analysis, using east–west gradient as a measure of periphery, showed an increasing proportion of males from east towards west (Kojola and Heikkinen 2006). Our analysis using the actual distance from the Russian border did not, however, provide evidence that the degree of periphery was related to the distance.

Another unexpected result was the lack of interaction between the distance from the Russian border and year on the sex of the harvested bear, which provided evidence that population structure stayed spatially unchanged during our 21-year study period. This stability might be due to harvest, because the proportion of females in the hunter harvest may have increased with increasing harvest rates (Kojola et al. unpubl.).

The peripheral nature of Finland's brown bear population and potential gradual increase of the proportion of females were both probably influenced by translocations from easternmost Finland to mid-Finland, which were carried out in the 1980s. Genetic analyses linked all females shot in mid-Finland during 1990s and early 2000s to one particular female (Saarma and Kojola 2007), which was released in 1982 and shot in 1995 (Kojola, an unpublished record).

The borders between management areas in Finland are mostly based on province borders. The structure of management, where provincial offices of the Finnish Wildlife Agency have a key role in implementation of brown bear management policy in their provinces, is a major reason for this situation. There are two exceptions from this rule, both from eastern Finland (Fig. 1). In two of the management areas, the westernmost DEVE and the reindeer area (REIN), females seemed to concentrate in the eastern segment of the management area (Fig. 1). Such a heterogeneity within the management zone doubts biological relevance of the current division into management areas.

The low proportion of females in the reindeer husbandry area might be due to a higher relative immigration rate by subadult males, because population density was higher in two neighboring areas, northwestern Russia in the east (Kojola et al. 2003) and northern Sweden in the west (Kindberg et al. 2011). Female philopatry is likely to be one major reason why the most heavily harvested northern bear population in Finland is more male-biased than in other management zones. Population density in the reindeer area (REIN) is much lower compared to the locations with the similar distance from Russian border in the south (Kojola et al. unpubl.). The combination of female philopatry, low population density and the low proportion of females is very likely leading to the situation where migration to Scandinavia is very low (cf. Kopatz et al. 2019).

Our analyses documented stability in local population structure of a hunted brown bear population. Finland's bear population is commonly regarded as a peripheral northwestern outskirt of the Eurasian meta-population (Pulliainen 1990, 1997, Kojola and Laitala 2000, Kojola et al. 2003, Kojola and Heikkinen 2006, Tammeleht et al. 2010, Kopatz et al. 2012, Schregel et al. 2015). We found no significant temporal autocorrelation in sex ratios over the years; only a slight, statistically insignificant relationship existed (Fig. 2). This seems to be at least partially due to management interventions in the form of female translocation, the changes we expected to occur in local population structure in the peripheral brown bear population in Finland during 21 years did not occur. Another potential reason, increased natal dispersal by females owing to increased population density remained unlikely because the proportion of females did not change although population size increased.

Based on the lacking sex ratio gradient by the distance from Russian border, the peripheral nature of Finland's brown bear population does not seem as obvious as has been formerly thought (Pulliainen 1997, Kojola and Laitala 2000, Kojola et al. 2003, Kojola and Heikkinen 2006). Bear densities in northwestern Russia have been much higher than in Finland (Kojola et al. 2003) which probably enhances immigration into Finland. The increased genetic diversity south from reindeer husbandry area (Hagen et al. 2015) supports this expectation but may also decrease the difference in sex ratio between eastern and western Finland because usually most immigrants are sub-adult males.