Large ungulates have a significant impact on vegetation in many parts of the world (Rooney and Waller 2003; Côté et al. 2004) and cause substantial agricultural and forestry damage (Conover 2001), thus requiring control by land managers and the government. Although there has been substantial discussion regarding the impact of large ungulates on vegetation and the appropriate density of ungulates (Mysterud 2006; Suzuki et al. 2013; Itô et al. 2014), population control by capture is necessary to avoid irreversible vegetation changes (Tanentzap et al. 2012). Therefore, if the government plans to reduce ungulate populations, it is necessary to estimate their abundance for management. The current population size, culling pressure to be applied to a target species, the duration of culling, and the number of animals to be culled to reduce the population to the target level are some of the factors to be considered when making predictions.

In recent years, the distribution area of sika deer (Cervus nippon) in Japan has expanded broadly (Yamane 2019), by approximately 2.7 times from 1978 to 2008 (Ministry of the Environment, Japan 2021a, b). Impacts of deer such as debarking on standing trees, decline of understory vegetation (e.g., Fujiki et al. 2010; Iijima and Nagaike 2015), and changes in plant species composition in semi-natural grasslands (Otsu et al. 2019) caused by abundant sika deer have been reported. As a countermeasure to mitigate the impacts of deer on vegetation, population control of sika deer is being conducted in Japan. For successful population management, the abundance of sika deer and its spatial and temporal changes during hunting and culling (hereafter, “capturing”) should be evaluated.

The harvest-based model (HBM) is used to estimate animal abundance under hunting pressure based on the response of the abundance index to the number of hunted animals (Iijima 2020). HBM was first used to estimate the abundance of sika deer in Hokkaido (Yamamura et al. 2008; Kaji et al. 2010). It is currently used in many prefectures in Japan (Iijima 2018). Because the HBM utilizes several abundance indices to estimate latent abundance and incorporates external data such as landscape components of demography to estimate the population growth rate, it is an integrated population model that estimates abundance and demographic parameters (Schaub et al. 2007; Kéry and Schaub 2011; Schaub and Abadi 2011). The merits of HBM in the abundance estimation of sika deer are that 1) the abundance of sika deer can be estimated by considering the effect of population control and monitoring data such as the abundance index, and 2) the model structure offers flexibility depending on the available datasets and their overdispersion terms due to the effects of unobserved factors. However, such flexibility requires a substantial understanding of the statistical modeling by the analyst to construct an appropriate HBM. Without such an understanding, the analyst may make an incorrect conclusion based on a biased estimation. For example, the use of the number of hunted animals as an abundance index in the HBM causes a seriously biased estimation of abundance (Fukasawa et al. 2020); however, such a model has been applied in the abundance estimation of sika deer in Japan (Iijima 2018). Yamamura (2016) also pointed out the possibility of biased estimates in the HBM using an inappropriate prior distribution depending on the available data.

The more serious problem is that government officials and planners accept only the results (especially the median or mean) without considering the structure of the model that produces the estimate, quantity, and quality of data used. Unfortunately, the significant impact of differences in the model structure, quantity, and quality of data used for abundance estimation are rarely understood by most government officials and planners in Japan. Sometimes, they simply recognize that “the abundance was estimated using an HBM.” Planning by trusting only the estimates obtained without understanding the structure of the model or mechanism by which the estimates are obtained involves the risk of misunderstanding the mechanism of population fluctuations in nature and choosing the wrong measures for population control. To identify the problems of using the HBM, it is necessary to understand how the structure of the model affects the parameters of the model and the results of the population estimates. Additionally, it is common for the monitoring data available at a given time to be insufficient. However, to promote wildlife management, wildlife managers must make the best use of all available data to obtain the most reasonable population estimates possible.

Another problem is that in Japan, there have been few cases of HBM with spatially fine-scale structure. Iijima et al. (2013) achieved a spatially fine-scale HBM based on hunting mesh (ca. 5 km × 5 km) to estimate the sika deer population in Yamanashi Prefecture. However, only Gifu Prefecture (2021) has used this model to plan prefectural management of sika deer. Other prefectures, such as Hyogo (Takagi 2019; Hyogo Prefecture 2022), Nagano (Nagano Prefecture 2021), and Chiba (Chiba Prefecture 2017), have also used different model structures to estimate sika deer abundance at a spatially fine scale; however, few prefectures have used such models. In addition, the Ministry of the Environment estimated sika deer abundance throughout Japan (excluding Hokkaido) without considering the spatial structure (Ministry of the Environment, Japan 2021a).

In this study, we clarified the impact of model structure and data quality on HBM estimates. To accomplish this objective, we estimated sika deer abundance using several spatially fine-scale HBMs with different structures and various deer abundance indices (hunters' report, number of captured deer, number of pellet groups, block count, and the relative abundance index based on camera trapping) collected over several years in Gifu Prefecture, Japan. We examined the effects of the model structure (primarily considering overdispersion in the observation models) on the estimated values with the aim of obtaining plausible population estimates that could contribute to wildlife management from currently available data. Based on the results, we discuss the structure of the model for a reasonable estimation. Furthermore, we evaluate the superiority of the results of a spatially fine-scale HBM compared to those of a traditional HBM.

Materials and methods

Process model

The process model of the spatially fine-scale HBM used in this study is as follows.

x_t,c is the sika deer abundance on a logarithmic scale in the c^th cell of the t^th year; hr_t,c is the hunting ratio in the c^th cell of the t^th year; xinit is the mean of sika deer abundance on a logarithmic scale in the first year (t = 1); ϕ_c is Gaussian conditional autoregressive (CAR) term; yhr_t is the yearly mean of the hunting ratio in the t^th year; r_c is the intrinsic population growth rate on a logarithmic scale in c^th cell; α is an intercept of r_c; β1 is a coefficient of the percentage of evergreen forests (Evergreen Forest Ratio; EFR) in the c^th cell; β2 is a coefficient of the percentage of deciduous forests (Deciduous Forest Ratio; DCR) in the c^th cell; β3 is a coefficient of the percentage of artificial grasslands (Grassland Ratio; GRR) in the c^th cell; and σs are the standard deviation of Gaussian distribution. The prior distribution of all σs in the process model is a vague uniform distribution, which is represented as U(0, 100) (Gelman 2006). The Gaussian CAR model is composed of a_cj, δ_c, and a_c+. a_cj is the indicator variable of the neighbors of cell c (a_cj = 1 if cells c and j share the same boundary and 0 otherwise). δ_c is the ID of cells that are neighbors of cell c and a_c+ is the total number of cells that are neighbors of c. We incorporated the spatial random effect into xinit because we predicted the spatial structure of deer abundance and had no prior information regarding deer abundance in previous years.

Because the population indices were obtained after autumn, it is possible that they had already been partially affected by the culling of the corresponding year. Therefore, it would be appropriate to delimit the year in the model from autumn to the next autumn if possible. Unfortunately, it was difficult to separate the old data in Gifu Prefecture into culling and hunting seasons. Therefore, the model in this study did not consider the change in population indices within a year and used various population indices obtained in autumn as representative values for the year.

Observation model

The observation model for the SPUE is shown below.

SPUE_t,c is the observation value of deer observed in the c^th cell of the t^th year; βSPUE is the coefficient of SPUE that links the true abundance (i.e., x_t,c) and the observed deer (i.e., SD_t,c); εSPUE_t,c is the overdispersion term in the c^th cell of the t^th year; Effort_t,c is the effort to obtain the SPUE in c^th cell of t^th year; and σ₆ is the standard deviation of the Gaussian distribution.

The observation model for the pellet group survey is shown below.

PG_t,c is the number of pellet groups in the c^th cell of the t^th year, βPG is the coefficient of the pellet group survey that links the true abundance (i.e., x_t,c) and the observed pellet group (i.e., PG_t,c), εPG_t,c is the overdispersion term in the c^th cell of the t^th year, Route_t,c is the route length of the pellet group survey in the c^th cell of the t^th year, and σ₇ is the standard deviation of the Gaussian distribution.

The observation model for the camera trapping survey is shown below.

CT_t,c is the number of photographed sika deer from September to October in the c^th cell of the t^th year; βCT is the coefficient of the camera trapping survey that links the true abundance (i.e., x_t,c) and the observed photographed sika deer (i.e., CT_t,c); εCT_t,c is the overdispersion term in the c^th cell of the t^th year; Day_t,c is the day length of the camera trapping survey from October to November in the c^th cell of the t^th year (usually 61 unless the camera malfunctions); and σ₈ is the standard deviation of the Gaussian distribution.

The observation model for the block count survey is shown below.

BC_t,c is the observation value of sika deer observed by a block count survey in the c^th cell of the t^th year, εBC_t,c is the overdispersion term in the c^th cell of t^th year, Area_t,c is the survey area of the block count survey in the c^th cell of the t^th year, and σ₉ is the standard deviation of the Gaussian distribution. As the results of the block count survey can be directly converted to deer density, the coefficient that links the true abundance was not included in the observation model in the block count survey, as done in Iijima et al. (2013).

The prior distribution of all σs in the observation model is a vague uniform distribution, represented as U(0, 100) (Gelman 2006).

The observation model for the number of captured deer is as shown below.

C_t,c is the number of captured deer in the c^th cell of the t^th year and D_t,c is the sika deer abundance in the c^th cell of the t^th year. It should be noted that the purpose of the observation model is to estimate hr_t,c, and C_t,c is not an abundance index.

Results

Data description

The number of captured deer differed significantly among cells (Fig. 1). Although the annual number of captured deer for the entire Gifu Prefecture was below 5000 in 2008, it increased to more than 17 000 in 2014 (Fig. 1). The increase in the number of captured deer was prolonged over the following years. The restriction of hunting in 2018 (partially in Gifu Prefecture) and 2019 (the whole of Gifu Prefecture) to prevent the spread of CSF decreased the number of captured deer compared to that of the previous year in each cell (Fig. 1).

Fig. 1.

Temporal change in the number of captured deer in each cell or throughout Gifu Prefecture. The thin lines mean the number of captured deer in each cell (scale on the left axis), and the thick line means the number of captured deer throughout the prefecture (scale on the right axis).

The relationships between the number of captured deer and each abundance index from monitoring, SPUE, pellet group survey, camera trapping survey, and block count survey are presented in Fig. 2. There was no clear pattern in the relationship between the number of captured deer and abundance indices. Compared with the number of captured deer, the relationship between abundance indices seemed positively correlated, although the relationship was weak.

Many parameters of the models with overdispersion terms in all observation models (i.e., models 1 and 3) did not converge ( of these parameters > 1.1). In contrast, all the model parameters considering overdispersion in observation models using the SPUE (i.e., models 0 and 2) and camera trap data (i.e., model 2) successfully converged ( of these parameters < 1.1). Furthermore, the median deer abundance in Gifu Prefecture was unrealistically high for models 1 and 3, and the 95% credible intervals of these models were relatively broad (Fig. 3), although these parameters did not converge. The median deer abundance in Gifu Prefecture in models 0 and 2 increased annually from 2008 to 2014 and then decreased (Fig. 3). In our four models, the total number of meshes as a state space from which sika deer abundance was estimated was 5184, of which 451 and 678 meshes (8.7% and 13.1%, respectively) exceeded densities of 100 deer/km² for models 1 and 3, respectively. In addition, 145 and 260 meshes (2.8% and 5.0%, respectively) exceeded densities of 200 deer/km², indicating extremely high values (see Discussion). Furthermore, the maximum estimated densities were 508.9 deer/km² and 706.8 deer/km² for models 1 and 3, respectively, which were extremely high and unrealistic. On the other hand, there were very few meshes with estimated densities exceeding 100 deer/km² for models 0 and 2 (two and seven meshes), and the maximum values were 112.1 deer/km² and 117.7 deer/km², respectively.

Here, we show the results of model 2, which provided realistic estimates for sika deer density. All parameters with a vague prior distribution converged and were identifiable in the model 2 (Appendix 1). Although the estimated abundance was generally correlated with all monitoring data, the plots in the scatter diagram showing the relationship between the estimated abundance and SPUE and camera trap data were scattered because the observation models of the SPUE and camera trap data considered overdispersion in model 2 (Fig. 4). The mean intrinsic population growth rate (logarithmic scale) was 0.075, and the 95-percentile interval ranged from 0.021 to 0.116. The estimated hunting ratios were spatially heterogeneous (Fig. 5), and the increase in deer abundance was suppressed by strong hunting pressure (Fig. 6). For example, there were two high-abundance regions in the center (around Gujo City and Gero City) and in the southwestern part (around Ogaki City and Ibigawa Town) of the Gifu Prefecture. The abundance in the central part decreased from 2014 because of the high hunting ratio, but that in the southwestern part did not decrease because of the relatively low hunting ratio. The results for model 0 were similar to those for model 2 (Appendix 2).

Fig. 2.

The relationship among each monitoring data in each cell in each year; the number of captured deer (Captured deer), sight per unit effort (SPUE) in hunters' reports, the number of pellet group by pellet group survey (Pellet group density), the number of deer found by a block count survey (Block count deer density), and the number of photographed deer by camera trap survey (Photographed deer). No pair of block count data and camera trap data in the same year and cell was obtained in this study period, so that the plot is omitted.

Fig. 3.

Estimated deer abundance of Gifu Prefecture. Left, estimates of all models; dark gray polygons indicate 95% credible interval of model 0 and 2, and light gray polygons indicate 95% credible interval of model 1 and 3. Right, estimates of models 0 and 3; the polygon with slanting lines indicates 95% credible interval of model 0, and the polygon with vertical lines indicates 95% credible interval of model 3. Model 0: An overdispersion was considered only in the observation model of SPUE, and there was no camera data; Model 1: An overdispersion was considered in all observation models, and there was no camera data; Model 2: An overdispersion was considered in the observation models of SPUE and camera trap survey, and there was camera data.; Model 3: An overdispersion was considered in all observation models, and there was camera data.

Fig. 4.

The relationship between the median of estimated abundance by model 2 and each monitoring data. The dotted line in the left-bottom panel indicates a diagonal line.

Fig. 5.

The posterior median of the hunting ratio (hrt,c) of model 2. Each cell indicates the hunting mesh (ca. 5 km × 5 km) of Gifu Prefecture.

Fig. 6.

The posterior median of the deer abundance of model 2. Each cell indicates the hunting mesh (ca. 5 km × 5 km) of Gifu Prefecture.

Discussion

The consideration of overdispersion terms for all observation models in models 1 and 3 caused an ill-convergence in the results. Although a previous HBM applied to Yamanashi Prefecture incorporated overdispersion terms into all observation models (Iijima et al. 2013), such an HBM could not be applied to Gifu Prefecture. The observational data other than SPUE in Gifu Prefecture were smaller than those of Iijima et al. (2013), both in terms of the number of years of survey conducted and the number of surveyed sites. Therefore, the estimates of models 1 and 3 might be strongly influenced by the trend of SPUE, which seemed to have lower precision data. As a result, population estimates could not be estimated appropriately. This study showed that we were able to successfully simplify the previous model to obtain estimates of HBM for sika deer management in a situation where we could not apply the previous model of Iijima et al. (2013) to currently available data. However, the risk that the simplification of the model may bias the estimates should be recognized. Although it is desirable to consider overdispersion in all observational models, as in Iijima et al. (2013), the models successfully estimated in this study did not take this structure, i.e., there is a possibility of overfitting observation data, except for SPUE. The model was simplified to obtain the most appropriate estimates possible from the currently available data in this study. However, this trial is not the best solution for a more appropriate estimation of long-term wildlife management. To estimate sika deer abundance with a more appropriate model structure, it would be desirable to improve SPUE data quality and obtain more observational data other than SPUE.

We constructed four models to estimate the population of sika deer inhabiting the Gifu Prefecture. Consequently, the estimates differed significantly from model to model. In the models (models 1 and 3) that considered overdispersion in all observation models, an unrealistically large number of deer was estimated, and none of the parameters converged. In previous studies, the maximum values of the carrying capacity of sika deer were reported to be 118 deer/km² by direct observation (Kaji et al. 2004) and 98.4 deer/km² by statistical modeling (Iijima and Ueno 2016). In models 1 and 3, 8.7% and 13.1% of the meshes exceeded 100 deer/km², and 2.8% and 5.0% exceeded 200 deer/km², respectively. Therefore, these values were difficult to determine realistically from previous studies. Furthermore, the maximum estimated densities were 508.9 deer/km² and 706.8 deer/km², which can be considered unrealistic. Therefore, we conclude that we could not obtain the appropriate estimation results using models 1 and 3. In contrast, models 0 and 2, which are observation models other than SPUE and camera trap data without overdispersion terms (i.e., strongly relying on the data from the pellet group survey and block count survey), produced realistic estimates for most meshes. There were very few meshes with estimated densities exceeding 100 deer/km² for models 0 and 2 (two and seven meshes, respectively), and the maximum values were 112.1 deer/km² and 117.7 deer/km², respectively. These values were realistic estimates that did not significantly exceed the maximum carrying capacity reported in previous studies. Furthermore, interannual trends in the estimated populations showed fluctuations that were consistent with some events in capturing sika deer in the prefecture. For example, there was an increasing trend in the estimated population until 2014, following which it began to decline in 2015. As shown in Fig. 1, the number of animals captured in Gifu Prefecture nearly doubled in that year compared to that in the previous year, with more than 17 000 sika deer captured. Although the downward trend in the population had strengthened in 2017, more than 17 000 deer were captured again that year. These facts suggest that the populations of sika deer in Gifu Prefecture estimated by models 0 and 2 are generally reasonable for understanding and examining the status of the current local population and the strategy of population adjustment for conservation and management.

In this study, the mean intrinsic population growth rate (logarithmic scale) was 0.075. This value was much lower than the value of 0.253 reported in a previous study using a similar model, in Yamanashi Prefecture during 2005 and 2010 (Iijima et al. 2013). There are two possible reasons for this. The first is the difference in the phase of the sika deer population in previous studies: the case in Yamanashi Prefecture (Iijima et al. 2013) was in an expansion phase in which the number of sika deer continued to increase yearly. On the other hand, according to our estimation, the abundance of sika deer in Gifu Prefecture has been decreasing since 2014, suggesting that the sika deer population phase may differ from that in Yamanashi Prefecture. However, Kaji et al. (2004) estimated the growth rates of sika deer populations on Nakanoshima Island and Cape Shiretoko in Hokkaido Prefecture to be 0.152 and 0.190, respectively, during a period when density-dependent effects might have affected the sika deer abundance. These values were larger than those estimated in this study. Therefore, it should be noted that there is a possibility that the growth rate estimated in our study was underestimated. Secondly, during the study period of 2008–2019, there were three snowy winter seasons: January 2014, December 2014–February 2015, and January 2017 (Takayama Station, Japan Meteorological Agency,  https://www.data.jma.go.jp/obd/stats/etrn/index.php?prec_no=52&block_no=47617, Accessed 29 Sep 2022). In particular, the winter of 2014 was one of the snowiest years and caused heavy snowfall damage to forests, houses, and infrastructure in Gifu Prefecture (Yamashita et al. 2016). In the spring, following these heavy snowfalls, local hunters found many deer carcasses in the mountains (Yusuke Seto, personal communication). As a result of multiple weather events occurring during the study period that might significantly impact survival rates, lower growth rates might have been estimated because of the impact of years in which actual growth rates were significantly lower.

The results of this study clearly show that changing the model structure and data can significantly affect the estimates. Although it is difficult to determine which model is correct, when there is no certainty about the mechanism or model structure in nature, it is desirable to conduct an estimation with multiple models, observe the estimation results and trends, and make decisions using the estimated results that seem optimal. This study indicates that in planning for appropriate population management based on estimates, the person in charge must have a good understanding of the structure of the model and the meaning and reliability of the estimated values. In this study, the model did not converge when data obtained from some scientific surveys were treated in the same manner as the SPUE. However, by setting up the structure of the model appropriately and considering the quality and quantity of the data, good estimation results were obtained. However, in the absence of such a comparative study and verification of the validity of the estimates, there is a risk that the persons in charge may not be able to set the target number of animals to be captured or may set incorrect targets.

It is desirable to implement wildlife protection and management methods based on the local populations of target animals. Even if the estimation results of the HBM, which considers the prefecture as one space, are available, they do not necessarily provide helpful information for protection and management. The size of the home range of sika deer was limited to a small spatial scale (50 ha, Miyashita et al. 2008; 19–217 ha, Takii et al. 2012), which would have caused the spatial difference in deer abundance in 5 km square grid cells in this study (Fig. 6). Furthermore, the number of hunted deer differed significantly locally (Fig. 1). These facts highlight the importance of local population management. This point of view is essential in regions such as Gifu Prefecture, which has a large area, significant differences in vegetation and climate between the north and south, and multiple local populations of sika deer. One of the achievements of this study is the estimation of the sika deer population with spatial heterogeneity in Gifu Prefecture. As a result, it became clear that there were significant differences in the local population trends. In particular, significant differences in trends were observed at sites with high estimated densities and catches (central and southwestern prefectures). Such spatial heterogeneity in local trends cannot be represented by the HBM, which is constructed by homogenizing the county area and averaging the data. Notably, the Ministry of the Environment conducted a spatial population estimation of sika deer in the Kanto region (Ministry of the Environment 2021,  https://www.env.go.jp/press/109597.html, Accessed 29 Sep 2022). Although this result appeared to take into account the spatial variation of the population, the total estimated population was from an HBM without a spatial structure. The number of pellet groups obtained in each region was used as an index to prorate the total estimated population and express the spatial density distribution. In this case, the HBM assumed that the population trend was the same for the entire area in the model and did not consider the trend of increasing or decreasing population in the local area. A similar method was used in estimating deer population in Tokyo Prefecture (Tokyo Prefecture 2022), which defines a fairly large area as the state space, so the process of population dynamics on spatially fine scales cannot be considered. Therefore, when understanding the spatial density distribution, it should be noted that the above-mentioned method of prorating the population for the whole state space using the local density index is entirely different from the spatially fine-scale HBM applied in this study.

In Gifu Prefecture, surveys on the decline of understory vegetation have been conducted several times in entire prefecture (Tsunoda et al. 2017; Gifu Prefecture 2021). The survey results indicated that the decline of the understory vegetation in the southwest region was significant in recent years compared to that in the past (Gifu Prefecture 2021). In this area, an adequate decrease in population has not been observed since the increase in population estimated by the results of this study (Fig. 6), which indicates that deer living at high densities in this area have decimated the understory vegetation. Even though trends of the understory vegetation decline and the results of this study were obtained from completely different data, the trends of the results are consistent, which indicates that the spatial estimation results of this study are reasonable. Therefore, a spatially fine-scale HBM, such as the one applied in this study, can contribute to policymaking decisions by providing an overview of the effects of the sika deer population on such local situations.

In this study, we constructed a spatially fine-scale HBM to estimate the population of sika deer in the Gifu Prefecture. We obtained reasonable estimates of the sika deer population using different local trends, thereby highlighting heterogeneous demographic events related to the sika deer population in the prefecture. The estimates obtained were spatially valid and consistent with the trends in the decline of understory vegetation in each region. Estimates obtained from the spatially fine-scale HBM are more helpful than those obtained using the entire prefecture as a single state space, thus enabling a more detailed understanding of the situation and guiding future actions. Our findings also suggest that reliable observation data as abundance indices might be needed for estimating wildlife abundance with the HBM, even if these indices were sparse over space and time. Such observation data provide a useful reference when estimating the parameters of the HBM in the face of spatiotemporal uncertainty. To manage the population of sika deer more realistically in the future, spatially fine-scale HBM should be applied, as data can be used adequately without factoring in estimations, while taking into account the spatial heterogeneity of the sika deer population. Each local government should accumulate the necessary data, understand the quality of the data, and develop an infrastructure for building a spatially fine-scale HBM. In addition, we could not construct a “spatially explicit model” that considers interactions between meshes (e.g., movements of Japanese deer) in the process model. This was due to limitations in the amount and type of data as well as in the capacity of computational resources. If sufficient data and computational power are available, estimating the sika deer population would be desirable by constructing a “spatially explicit model” that makes it possible to design the process model closer to reality.