Land abandonment and changes in snow cover period accelerate range expansions of sika deer

Abstract Ongoing climate change and land‐use change have the potential to substantially alter the distribution of large herbivores. This may result in drastic changes in ecosystems by changing plant–herbivore interactions. Here, we developed a model explaining sika deer persistence and colonization between 25 years in terms of neighborhood occupancy and habitat suitability. We used climatic, land‐use, and topographic variables to calculate the habitat suitability and evaluated the contributions of the variables to past range changes of sika deer. We used this model to predict the changes in the range of sika deer over the next 100 years under four scenario groups with the combination of land‐use change and climate change. Our results showed that both climate change and land‐use change had affected the range of sika deer in the past 25 years. Habitat suitability increased in northern or mountainous regions, which account for 71.6% of Japan, in line with a decrease in the snow cover period. Habitat suitability decreased in suburban areas, which account for 28.4% of Japan, corresponding to land‐use changes related to urbanization. In the next 100 years, the decrease in snow cover period and the increase in land abandonment were predicted to accelerate the range expansion of sika deer. Comparison of these two driving factors revealed that climate change will contribute more to range expansion, particularly from the 2070s onward. In scenarios that assumed the influence of both climate change and land‐use change, the total sika deer range increased by between +4.6% and +11.9% from the baseline scenario. Climate change and land‐use change will require additional efforts for future management of sika deer, particularly in the long term.

Climate change affects the distributions, population dynamics, and body condition of many species (Drever et al., 2012;Kausrud et al., 2008), including large herbivores (Büntgen et al., 2014). Some of the effects can emerge from the direct effects, such as thermoregulation or water limitation, and the indirect effects, such as access to forage or vegetation productivity and quality (Mysterud & Saether, 2011). In the high-latitudinal regions, large herbivores often enter negative energy balance in winter because of the deep snow, cold temperatures, and limited forage intake. Ongoing climate change will lead to warmer and shorter winters with less snow in these regions (IPCC, 2013). This would affect the survival and reproduction rates of herbivores by decreasing the energy loss during winter and increasing the length of the growing season for vegetation (i.e., food resources).

Department of Economic and Social Affairs, Population Division 2015).
This human depopulation causes land abandonment, particularly in rural areas. These abandoned fields provide additional habitats for wildlife (Bowen, McAlpine, House, & Smith, 2007;Moreira & Russo, 2007) and are now among the major causes of expansion of the distributions of large herbivores, particularly in developed countries.
The consequences of these complex interactions may not be straightforward and depend on the species and communities. Although there is growing recognition that the effect of altered species interactions sometimes exceeds the direct effect of change in environmental drivers (e.g., Ockendon et al., 2014), this perspective has until recently been overlooked.
Japan is one of the world's snowiest regions (Yano, 1993), particularly in the mountainous region along northwest coast. Because Japan's mountainous areas are located at relatively low latitudes (35-45°N) in a warm temperate zone, even slight climate changes are likely to affect snow conditions (Kawase et al., 2013;Yamaguchi, Abe, Nakai, & Sato, 2011). These changes may cause significant changes in snow-dominated ecosystem in Japan (Kudo, Amagai, Hoshino, & Kaneko, 2011;Tsuyama et al., 2011Tsuyama et al., , 2012. Moreover, snow is considered important for survival of the sika deer (Cervus nippon Temminck, Figure 1), a large herbivore that is causing serious negative impacts on vegetation and ecosystems throughout Japan (Ohashi et al., 2014;Takatsuki, 2009a). Until the 1970s, the distribution of the deer was biased toward the less-snowy areas with a maximum snow depth of less than about 50 cm (Takatsuki, 1992a). In the 1990s, expansions of the wintering ranges of the deer into snowy areas were reported, and this is suspected to have been facilitated by climate change (Li, Maruyama, Koganezawa, & Kanzaki, 1996).
On the other hand, several reports have emphasized the importance of the effect of land-use change on the recent range expansion of sika deer (e.g., Agetsuma, 2007;Miyashita et al., 2008). Sika deer prefers to live in landscapes with mixed grassland and forest. Logging increases the availability of forage for sika deer because of the improved light conditions, and these changes last for several decades after logging until the forest regrowth (Takatsuki, 1992b). In addition, a shift in Japan's industrial structure since the 1960s, together with the energy revolution, has caused depopulation of rural areas, and the land abandoned in these areas may have provided forage and shelter for sika deer. The country is now facing a rapid decline in human population, which is leading to the further abandonment of land in rural areas, creating an opportunity for sika deer's range expansion.
Climate change, together with the abandonment of land in rural areas, will likely add to the population growth and range expansions of sika deer, and accelerate the changes in ecosystem function. To minimize the negative effect on ecosystems, it is critical that we develop plans for the strategic management of large herbivores, and these plans are necessary to incorporate the predicted future changes in both climate and land use. Therefore, the assessment of the impacts of climate change and land-use changes on sika deer populations is essential.
Our first objective here was to statistically evaluate whether the two factors (climate change and land-use change) were related to the change in range of sika deer in previous decades. If this was the case, it was then important to determine which factor had contributed more F I G U R E 1 An adult female sika deer (Cervus nippon Temminck) and her cub to the changes. Our second objective was to use statistical models to predict the future range changes of sika deer under climate change and land-use change scenarios.

| Target species
Until the late 19th century, sika deer had been widely distributed across lowlands in Japan. In the 20th century, overexploitation and habitat loss caused by human activity reduced the deer's population and range. Since then, the deer's distribution has been skewed toward mountainous areas with little human activity. After the establishment of legal protection acts in the 1970s, the deer's population recovered, and its geographic range expanded by 170% in 25 years (Ministry of the Environment, 1979Environment, , 2004 Figure S1). The animal's negative impacts on forestry, agriculture, and ecosystems have now become problematic (Ohashi et al., 2014;Takatsuki, 2009a). Accordingly, in 2013, the Japanese government decided on a national goal to halve the sika deer population. In addition, in 2014, the Wildlife Protection and Hunting Management Law was revised (Ministry of the Environment, http://www.env.go.jp/nature/choju/index.html), and the government plan was to provide the financial supports to local governments to reduce the sika deer.
Since the 1990s, there were vigorous debates about the likely contributors to the recent range expansion of sika deer. Lack of predators (wild wolves-natural predators of these deer-became extinct in Japan 100 years ago), aging and a decrease in the numbers of hunters, climate change, and land-use change were the major suspected factors. Several studies have used a species distribution modeling approach to examine the determinant factors of sika deer distribution. In an analysis in Hokkaido, in northern Japan, Kaji, Miyaki, Saitoh, Ono, and Kaneko (2000) pointed out the importance of the effects of snow depth and forage quality. Later, Okumura, Shimizu, and Omasa (2009), in an analysis covering the other three main islands of Japan, revealed the importance of the effects of snow and human land use. However, a recent analysis that covered the whole of Japan found that not climate factors, but anthropogenic factors affected sika deer distribution (Saito, Momose, Inoue, Kurashima, & Matsuda, 2016). Although the last two of these studies considered the effects of the dispersal ability of sika deer by including their distance from their 1970s range, they did not consider the temporal changes in habitat conditions. Accordingly, the importance of temporal changes in habitat condition is still unclear.

| Dataset of sika deer distribution
Here, we used the presence-absence data collected across Japan for sika deer in 1978 and 2003 (Ministry of the Environment, 1979Environment, , 2004. These surveys were conducted by using questionnaires, interviews, and field survey. The presence or absence of sika deer was recorded in 3′45″ × 2′30″ grid cells (approximately 5 × 5 km 2 , hereafter called "5-km grid cells"). These data were collected regardless of the season. Respondents to the questionnaires and interviews were forest managers, forest researchers, park rangers, hunters, and members of forest owner cooperatives (Saito et al., 2016). Along with these surveys, an additional survey was conducted after primary data aggregation to prevent an imperfect detection (Ministry of the Environment, 1979Environment, , 2004. Therefore, we assumed that the probability of imperfect detection in these surveys was sufficiently small to prevent large biases in the model output. We used data from 15,256 grid cells on the four main islands of Japan (Hokkaido, Honshu, Shikoku, and Kyushu: Figure S1); the total coverage was about 351.8 × 10 3 km 2 . In 1978, the number of grid cells in which sika deer were recorded as present was 4079, and the number in which they were recorded as absent was

| Climatic variables
We used the snow cover period and the maximum snow depth as candidate climatic variables to explain sika deer distribution. We estimated the snow cover period and the maximum snow depth on 45″ × 30″ grid cells (approximately 1 × 1 km 2 ; hereafter referred to as 1-km grid cells) by using simple degree day method (Kominami, Tanaka, Endo, & Niwano, 2005). Firstly, a daily mean temperature and daily precipitation were calculated in 1-km grid cells by distance-weighted mean of three nearest Automated Meteorological Data Acquisition System (AMeDAS) station managed by Japan Meteorological Agency (data were derived from Agriculture, Forestry, and Fisheries Basic Numeric DataBase of AFFRIT, MAFF, Japan). Snow melt coefficient was optimized using satellite data (SPOT Vegetation) (Asaoka & Kominami, 2012. Daily SWE (snow water equivalent) and then existence of snow cover and snow depth were calculated in 1-km grid cells; we calculated the mean value in 5-km grid cells for the analysis, and the mean values between either 1973 and 1977 or 1998 and 2002 were used. Because the Hokkaido subspecies of the deer (Cervus nippon yesoensis), which has a large body, seems to have greater tolerance to deep snow (Tokida, Maruyama, Ito, Furubayashi, & Abe, 1980), regional differences between Hokkaido and the other three islands were considered as an interaction term between the climatic variables.  1976, 1987, 1991, 1997, and 2006 Information, Ministry of Land, Infrastructure, Transport and Tourism, http://nlftp.mlit.go.jp/ksj-e/gml/datalist/KsjTmplt-L03-a.html.). We used the inverse distance weight method for interpolation. These data were originally calculated in 1-km grid cells; we calculated the mean values in 5-km grid cells for the analysis.

| Topographic variables
For other environmental factors, we used the mean slope inclination to consider topographic accessibility, which could affect sika deer distribution by influencing the frequency of human activity (Takii, Izumiyama, Mochizuki, Okumura, & Sato, 2012). These data were originally calculated in 1-km grid cells; we calculated the mean values in 5-km grid cells for analysis.

| Model construction
We analyzed whether sika deer changed its range according to the estimated changes in suitability surfaces ( Figure 2). Sika deer occupancy in 2003 was assumed to be the result of the processes of persistence (a previously occupied cell may still be occupied) and colonization (a previously unoccupied cell may be newly occupied).
The occupancy probability of a cell in 2003 was defined as conditional upon the state of occupancy of the cell in 1978, as in the other dynamic occupancy models (Bled, Nichols, & Altwegg, 2013;Yackulic et al., 2012).
If we consider i = 1, 2, …, N spatial units (i.e., grid cells), the occupancy X i,yr (X i,yr = 1 for presence, 0 for absence) in cell i during the survey year (yr), by a Bernoulli distribution with parameter q i,yr can be denoted as: The probability of occurrence of sika deer in 2003 was derived from the occupancy status in 1978 and by a dynamic process of colonization and extinction as follows: where φ i and γ i are the persistence and colonization probabilities, respectively, for cell i between 1978 and 2003. We defined these probabilities by using two components: (1) the probability of immigration from neighboring cells during one time period (probability of dispersal, D i ) and (2) the species' preference for the cell's environment (habitat suitability, H i ) (Lawes, Mealin, & Piper, 2000). Therefore, we defined the persistence (φ i ) and colonization (γ i ) probabilities as: where φ 0 and γ 0 are intercepts, φ d and γ d are slopes for the probability of dispersal D i , and φ h and γ h are slopes for habitat suitability H i .
We assumed that the value of D i was determined by neighborhood occupancy during the previous time period. A cell that has a large number of occupied neighbors is more likely to stay occupied, or to be colonized, than to become unoccupied. We tested the combinations of five patterns of neighborhood radius (10, 25, 50, 75, and 100 km) and three patterns of assignment of weight to grid cells within the neighborhood (equal weight, weighted by inverse distance, and weighted by square of inverse distance), which gave 15 candidate neighborhood occupancy values. Although there are several other methods to separate the dispersal from habitat effects (Broms, Hooten, Johnson, Altwegg, & Conquest, 2016;Sutherland, Elston, & Lambin, 2012), this method has the advantage of low computational intensity (Dormann et al., 2007) and allows the application of the model to a large dataset that covers a large spatial extent with fine resolution.
To obtain an optimal habitat suitability surface to explain the change in range of sika deer between 1978 and 2003, we used the distribution data for sika deer and environmental data in 1978 to train the multiple candidate models to obtain the habitat suitability surfaces: where α is intercept and β j are coefficient parameter of jth environmental variable, x j. These candidate models include all possible combinations of climatic variables, land-use variables, and topographic variables. To include the nonlinear relationship between sika deer occurrence and habitat conditions, we checked the generalized additive model fitting curves and decided to include a second-order term for the proportion of agricultural land. We also compared the predictive performances of two incompatible climatic variables, namely the snow cover period and the maximum snow depth. From 16 candidate combinations of all variables, we excluded the model with both the snow cover period and the maximum snow depth to avoid the effect of multicollinearity and nominated 12 candidate habitat models for further analysis (Table S1).
Because the presence and absence of sika deer may not be at equilibrium as a result of dispersal limitation, we also considered the models with autocovariates. We considered two patterns of neighborhood radius (10 and 25 km) and three patterns of assignment of weight to grid cells within the neighborhood for autocovariates, which resulted in 84 candidate models.
Subsequently, we used the environmental variables in 1978 and 2003 to calculate a habitat suitability surface for 84 habitat suitability models in each survey period, namely habitat suitability surface in 1978 (H 1978,i ) and 2003 (H 2003,i ). A grid cell with a relatively high suitability value was likely to stay occupied or to be preferentially colonized by sika deer. We defined habitat suitabilities as the log odds (probability on a logit scale) of predictions from the fitted models: Combination of the two survey periods and 84 habitat suitability models resulted in 168 candidate habitat suitability surfaces.
By using the program WinBUGS (Lunn, Thomas, Best, & Spiegelhalter, 2000), we used Bayesian modeling and Markov chain Monte Carlo posterior sampling to estimate the parameters of the dynamic-range change model to explain the sika deer presence in 2003. We ran three chains, using noninformative priors, for 1,000 iterations after a 100-iteration burn-in period. Overall, we compared the Watanabe-Akaike Information Criteria (WAIC) (Watanabe, 2010) of 2,520 competing models and selected the model with the minimum WAIC value, which indicated the best combination of neighborhood occupancy and habitat suitability to explain the range dynamics be-

| Future climate change and land-use change scenario for the future
For future predictions, we used four environmental scenario groups: (1) baseline, (2) land-use change only, (3) climate change only, and (4) land-use change and climate change, from 2004 to 2103 (Table S2).
Baseline is a scenario of no change in climate and land use.
Land-use change only scenario group was made from nine landuse change scenarios (Hanasaki, Takahashi, & Hijioka, 2012) based on different patterns of future human population projection. The human population projections were made by combining three patterns of birth rate and death rate-high birth rate and low death rate, medium birth rate and medium death rate, and low birth rate and high death rate-and three patterns of spatial distribution pattern, namely centralization, baseline, and decentralization (Ariga & Matsuhashi, 2012

| Drivers of change in the range of sika deer in the last 25 years
In the top five models with the lowest WAIC values, we used landuse variables, climatic variables, and topographic variables to calculate habitat suitability, excluding any autocovariate (Table 1). Within the climatic variables, the snow cover period was selected in the top five models, but the maximum snow depth was not selected in any of these models. In addition, when comparing the models with same distance and weights, the models that used habitat suitability surface Within the parameters of the range change model, the coefficients for probability of dispersal and habitat suitability were positive in both the persistence model and the colonization model ( Table 2).
The coefficient for the probability of dispersal was larger in the colonization model than in the persistence model, whereas the coefficient of habitat suitability did not differ markedly between persistence and colonization.
In the best habitat suitability model, sika deer showed a negative response to snow cover period; this was less significant (i.e., more weakly negatively correlated) in Honshu, Shikoku, and Kyushu, as indicated by the interaction term (Table 3). Sika deer showed positive responses to slope inclination, forest cover, and wasteland cover and negative responses to built-up area (Table 3). Unimodal responses to agricultural land were also confirmed.
From 1978 to 2003, habitat suitability for sika deer increased in 10,925 grid cells, accounting for 71.6% of the total area of Japan, mainly in the northern and mountainous regions (Table 4, Figure 3a-c).
In contrast, habitat suitability decreased in 4,331 grid cells, accounting for 28.4% of the total area, mainly in suburban areas (Table 4,

| Predicted future range changes of sika deer in the next 100 years
The result of our simulation of future changes showed the considerable range expansion of sika deer in the coming 100 years. In 2028, the effects of climate change and land-use change on deer range would be relatively small ( Figure 4). However, the differences among the four groups of future climate change and land-use change scenarios would become apparent from 2078 onward. In the baseline T A B L E 4 Trend of change in habitat suitability between 1978 and 2003 (the numbers of cells in which it increased or decreased), and the relative contribution of land-use change and climate change (LU > CL: the absolute value of the changes in habitat suitability caused by land-use change was larger than that caused by climate change; LU < CL: the absolute value of the change in habitat suitability caused by climate change was larger than that caused by land-use change). χ 2 -value = 12,697, df = 1, p-value < .001 Habitat suitability (no. of grid cells)

Increased Decreased
Absolute value of difference in contribution LU > CL 220 4011 LU < CL 10,705 320 Total 10,925 4331 scenario, sika deer would expand its range by 304.2 × 10 3 km 2 in 2103 ( Figure 5a). In the "land-use change only" scenario, the total range would increase by 308.5 × 10 3 to 312.1 × 10 3 km 2 , a + 1.4% to +2.6% increase over the baseline scenario ( Figure 5b). Compared with the baseline scenario, the deer range would expand around suburban areas (Figure 5e). In the "climate change only" scenario, the total sika deer range would increase by 313.7 × 10 3 to 334.0 × 10 3 km 2 , resulting in a + 3.1% to +9.8% increase over the baseline scenario ( Figure 5c). Compared with the baseline scenario, the deer range would expand in northern Honshu (Figure 5f). In the "land-use and climate change" scenario, the total sika deer range would increase by 318.3 × 10 3 to 340.4 × 10 3 km 2 , which result in +4.6% to +11.9% increase over baseline scenario (Figure 5d). Range expansion was predicted to occur both in areas surrounding the suburbs and in northern Honshu (Figure 5g).

| DISCUSSION
4.1 | Which of these two factors-land-use change or climate change-has contributed more to the range change of sika deer in recent decades?
We found here that habitat suitability of sika deer was explained by a combination of climatic variables, land-use variables, and topographic variables (Table 1). Our results revealed the importance of the snow cover period compared with the maximum snow depth, which had been used in previous studies (Kaji et al., 2000;Okumura et al., 2009;Saito et al., 2016). In snowy regions in Japan, evergreen dwarf bamboo species dominate on the forest floor (Tsuyama et al., 2011(Tsuyama et al., , 2012, and snow cover prevents the access to this forage. Therefore, the snow cover period directly affects the amount of forage available in winter (Minamino & Akashi, 2011). This indirect effect of climate might be a more important determinant of the distribution of sika deer than are direct effects of climate such as the energy cost of moving in deep snow (Parker & Robbins, 1984).
The response of sika deer distribution to land-use variables (Table 3) was consistent with the other previous studies. The positive responses to forest and wasteland and the unimodal response to agricultural land may indicate the preference for forest, grassland, and forest-agriculture ecotones (Saito et al., 2016). Sika deer tend to avoid built-up areas, owing to the frequent human disturbance (Okumura et al., 2009;Saito et al., 2016). In our study, slope inclination was selected as an important variable. In general, deer avoid using steep slopes to minimize unnecessary energy loss (Ganskopp & Vavra, 1987;Moen, 1976). However, in the mountainous areas of Japan, deer find steep slope safe, with little human disturbances, during the hunting season (Takii et al., 2012). In this regard, our result agrees with those of this earlier study and may reflect the indirect effect of human disturbance, which had been overlooked in previous macroscale studies.
We found here that both climate change and land-use change between 1978 and 2003 have affected the changes in deer range ( Figure 3). Because sika deer have high mobility, they can change their distribution correspondingly to the spatiotemporal change in habitat condition, rather than stick to previous habitat condition. Our result indicated that a decrease in snow cover period in the last few decades has increased the habitat suitability of sika deer throughout most of the Japanese islands-especially in northwestern Japan and in mountainous regions. In the snowy area, sika deer migrate between separate summer and winter ranges to avoid unsuitable climatic condition, while it becomes sedentary in the area with no snow (Takatsuki, 2009b;Yabe & Takatsuki, 2009). Under a favorable climate, there is no need for seasonal migrants to return to their winter ranges, so they tend to stay in their summer ranges. These individuals become the founders of new populations, and this leads to further range expansion (Maruyama, 1981). The short snow cover period in recent years may have contributed to the expansion of populations in habitats that were previously unsuitable because of climatic conditions.
In contrast, changes in land use resulted in a decrease in habitat suitability of sika deer at suburban area and in a slight increase in mountainous area. Decrease in habitat suitability in suburban area was mainly due to the decrease in forest and the increase in built-up area.
Slight increase in habitat suitability in mountainous area was mainly due to the decrease in wasteland and the increase in forest.
Our results also revealed that the response to changes in snow cover period differed between Hokkaido and the other islands.
Habitat suitability was greater in Hokkaido than on the other three islands (Table 3). Additionally, we found a more significant effect of the snow cover period in Hokkaido, as indicated by the interaction term (SCP × RGN in Table 3). Both ecological and historical differences between two regions are expected to cause the difference in the sensitivity of populations to changes in the snow cover period. As ecological factor, Hokkaido population has the greater tolerance to deeper snow due to its large body size, compared with the population in other three islands. Release from suppression by a long snow cover period might result in explosive range expansion as a reaction, as we have observed in past decades (Kaji et al., 2000).  (Dawe et al., 2014). However, process and mechanism to cause the difference between two region, and the relative importance of these two factors, are still unknown. More detailed research is required to clarify the factors responsible for the population differences between Hokkaido and the other islands.

| Effects of future land-use change and climate change on the range of sika deer
We found here that the effects of climate change and land-use change were minimal in 2028, although we had predicted that the sika deer range would expand (Figures 4 and 5). However, the effects of landuse change and climate change became more substantial from 2078 onward (Figures 4 and 5). Both land abandonment caused by human depopulation and the decrease in snow cover period caused by climate change had predicted to accelerate the deer's range expansion, and combining the two drivers gave the largest rate of range expansion. In one study in Europe (Rondinini & Visconti, 2015), populations of red deer (Cervus elaphus Linnaeus)-a species that is closely related to sika deer and is distributed over the Eurasian continent-were also predicted to increase under the predicted climate change and land-use change, unlike the case with other large herbivores, such as moose (Alces alces Linnaeus) and reindeer (Rangifer tarandus Linnaeus), which inhabits in more northern biome. These differences in response to climate change relate to tolerance to heat condition. Previously, the southern limit of sika deer was at northern Vietnam in the tropical environment, which has now become extinct in the wild (McCullough, 2009). Sika deer have relatively high tolerance to high temperature and may receive the positive effect after climate change, at least in Japan.
There are many predictive studies that have reported the risks of species range loss and extinctions as a result of climate change and land-use change (Rondinini & Visconti, 2015). However, most of these works only consider the direct effect of change in driver and fail to notice the indirect effect of change in biotic interactions, such as the intensification of deer impact on vegetation. For example, Siebold's beech (Fagus crenata Blume)-the dominant tree in the cool temperate zone of Japan-is projected by some models to decrease its potential habitats widely according to the future climate change (Koide et al., 2016;Matsui et al., 2004Matsui et al., , 2009Nakao et al., 2013). Additionally, the regeneration of Fagus crenata would become suppressed under heavy browsing by sika deer (Takatsuki & Gorai, 1994 ing on their vulnerability to deer browsing. To predict the future vegetation, we need to take account for the synergistic effect of changes in climate, land-use, and intense impact of deer browsing on vegetation.
Our results revealed that the effects of climate change and landuse change may strengthen after the 2050s. Although the effects of both drivers are not likely to result in short-term consequences, if we overlook their effects, we may underestimate the effort required to control deer populations, particularly over the long term. The sika deer range in northern Honshu and in mountainous areas may particularly be strongly influenced by climate change in the next 100 years.
Building our capacity for monitoring and population control in these areas may help with early detection of range expansion in the front line and thus with the rapid population control.