Prediction of the visit and occupy of the sika deer (Cervus nippon) during the summer season using a virtual ecological approach

Study area and camera trapping system

The study area included the northern region of Tochigi Prefecture, Japan (Fig. 2). In Tochigi Prefecture, 54.4% of the land was covered by forest, 19.1% was covered by agricultural land in 2019 (Tochigi Prefecture 2021, https://www.pref.tochigi.lg.jp/a03/documents/keikakusho2267.pdf, accessed on Feb. 10, 2023). The northern region of Tochigi Prefecture has a relatively large area of forest. This area was the home range of the highest density of sika deer in Tochigi Prefecture in 2021. The camera trapping system consisted of 14 cameras (model no. 6210; Ltl-Acorn, Des Moines, IA, USA) that were placed in late April 2018 at 12 sites within the forest interior with two camera sets, namely ID 10–11, and ID 12–13 in neighboring areas (Fig. 2). The 12 sites spanned 84 km from west to east and 39 km from north to south (Fig. 2). The elevation of the sites ranged from 349 to 1033 m. The cameras were set horizontally at 50 cm above the ground and were operated until late November 2018. The cameras were checked every 1 or 2 months and the batteries and memory cards were replaced when necessary. Movements of the sika deer were reordered monthly from May to November. The month of April was excluded because the cameras were placed in late April. The virtual ecological model required the presence/absence of records for validation (described below), thus the number of deer captured in the photos was not considered. Finally, the visit and occupy of sika deer were recorded at 14 sites each month.

A grid size of approximately 1 km (termed “1-km mesh” hereafter) was used a as the study unit (Fig. 2). The 1-km mesh grid system is a standard Japanese unit used for several types of statistics (https://www.stat.go.jp/english/data/mesh/02.html, accessed on Feb. 10, 2023). To determine the appropriate number of 1-km mesh grids for the simulation study, a 10-km mesh grid, which is the high-order standard Japanese unit (i.e., one 10-km mesh includes 100 1-km meshes), was divided into the minimum number of areas to cover all 14 camera sites as the simulation target area to avoid arbitrary (Fig. 2). Finally, 4200 1-km mesh areas were included for the simulation (Fig. 2).

Virtual ecological model

A simple cellular automaton (CA) model can predict the visit and occupy of a target species based on candidate habitats in consideration of the proximity to food resources³². The grid was set to the same size as the unit of the predicted ranges. The model yields a theoretical number of visits (described below) to each cell, which serves as an area preference of the target species. Each cell has two parameters: cell identification (ID) and movement path vector (Fig. 3a). The cell ID indicates the spatial location of the cell within the study area. The movement path involves four variables representing the four directional vectors into adjacent cells (i.e., top, left, bottom, and right) (Fig. 3b). Each variable is a probability value (i.e., 0 to 1) independent of the other three variables that indicates the probability of movement success to the adjacent cells. In this study, the probability value was based on the proximity to availability food resources.

A group of sika deer was used as the unit for analysis. The model simulates the capability of movement within the target area. Thus, if a virtual population visited a neighboring cell, the number of visits to the cell is increased without disappearance of the starting cell. The virtual population moves in accordance with the movement probability values.

Movement probability between cells

The term “movement probability” is defined as the probability of movement success into an adjacent cell to the top, left, bottom, or right (Fig. 3b) with four probability values:

where m1, m2, m3, and m4 indicate the probability of movement success into the top, left, bottom, and right cells, respectively (Fig. 3b). Since these values are independent of one another, the maximum and minimum sums of m1, m2, m3, and m4 are theoretically 4 and 0, respectively. If all probability of movement success values are 0, the sika deer population in this cell cannot move to any other cell. Moreover, if all probability of movement success values are 1, the population in the cell can move to all adjacent cells.

The amount of food resources of deer was acquired from remote sensing measurements^35,36. Thus, two variables were used to represent food resource availability: the kernel normalized difference vegetation index (kNDVI)⁴¹ and landscape structure (Supplementary Fig. 1).

The kNDVI uses remote sensing measurements to assess the components of green vegetation⁴¹. As compared to the ordinal NDVI, which is the most widely used index of the condition of vegetation on terrestrial surfaces, the kNDVI has greater resistance to saturation, bias, and complex phenological cycles, and exhibits enhanced robustness to noise and stability across spatial and temporal scales⁴¹. The kNDVI appropriately represents the condition of vegetation to reflect the food resource availability for sika deer. The kNDVI was analyzed from the atmospherically corrected surface reflectance observed with the Landsat 8 Operational Land Imager and Thermal Infrared Sensor instruments at approximately 16-day intervals with a spatial resolution of 30 m (data collected in 2018). The mean kNDVI was calculated monthly for each 1-km mesh within the study area. The probability values (m1, m2, m3, and m4) were defined as the proximity to available food resources in a destination cell divided by the maximum value of the target area as relative values throughout the study area. These values reflect the spatial positions of the available food resources in the study area. If the food resources are continuously available, then the sika deer population tend to visit and occupy linearly.

The landscape structure is defined as a mixture of forests and grasslands because previous studies suggest that the forest edge has high availability of food resources for sika deer^37,38,42,43. The dataset was generated from a current vegetation map that classified the dominant plant species provided by the Biodiversity Center of Japan (Ministry of the Environment, https://www.biodic.go.jp/index_e.html, accessed on Feb. 10, 2023). The types of vegetation of the forests and grasslands were retrieved from the literature, then the original vegetation classes were re-classified⁴⁴ and overlayed on the 1-km mesh map. In this study, agricultural land types were classified as grassland. For a mesh with both forests and grasslands, the probability of movement was assigned a value of 1, while a mesh with either a forest or grassland was assigned a value of 0.5, because to treat these 2 components fairly. Every mesh of the study area included either a forest or grassland.

Movement simulation

First, simulations were conducted using two independent variables: kNDVI and landscape structure. Each simulation was initiated from one cell with the month, which is referred to as a “trial.” One step is defined as one day, thus the trial conducted in May consisted of 31 steps. A previous study reported that sika deer can travel about 50 km every 2 weeks³⁴. Thus, one step (movement of 1 km) in one day was considered a reasonable distance. Each trial was repeated for all cells i.e., all cells was used as the starting cell of “trial”. The sum of all trials is termed a “run.” Thus, each “run” consisted of n trials, where n is the number of cells in the CA field. In this study, there were 4200 cells. At each step, each attempt to visit a neighboring cell (top, left, bottom, and right) was based on movement probabilities. For each successful movement, the presence/absence value assigned to the cell was increased from 0 to 1, i.e., change from absence to presence. The next step was then initiated from any newly visited cell and the previously visited cells. Cells with high values indicated the possibility of visitation by a virtual population from several other cells. The assigned value was used as a metric of the preference of the visited cell. In this study, 100 runs were conducted each month from May to November.

Second, simulations were conducted using a combination of movement-related variables with two types of combination models: kNDVI AND landscape structure and kNDVI OR landscape structure. With both the logical AND and OR models, each step has two processes: probability approach with the kNDVI and landscape structure. With the AND model, if the virtual population passes the probability of the kNDVI to move to a neighboring cell, then the probability of movement to a neighboring cell is based on the landscape structure. In the logical AND model, we used kNDVI first because that could reflect a seasonal change in the availability of food resources. With the OR model, if the virtual population passes the probability of the kNDVI, or passes that of the landscape structure, the virtual population can move to any neighboring cell.

Additionally, equivalence model simulation was conducted with all probability values (m1, m2, m3, and m4) set to 0.5.

Validation of the simulation results using the camera trap data

The results of the CA model simulation were validated by the presence/absence of the monthly records of sika deer collected with the cameras. The occurrence of a visit to a camera was determined using a generalized linear model with a binomial distribution (log link) and model selection based on Akaike’s information criterion (AIC). The explanatory variable was the theoretical number of simulated visits to a 1-km cell with a camera trap. If the AIC value of the model was > 2 points lower than that of the null model⁴⁵ (i.e., with no explanatory variable), the run was considered “correct”. The data from the kNDVI, landscape structure, AND/OR, and null/equivalence models were used. The number of “correct” runs of every 100 runs with each model was calculated. Therefore, all values could theoretically be 100.

Then, the predictive ability of the model was evaluated using the results considered as “correct” with the AIC. The AIC values of all runs were compared, where one simulation set used four variables. If the four models (i.e., kNDVI, landscape, AND, and OR models) were all “correct” in one run, the AIC values were compared and the lowest AIC value of the model was recorded. Notably, differences among the AIC values were not considered because the effectiveness of the model was already evaluated in the first validation procedure. Calculations for all months were conducted. Therefore, the maximum value among the four models was 100, assuming that the run was “correct” with the lowest AIC.

Finally, a map was generated of the theoretical number of visits by sika deer in each month based on the best performance among the four simulations. The map included the average number of theoretical visits over 100 runs. The results considered incorrect were not excluded because in real-world applications, simulated results are not evaluated.

All statistical analyses were performed using R software (ver. 4.0.2; https://www.r-project.org/, accessed on Feb. 10, 2023).

Source: Ecology - nature.com