Number, sex, and age data of brown bear kills used in the simulation

Hunting statistics spanning from 1969 to 1982 were compiled using records from the Hokkaido Nature Conservation Division; forestry divisions of the Oshima, Hiyama, and Shiribeshi subprefectural offices; and forestry sections of the municipalities in the study area. Starting in 1983, in addition to hunting statistical data, we obtained sex and age information by counting cementum layers of teeth (Mano 1987, 1995; HIES 2000), as depicted in Figures 2 and 3. Before 1982, only the total number of harvested bears was available because sex was not recorded.

Fig. 3.

Brown bear (Ursus arctos) harvest in the Oshima Peninsula, Hokkaido, Japan, for 1969–2021. Only total kills were available from 1969 to 1982.

Between the late 1960s and 1986, when the Spring Bear Removal program was implemented, the average number of bear deaths was approximately 100, with annual fluctuations. In 1986, the number of bears killed rose to 136, but decreased in 1987, with an average of 60 kills annually from 1987 to 2000. Subsequently, it increased again, reaching 183 in 2005. Further, the number of kills continued to increase, albeit with annual fluctuations, and in 9 of the 13 years from 2009 to 2021, annual kills exceeded 150. In 2020, 211 bears were killed, marking the highest number recorded (Fig. 3). Detailed harvest counts are shown in Appendix 1 (  Supplemental material (URSUS-D-23-00016R1._Mano_et al._SUPPLEMENTAL_MATERIAL_AND_APPENDICES.docx)).

From 1983 to 1989, when Spring Bear Removal was allowed, 56% of the bears hunted were male and 44% were female. From 1990 to 2021, after the abolition of Spring Bear Removals, males accounted for 67% of the kill, double that of females (Fig. 3).

We generated sex and age data for bears killed each year for computer simulations through the following process. Since 1983, surveys on the sex and age of killed bears were continuously conducted. Although the 0.4% of the killed bears with unknown sex, we randomly assigned sex with a 50% probability for these cases. The timing and methods of killing changed before and after the abolition of the Spring Bear Removals in 1990, and the male sex ratio in killing increased (Mano 1998); therefore, we divided the study period into before and after 1989–1990 and generated sex–age data for the killed bears.

Although the percentage of bears that were age-determined varied from 36% to 98% annually, the overall determination rate for 1983–2020 was 69%. We determined the ages of bears whose teeth were not recovered and whose ages were unknown randomly based on the distribution of age-determined bears in each of the 2 periods (Appendix 2,   Supplemental material (URSUS-D-23-00016R1._Mano_et al._SUPPLEMENTAL_MATERIAL_AND_APPENDICES.docx)). The sex ratio and age structure of the harvest before 1989, when the Spring Bear Removal was conducted, did not appear to have changed significantly (Mano 1987), so we determined the age structure for 1969–1982 harvest randomly based on the pooled sex and age structure for 1983–1989 (Appendix 3,   Supplemental material (URSUS-D-23-00016R1._Mano_et al._SUPPLEMENTAL_MATERIAL_AND_APPENDICES.docx)).

Population trend

Spring Bear Removal from 1966 to 1989 caused the Hokkaido population size to decrease, with larger declines in the northern and western areas (Aoi 1987, 1990; Mano 1991; Mano and Moll 1999). In the 1980s, nearly all hunters of bears in spring in the Oshima Peninsula observed a decline in bear signs during the spring hunting season and agreed that the population declined as a consequence of the Spring Bear Removal (Mano 1998). Interviews with a skilled hunter who had hunted brown bears for >30 years in the northern part of the study area near the town of Oshamambe, indicated that the brown bear population had declined as of 1973 compared with earlier years (Niizuma 1973). Collectively, reports suggest that the brown bear population in the Oshima Peninsula declined during the Spring Bear Removal. However, this does not necessarily imply a decrease as pronounced as in the northern areas. Thus, we assumed that the population size in 1990, the year the Spring Bear Removal was abolished, did not exceed that in 1968, the first year of the simulation (Fig. 2).

After 1990, the number of bears killed as a result of conflicts increased and the bear distribution range expanded (HIES 2000, Japan Bear Network 2014). On the Oshima Peninsula, national forests cover 51% of the forest area, or approximately 2,880 km² (Fig. 1), and almost all of the area is brown bear habitat. In the national forest area, surveys were conducted between June and November in 1996–1998, 2002–2004, and 2012–2020 by forestry workers who recorded whether brown bears, related sign, or feces were noted. The number of days with brown bear information and the number of feces per day reported by forestry workers were higher after 2012 than before 2004 (Table 1, Fig. 4; Welch's t-test; P < 0.01 for days with confirmed information and P < 0.01 for the number of confirmed feces). Furthermore, awareness surveys of hunters conducted in the 1992 and 2012 suggested an increase in brown bear abundance in 2012 (Hokkaido Research Organization, Institute of Environmental Sciences 2013). Thus, we assumed that the population increased from 1990 to 2012 (Fig. 2).

Table 1.

Trends in brown bear (Ursus arctos) density indices based on bear signs and bear observations by forest worker crews of the national forests in the Oshima Peninsula region Hokkaido, Japan, during 1996–1998, 2002–2004, and 2012–2020.

Fig. 4.

Trends of the brown bear (Ursus arctos) density indices based on bear signs found and bear observations by forest worker crews of the national forests in the Oshima Peninsula (Hokkaido, Japan) region from 1996–1998, 2002–2004, and 2012–2020. (A) is bear information confirmation rate per 100 work crews and (B) is number of feces confirmed per 100 work crews, respectively.

Population dynamics model

The birth period of brown bears in Hokkaido ranges from late January to early February (Maeda and Ohdachi 1994). In the model, we assumed that all births occurred on 1 January and generated the population size just before birth (i.e., at the end of the year). The values for the number of females at age x, adult females, males at age x, and adult males on 31 December in year t are represented as N_f(t, x), N_f(t, k), N_m(t, x), and N_m(t, k), respectively (x = 0, 1. . . k – 1, where k is the age at maturity). The mean birth rate, harmonic mean of the birth interval in years, average litter size, and average sex ratio at birth were denoted as B, I*, M*, and S*, respectively. The average sex ratios at birth were defined as N_m(t, 0)/[N_m(t, 0) + N_f (t, 0)].

The following equations further define inputs to the model:

Let B_m be the mean male birth rate (mean no. of male births per year). Then,

Let B(t) be the birth rate in year t.

where f is the magnitude of annual demographic stochasticity; and Z(t) is a uniform random distribution between –1 and 1.

Next, let L_c, L_fs, L_fa, L_ms, and L_ma be the survival rates of the cub-of-the-year (COY), subadult (1–5 yr old) female, adult (≥6 yr old) female, subadult male, and adult male, respectively. The survival rate was defined as the natural mortality rate caused by factors other than human-caused mortality, such as control killing, hunting, and/or traffic accidents. C_m(t,x), C_ma(t), C_f(t,x), and C_fa(t) were the number of cubs and subadult males of x years of age (<6 yr), adult males (≥6 yr), cubs and subadult females of x years of age (<6 yr), and adult females (≥6 yr) killed in year t, respectively. For simplicity, we assumed that there was no correlation between susceptibility to control kill, hunting, or traffic accidents in a given year and the probability of death due to natural causes, such as disease or starvation, in the following year.

Assuming that animals k years old or older are adults and those between 1 and k years old are subadults, the number of females in each age class in a given year t is

In addition, the number of males by age class is

The total number of individuals is expressed as

Demographic parameters

There has been no observational or collateral evidence indicating a change in the rate of natural increase. Furthermore, there have been no significant differences in brown bear breeding parameters in Hokkaido since the 1980s (Mano and Tsubota 2002, Shimozuru et al. 2017, Moriwaki et al. 2017). However, to consider various possibilities, we conducted simulations assuming that demographic parameters change annually. Therefore, our analysis encompasses scenarios with both interdecadal changes and density dependence in the rate of natural increase.

We defined birth rate following Mano and Tsubota (2002) and Shimozuru et al. (2017). The brown bear population on the Oshima Peninsula has a mean litter size (M*) of 1.80 ± 0.26 cubs (mean ± standard deviation [SD]), with a 1:1 sex ratio (S*), and the birth interval was 2.3–3 years (Mano and Tsubota 2002). Although the age at first birth was estimated to be 4.0 years, cub survival rates for first litters were low. In contrast, litter size from continuous direct observation of females with young in the Shiretoko Peninsula was estimated at 1.76 ± 0.08 cubs, with a birth interval of 2.53 years and mean age at first birth of 5.3 years (Shimozuru et al. 2017). For the model, we assumed age at successful parturition (a) as 6 years and used 2.6 years as a mean reproductive interval. There were no significant differences between the Oshima Peninsula and Shiretoko Peninsula regions in litter size, birth interval, and age at first birth, so we employed a litter size of 1.8 cubs, birth interval of 2.6 years, and age at first birth of 6 years (k = 6) for the Oshima Peninsula region. We set the sex ratio to 0.5.

The major food items in autumn are nuts and berries (Sato et al. 2005), which show significant annual fluctuations (Masaki and Abe 2008). These fluctuations could cause annual fluctuations in birth rate, as reported in American black bears (Ursus americanus; e.g., Rogers 1976, Eiler et al. 1989, Elowe and Dodge 1989). However, the fluctuation magnitude is uncertain so we assumed a degree of birth rate fluctuation among years of 50% for convenience.

Sex- and age-based survival rates are difficult to estimate for bears, especially cubs. Bunnell and Tait (1985), Craighead et al. (1995), and Harris et al. (2007) estimated the cub survival rate for the first year after birth to be 65%. For the hunted brown bear population in Scandinavia, cub's survival was estimated as 59% (Gosselin et al. 2015). In contrast, Shimozuru et al. (2017) estimated the survival rate from 0.5 to 1.5 years of age to be 60–73% under conditions of no hunting or control killing in the Shiretoko National Park, Hokkaido. Using these findings as a reference, we assumed a 65% survival rate in the first year after birth.

The natural mortality of subadults is considered to be higher than that of adults and that of females to be higher than that of males (Bunnell and Tait 1985). Survival estimates from radiotracking adult and subadult brown bears approximated natural mortality rate at a few percent, assuming that human-caused mortality accounts for approximately half of the total mortality (McLellan et al. 2000: table 1).

For the Oshima Peninsula brown bear population, natural mortality, including of subadults and older, was estimated to be 6% for females and 9% for males (HIES 2000) but with large confidence intervals (females: 0–18%, males: 0–25%). Prioritizing simplicity in the calculation assumptions, we assumed that there was no significant difference in the natural survival rate in the age classes older than subadults or between males and females. Thus, we assumed the mean natural mortality rate to be the same and constant for both sexes from subadults to adults >1 year of age. We examined the natural mortality rate in 1% increments from 3% to 8%. We considered uncertainties by incorporating the interannual variation in the survival rate L(t) for year t, which was the same for subadults, adults, males, and females.

Initial population size and age structure

We conducted a computer simulation of the brown bear population starting from 1968. However, information regarding the initial population size was not available, so we assumed the initial number to be a uniformly distributed random number between 1 and 1,809 individuals. This upper limit of 1,809 was observed as of 2012, as described below. For the initial age structure, we assumed that 65% of the population consisted of adults, and we set the initial sex ratio to N_m (1968)/N (1968) = 7/17 (i.e., males are 70% of the female population for convenience).

Conditions for valid simulation results

Uncertainty and/or annual fluctuations of model parameters can result in a range of population dynamic patterns, so we excluded unrealistic cases (e.g., extinction or irruptive growth) and retained only cases that met the following conditions:

Stochastic population dynamics models tend to produce population irruptions that are unrealistic (Lewontin and Cohen 1969, Boyce et al. 2001). We avoided unrealistic population growth by setting an upper bound obtained by extrapolating the density estimate in 2012 obtained for an area of approximately 300-km² within the Oshima Peninsula region. This forested area was considered to have the highest population density in the entire region from the results of a sign survey in the 1990s (HIES, unpublished data). In the study area, the density of females in 2012, estimated using the spatially explicit mark–recapture method identifying individuals by genetic analysis of samples obtained from hair traps, was 0.215 bears/ km² (95% confidence interval [CI] = 0.141–0.327 bears/ km²; Tsuruga et al. 2014). We extrapolated the 95% confidence upper limit for this density to a distribution area of 5,531 km² for brown bears in the Oshima Peninsula region, and used the resulting number of 1,809 bears as the upper limit for N_f (2012) in simulations (Fig. 2). We calculated 90% credible intervals (CIs), assuming that the estimated number of individuals followed a gamma distribution.

Sensitivity analysis of parameters' uncertainty

To examine the sensitivity of demographic parameters, we varied subadult and adult survival rates by 1%, ranging from 92% to 98%, and considered COY survival rates at 60%, 65%, and 70%. Regarding reproductive parameters, we investigated birth intervals of 2.3, 2.6, and 3.0 years. Table 2 provides an overview of the parameter combinations investigated.

Table 2.

Combinations of demographic parameters (with uncertainties) employed in the computer simulations used to build a population dynamic model constructed using long-term harvest data and estimated population dynamics of brown bears (Ursus arctos) on the Oshima Peninsula, Hokkaido, Japan, from 1968 to 2021. COY is cub-of-the-year.

Population and its trends

Various population dynamics patterns were obtained depending on the combination of parameters (Table 2). The number of trials needed to achieve 1,000 trials that met all the above 3 conditions is denoted as T in Tables 3 and 4. First, in general, under conditions of high birth rates, the higher the COY survival rates and subadult and adult survival rates, the lower the initial population in 1968 and the higher the population in the final year of the simulation in 2021. Lower survival and birth rates tended to cause milder population fluctuations over time, whereas higher survival and birth rates caused the population to increase rapidly after the 1990s (Figs. 5 and 6). In almost all combinations, the population decreased from 1968 to 1986, then increased, and then decreased in 2005. Thereafter, some high-survival and high-birth rate combinations increased again, whereas others decreased (Figs. 5 and 6).

Table 3.

Estimated population sizes for 1968, 1990, 2012, and 2021 for various brown bear (Ursus arctos) subadult and adult survival rate (L) scenarios (calculated from data collected on the Oshima Peninsula, Hokkaido, Japan, from 1969 to 2021) under conditions of 65% COY survival (L_c) and 2.6-year birth interval (I). N_f represents females, N_m represents males, and N represents the total population, respectively. Calculation results were truncated to the decimal point.

Table 4.

Estimated population sizes for 1968, 1990, 2012, and 2021 for various combinations of brown bear (Ursus arctos) cub-of-the-year survival rate (L_c) and birth interval (I) scenarios (calculated from data collected on the Oshima Peninsula, Hokkaido, Japan, from 1969 to 2021) under conditions of 94% subadult and adult survivals (L). N_f represents females, N_m represents males, and N represents the total population, respectively. Calculation results were truncated to the decimal point.

Fig. 5.

Population dynamics of brown bears (Ursus arctos) on the Oshima Peninsula (Hokkaido, Japan) from 1968 to 2021 under various subadult and adult survival rate scenarios with a cub-of-the-year (COY) survival rate of 0.65 and birth interval of 2.6 years obtained by simulations. The gray area indicates the range of the 90% credible interval (CI).

Fig. 6.

Population dynamics of brown bears (Ursus arctos) on the Oshima Peninsula (Hokkaido, Japan) from 1968 to 2021 under various cub-of-the-year (COY) survival rate (cs) and birth interval (bi) scenarios with subadult and adult survival rates of 0.94. The gray area indicates the range of the 90% credible interval (CI).

The computer simulation was conducted with the initial values of females as of 1968, ranging from 1 to 1,809, but the estimation results were within a certain range for each parameter combination. The credible intervals (CIs) for high survival and birth rates remained small during the decline from 1968 to the 1980s. However, they expanded rapidly after the 1990s. In contrast, the lower survival and birth rates the milder the CI expansion after the 1990s (Figs. 5 and 6).

Effects of subadult and adult survival rates

To examine the effect of subadult and adult survival rates on population estimation, we compared the populations in 1968, 1990, and 2012, when the model conditions were set with 2021 as the final year. We compared subadult and adult survival rates ranging from 92% to 97% in 1% increments with 65% COY survival and a 2.6-year birth interval (Table 3, Fig. 5). At the median trend with 94% subadult and adult survival, the total of 1,107 bears (634 females, 473 males) in 1968, decreased to 916 (594 females, 322 males) in 1990, respectively, then increased to 1,715 bears (1,103 females, 612 males) in 2012, and to 2,030 (1,368 females, 662 males) in 2021 (Table 3 and Fig. 5; L = 0.94). Furthermore, at 97% subadult and adult survival total population was 781 (447 females, 334 males) in 1968, 580 (412 females, 168 males) in 1990, and increased to 1,508 (1,045 females, 463 males) in 2012 and 2,475 (1,681 females, 794 males) in 2021 (Table 3 and Fig. 5; L = 0.97). At the median trends with 92% subadult and adult survival rates, the total population size decreased from 894 females and 668 males for a total of 1,562 bears in 1968 to 832 females, 600 males, and 1,432 in total in 1990. The population increased to 1,963 bears (1,187 females, 776 males) in 2012, and then decreased to 1,803 bears (1,206 females, 597 males) by 2021 (Table 3 and Fig. 5; L = 0.92). Thus, the abundance of each age group differed depending on the subadult and adult survival rates; however, the differences in the estimates were small around 2016 (Fig. 7). Specifically, the 95th percentile values in 2014 and the median in 2016 are nearly identical for all the survival rates. The 5th percentile values are also identical in 2019 (Fig. 7).

Fig. 7.

Comparisons of 95th percentile, median, and 5th percentile of the estimated brown bear (Ursus arctos) population dynamics on the Oshima Peninsula from 1968 to 2021 for various subadult and adult survival rate scenarios under conditions of cub-of-the-year (COY) survival rate of 0.65 and birth interval of 2.6 years. Auxiliary scales are indicated by dashed lines in 2014.

Effects of COY survival and birth rates

Similarly, we compared the impact of COY survival and birth rate on the estimation. Estimated results for 1968, 1990, 2012, and 2021 were compared for the subadult and adult survival rates of 94% (Table 4). The female, male, and total median values under conditions of 60% COY survival and a 3-year birth interval scenario were 1,468 (females 840, males 628) in 1968; 1,289 (f. 791, m. 498) in 1990; 1,894 (f. 1,197, m. 697) in 2012; and 1,871 (f. 1,294, m. 577) in 2021. Furthermore, under the high survival and birth rate scenario conditions of 70% COY survival and a 2.3-year birth interval, the median females, males, and total were 520, 387, and 907 in 1968; 481, 248, and 729 in 1990; 1,062, 583, and 1,645 in 2012; and 1,510, 803, and 2,313 in 2021 (Table 4).

Comparing the population estimation based on COY survival and birth rates, the 95th percentile had the smallest difference in estimates around 2013, the median around 2016, and the 5th percentile around 2018 (Fig. 8).

Fig. 8.

Comparisons of 95th percentile, median, and 5th percentile of the estimated brown bear (Ursus arctos) population dynamics on the Oshima Peninsula (Hokkaido, Japan) from 1968 to 2021 for various cub-of-the-year (COY) survival rate and birth interval scenarios under conditions of subadult and adult survival rates of 0.94. Auxiliary scales are indicated by dashed lines in 2014.

Characteristics of parameters meeting the conditions

The intrinsic rate of natural increase r(t) for each year is calculated backward from ln{N_f(t + 1)/[N_f(t) – C(t)]}. We proceeded to analyze the characteristic of successful simulation trials. The intrinsic rate of natural increase, denoted as r(t), tended to be higher values in 1968 and 1990 across 1,000 trials (Fig. S1,   Supplemental material (URSUS-D-23-00016R1._Mano_et al._SUPPLEMENTAL_MATERIAL_AND_APPENDICES.docx)); however, it is essential to consider the possibility that this observation could be an artifact.

In just 503 trials (50.3%), we observed a negative correlation between the number of individuals (N(t)) in successful trials and r(t). However, even when exclusively considering trials with this negative correlation, the results of population estimation remained largely consistent (refer to Fig. S2,   Supplemental material (URSUS-D-23-00016R1._Mano_et al._SUPPLEMENTAL_MATERIAL_AND_APPENDICES.docx)). The range of population fluctuations tended to be smaller.

Furthermore, we assumed that population trends after 2012 are unknown. However, we explored a scenario where we introduced a condition indicating a population increase in 2021. Among the 1,000 successful trials, the population increased in 883 trials. Under this condition, the population estimates focused on the median, 5th-percentile, and 95th-percentile values in 2021 were 2,033, 1,416, and 3,165, respectively. Notably, the lower limit was adjusted upward for the increasing trend (refer to Fig. S3A,   Supplemental material (URSUS-D-23-00016R1._Mano_et al._SUPPLEMENTAL_MATERIAL_AND_APPENDICES.docx)).

For the 117 scenarios that did not increase, the median, 5th percentile, and 95th percentile values were 1,374, 1,024, and 1,772, respectively. Under the assumption of no population increase, the lower limit only saw a slight decrease, while the upper limit was significantly revised downward (Fig. S3B,   Supplemental material (URSUS-D-23-00016R1._Mano_et al._SUPPLEMENTAL_MATERIAL_AND_APPENDICES.docx)). This suggests that population estimates obtained through this method may be subject to revision if additional information confirming a population trend post-2012 becomes available.

Characteristics of simulations

We conducted simulations using independent random numbers for 5 parameters: number of offspring, birth interval, survival rate of COY (Cubs of the Year), subadults aged 1–5 years old, and adults. The uncertainty for each parameter was assumed to be ±5%, and the overall magnitude of process uncertainty per year, considering the square root of the number of independent random variables, was approximately ±14%. This process was repeated for 50 years. Many simulations resulted in either too high or too low rates of natural increase. As a result, <0.6% of all trials met the criteria for the population trend (T^a >165,000 in Table 3).

We imposed qualitative conditions regarding population trends and selected only trials that met the 3 specified conditions. Similar to Approximate Bayesian Calculation (ABC), this method lacks knowledge of the likelihood of each trial, rendering it unable to efficiently perform a parameter search, as seen in the Monte Carlo Markov Chain Method (MCMC; Frost et al. 2023). In ABC, a prior distribution of all parameters and the initial population size is provided. However, in this study, average values of the parameters were assumed, and the robustness of the results was examined through sensitivity analysis.

Among the 1,000 successful trials, the lowest female population size in 1968 was 471. This implies that, despite assuming an initial female population size between 0 and 1,809, there was no need to explore initial population sizes of <471, which constituted more than one-quarter of all trials. Although there are numerous other parameter combinations where results may not be successful, searching within a broad range enables the elimination of preconceptions.

Different combinations of demographic parameters yielded different population estimates and patterns of increases or decreases over 40 years. However, despite the 3 conditions (i.e., the population in 1990 must not exceed the initial population in 1968, adults must not become extinct, and the population in 2012 must not exceed the upper limit obtained by the independent population estimation), we were able to narrow down the range of population size meeting these criteria in the final years of the computer simulation. The median, 5th percentile, and 95th percentile estimates all converged after the 2010s, at the end of the simulations compared with earlier periods (Figs. 7 and 8). These results demonstrate that, even when accounting for parameter uncertainties, certain inferences on population estimation can be drawn from the fact that harvesting lasted for >40 years and from the knowledge of qualitative population trends.

For the condition that the population declined as a result of the Spring Bear Removal, the upper limit was regulated by the initial population (1968) and the population trend was considered to have remained unchanged until the 1980s, without an increase in the margin of error. However, the population increased until 2012 when a new upper population limit was established and the credible interval (CI) continued to increase. After 2012, when no new upper limit existed, the CI widened (Figs. 5 and 6). Thus, the upward swing in estimates after 2012 would be more pronounced than before.

We set the initial sex ratio in this simulation for convenience. In this case, the initial female population size must be essentially important because the subsequent population dynamics would be mostly influenced by the initial female population size.

Brown bear population dynamics

Using only the results of computer simulations, it was difficult to select the most probable scenario among the valid population trends obtained. Although we were able to estimate the population size around 2014 within a relatively narrow range regardless of demographic parameters, we were unable to narrow down the specific values of the demographic parameters. As shown in Table 4, there is a negative correlation between the initial population size and survival rate and birth rate, but it remains unclear which combination is more realistic.

Bear-sign indicators in national forests increased in the 2010s compared with the pre-2000s; however, this alone made it difficult to estimate the rate of increase (Fig. 4). A series of results with scenarios of high subadult and adult survival rates suggest that the population increased rapidly in the 1990s to levels exceeding those of the pre-1970s. Certainly the population would have increased since 1990, but it has been growing unrealistically fast, and the population level around 2020 has been several times larger than that in the 1990s. In contrast, scenarios of low subadult and adult survival rates do not show extreme increases or decreases with a high median population size of approximately 2,000 (Fig. 5). It is uncertain how such gentle population fluctuations can be recognized by hunters at the Spring Bear Removal or forest workers. In other words, it is difficult to identify actual population dynamics based on these results alone.

Among the parameters examined, reproductive parameters, such as first birth age and birth rate, are easier to estimate than the survival rate, which is difficult to estimate with high precision (Eberhardt 1990). We considered default conditions of 65% COY survival and a 2.6-year birth interval, which are also close to the conditions found in the Shiretoko National Park (Shimozuru et al. 2017). Under these conditions, the high subadult and adult survival scenario resulted in a rapid increase from low abundance levels in the 1960s to the 1990s and beyond. In contrast, in the low subadult and adult survival scenario, the population remained high since the 1960s, and the increase or decrease was moderate (Fig. 5).

The most probable value for subadult and adult survival is important to determine. Natural mortality rates (N, 95% confidence interval) for subadult and adult female bears obtained by the radiotracking surveys conducted in the southwestern part of the study area were 6% (13, 0–18%) for 1987–1996 and 5% (11, 0–16%) for 1999–2003 (HIES 2000, 2004). The estimates remained virtually unchanged throughout the period. Adopting a subadult and adult mortality rate of 6% yielded the moderate estimates (Table 3; Fig. 5 [L = 0.94]). The 5th percentile of estimation suggests that the number of animals decreased from 962 in 1968 to 740 in 1990, then increased to 1,149 in 2012, and decreased to 1,124 in 2021. The number of animals has remained constant or even slightly decreased since the 2010s (Table 3; Fig. 5). Dividing these populations by the 5,531-km² forest area yields 36.7 bears/100 km² in median and a minimum density of 20.3 bears/100 km² as of 2021.

Challenges in estimation methods

Density estimation of bears using spatially explicit mark–recapture models has been widely practiced in Japan in recent years; however, its effective use in population estimation has been challenging (Tsuruga 2008, Ueuma and Nakayachi 2006, Morimitsu 2008, Yamauchi and Saito 2008). This is because density estimation surveys can only be conducted over a limited area (from tens to hundreds of km²) and it has not been clarified how these results can be used to estimate the entire population, which has a distribution area of several thousand square kilometers.

In this study, we proposed a method for utilizing the results to calculate the upper limits of computer simulations. The estimated density in areas that are higher-density than the average density in the population distribution area can be used as the upper limit for the estimation. However, it is necessary to obtain a new upper limit at regular time intervals to control the expansion of estimation errors.

Although the model incorporates process uncertainty through random fluctuations in demographic parameters, it did not consciously account for long-term changes in time series parameters. In general, for K-selected species, such as bears, the demographic parameters are expected to be stable (Bunnell and Tait 1981); but for the analysis of long-term trends over several decades, the demographic parameters may change with environmental changes. For example, a decline in the reproductive performance of grizzly bears (U. arctos) associated with the closure of garbage dumps was observed in Yellowstone National Park (Stringham 1986). Although there are currently no observations to date suggesting that such significant fluctuations occur in the brown bear population of the Oshima Peninsula region, monitoring is necessary. Moreover, given demographic fluctuations, although we employed uniform distribution to the number of births and deaths, it may be better to apply a binomial distribution. However, the number of individuals was sufficiently large in all trials and the effect of population fluctuations was considered negligible.

Furthermore, although the model does not explicitly account for density effects, long-term observational data suggest that litter size may be influenced by population density (Schwartz et al. 2006) or cubs' survival could be reduced by infanticide (Gosselin et al. 2015). Density effects were not considered in the model, so an unlikely population irruption could have occurred.

Would the accuracy of population estimates improve if the trend of population could be determined? The estimated range increased under increasing conditions from the 1990s to 2012 (Figs. 5 and 6); therefore, population trend information would be useful to obtain reliable population estimation (Fig. S3). Also, it would be necessary to periodically obtain an upper limit to remove excessive increase results by the stochastic model. In this calculation, there occurred a 9-year period during which the population could increase without an upper limit after the upper limit was set in 2012. It would be considered necessary to obtain a new upper limit with a 10-year interval, considering the duration of the bear management plan.

It should be noted that hunting statistics generally do not include records of poaching, and some kills are unreported (HIES 2000). However, this implies that the actual bear harvest number is likely higher than that used in our model, introducing a downward bias to our population estimates compared with the actual situation. Assuming a conservative approach to population management, such a bias is not necessarily a problem. The method we used did not require estimation of the harvest and natural mortality rates, which are difficult to obtain with high accuracy. Such mortality rate information would be useful for analysis of population dynamics where large numbers of bears have been harvested over long periods; however, intensive studies, such as large-scale mark–recapture or radiotracking methods, have not been conducted. In addition, it may also be applied to bears other than brown bears, for which information is limited to hunting statistics.

Management implications

Takinami et al. (2021) demonstrated a significant decline in the brown bear population in northern Hokkaido by the 1980s, followed by a gradual increase after the 1990s. In contrast, the decline estimated in this paper for the Oshima Peninsula area was relatively mild, indicating a rapid increase after the discontinuation of the Spring Bear Removal. This discrepancy is likely attributable to regional differences in the impact of the Spring Bear Removal.

Although the Oshima Peninsula has limited space, there are likely refugia in various locations where hunters' access is physically restricted by the steep mountainous terrain, making it impractical to eliminate all bears. On the other hand, the northern Hokkaido region reported by Takinami et al. (2021) features a relatively gentle hilly terrain with ample snow cover, facilitating efficient bear hunting. Additionally, the habitat in this region is probably less productive compared with that of the Oshima Peninsula, resulting in a lower population density. These observations align with the regional population estimates outlined in the Hokkaido Brown Bear Management Plan (Hokkaido Government 2022).

HIES (2000) estimated the population size, excluding COY, on the Oshima Peninsula by calculating harvest rates using radiotracking surveys and the number of kills since 1992. They estimated 380 females and 142 males, for a total of 522 individuals. In addition, according to a questionnaire survey of hunters in Hokkaido, the brown bear population size was estimated at 281–544 in 1992 and 379–685 in 2012. Inukai et al. (1985) estimated that the brown bear population in the region in 1980 ranged from 490 to 638. In this estimate, if we adopt the 5th percentile obtained under the scenario of subadult and adult survival rates of 94%, COY survival rate of 65%, and birth interval of 2.6 years, then the 5th percentile estimate in 1990, when the population declined from 1968, was 740 as compared with 1,149 in 2012 and 1,124 in 2021. Population size estimations from previous studies may have been underestimated.

We assumed that the population increased in 2012 compared with that in 1990. However, although the number of deaths has increased since then, it is impossible to determine whether the population has continued to increase or has declined since 2012, because there are no reliable results from longitudinal surveys on the population size or its index. Except for L = 0.92, in Figure 5, the upper limit values, including Figure 6, continue to increase; and except for L = 0.96 and 0.97 in Figure 5, the lower limits, including Figure 6, start to decrease. Interval estimates of populations resulted in lower bounds of about 1,050–1,250 individuals and upper bounds of about 2,900–4,500 individuals in all scenarios.

Nonetheless owing to differences in the natural increase rates, the lower bounds on the number of kills needed to reduce the population and the upper bound on the number of kills needed to avoid reducing the population differ considerably. Bayesian estimates of population size based on population dynamics models that consider catches usually show a negative correlation between population estimates and natural increase rates (Matsuda et al. 2002: fig. 7c). However, in the estimation method used in our study, the natural increase rate was high in the scenario that gives an upper limit of the population estimation, whereas the lower limit was low. This implies that there is a high degree of uncertainty in the target setting for the number of bears to be culled. Although increasing social pressure to enforce population control by killing to reduce the growing human–bear conflict has been noted in recent years (Uchiyama 2023), careful management based on harvesting and monitoring of the population must be necessary.