We modeled the distribution of northern pika in central Hokkaido, Japan, because detailed habitat thermal measurements have been conducted in this area (Sakiyama et al. 2021; Sakiyama and García Molinos 2024). The northern pika in Hokkaido is located at the southern margin of its distribution and therefore its population is likely among the most vulnerable of the species to climate warming. Although northerly neighboring populations are found across the Sea of Japan in Sakhalin and Primorsky, the Hokkaido population is essentially a disjunct population. As such, the accessibility and availability of climate change refugia within Hokkaido is of significant importance for the Hokkaido population, while it is likely to respond to climate change by shifting their distribution upslope or by utilizing cool microclimates. However, previous studies have not predicted its future distributions.

To model the distribution of northern pikas, we compiled occurrence data in central Hokkaido from previous literature. Given almost all the compiled data comprised presence (but not absence) locations, we decided to apply a presence-only algorithm for distribution modeling (see Model development section). We only used presence locations reported on high resolution maps (i.e. 1/25,000 and 1/50,000 scale) for analysis since we were interested in modeling at fine scales and those surveyed between 1961–2010 to match the timeframe with the environmental data. Overall, a total of 305 presence points encompassing a wide elevational gradient (50–2,210 m) were obtained from the literature (Suppl. material 1: table S1, fig. S1). Survey years spanned between 1963 and 2009 (Suppl. material 1: fig. S1a, b). Detection signs of northern pika presence were clearly mentioned for almost all points (298 points, 97.7%), with the vocalization representing the main method enabling detection of 241 points (79.0%) while other methods were also used for detection: direct sightings (78 points, 25.6%), haypiles (47 points, 15.4%), feces (35 points, 11.5%), and eatmarks (9 points, 3.0%).

As calibration extent for this study, we considered the area corresponding to the known geographical range of the species comprising the continuous mountainous areas in central Hokkaido delimited by the flatlands (i.e., river valleys) as the extent borders. Specifically, we considered the Abira and Ishikari River systems as the western, the Teshio River system as the northwestern, the Okoppe River system and the Sea of Okhotsk as the northeastern, the Tokoro and the Tokachi River systems as the eastern, and the Pacific Ocean as the southern borders of the extent (Suppl. material 1: fig. S1a). Next, we removed irrelevant land cover types from within the study extent area such as urban land, agricultural fields, and aquatic area (e.g., lakes) since northern pikas require pristine montane conditions for inhabitation. The resulting filtered study extent area covered approximately 21,000 km² in central Hokkaido (Suppl. material 1: fig. S1a).

To understand how differences in spatial resolutions of climate data affect model performances and predictions of species distribution, we constructed coarse-resolution and fine-resolution models (Fig. 1; Suppl. material 1: table S2). In the coarse-resolution model, we used environmental predictors at a resolution of 30 arcsecs, which converts to a grid cell of 691 m by 926 m (longitude × latitude) and size of 639,682 m². We consider this setting to be a coarse resolution for the northern pika because their home range is much smaller, typically within the range 1,822–11,530 m² (n = 4; Onoyama et al. 1991). We also note that the mountainous terrains of our study area are significantly smoothed at this resolution, with median and maximum elevational differences within one 30-arcsec grid cell as large as 182 m and 770 m, respectively, when computed based on 10-times finer (3 arcsecs) elevation data (Suppl. material 1: fig. S2). Since species occurrence-environment relationships might be confounded by such a large degree of simplification, we developed the fine-resolution model using environmental data at 3 arcsecs, which converts to a grid cell of 69.1 m by 92.6 m and size of 6,397 m². We consider this setting to reflect the local ambient air conditions experienced by pika individuals. Based on the fine-resolution model, we then generated two habitat suitability predictions to account for the duality in thermal conditions experienced by the northern pika. The first corresponds to the raw predictions based on the fine-resolution model, which represents the above-ground habitat suitability considering the effect of local ambient air. The second prediction incorporates information on rock-interstice microclimates that are buffered from the ambient air, to further predict sub-ground habitat suitability.

Figure 1. 

Schematic diagram describing the different model settings of the study. We considered two spatial resolutions for analyzing the northern pika distribution: coarse-resolution (30 arcsecs), representing conventional, temperature data; and fine-resolution (3 arcsecs), which is more adequate for the species due to their small home range size. We used widely available temperature data for developing the coarse-resolution model, while downscaled temperature data, which we considered to reflect the ambient air temperature, was used for the fine-resolution model. We further used the fine-resolution model to evaluate the habitat suitability in areas with rocky landforms (cells surrounded by thick white lines in the microclimate layer) by incorporating the rock-interstice microclimate, which is known to be buffered from the ambient air temperature. Finally, we conducted hindcasting and future projection analyses to predict distribution changes over time.

We used the terrain slope angle as a topographical predictor in the SDMs because it has been shown to have a positive relationship with the northern pika distribution (Sakiyama et al. 2021). We computed the slope angle raster in QGIS 3.34 (QGIS Development Team 2023) using the MERIT DEM (Yamazaki et al. 2017), a digital elevation model provided at a spatial resolution of 3 arcsecs. Slope angles for the coarse-resolution setting were calculated using a DEM aggregated to 30 arcsecs (Suppl. material 1: table S2). We also considered surface geology as a categorical binary predictor because 99% of the known distribution of the northern pika lies on igneous and metamorphic rocks and moraines (Kawabe 1992). To do so, we assigned whether any of these geological types occurred per grid cell within the study area based on the Seamless Digital Geological Map of Japan (1:200,000) V2 (Geological Survey of Japan, National Institute of Advanced Industrial Science and Technology 2017) (Suppl. material 1: table S2).

To model and predict the distribution of northern pikas, we computed the mean daily temperature and mean diurnal range over July and August to account for the effects of thermal magnitude and variation during the summer. We focused on summer conditions because the effects of climate change on cold-adapted species are likely to be particularly strong during this season. We used the average following previous studies exploring the effects of thermal conditions on the northern pika distribution (Sakiyama et al. 2021; Sakiyama and García Molinos 2024). Briefly (see Suppl. material 1: appendix S1 for a detailed explanation), we first obtained temperature data at 30 arcsecs using the chelsa-cmip6 R package, which provides access to downscaled Global circulation model (GCM) projections (Karger et al. 2023). Four GCMs that are known to reproduce well the climatic conditions of our study area were selected (i.e. ACCESS-CM2, IPSL-CM6A-LR, MIROC6, and MRI-ESM2-0) (Shiogama et al. 2021). We used temperature data for 1981–2010 to model the baseline distribution. To create the bioclimatic variables, we adapted the computation formulas used in ANUCLIM, a software package used to build spatial climatic data (Xu and Hutchinson 2011, 2013), to the summer season, and calculated mean summer temperature and mean summer diurnal range. The resulting variables were then used in the coarse-resolution model (Suppl. material 1: table S2).

For the fine-resolution model, we downscaled the temperature data obtained at 30 arcsecs to 3 arcsecs through a kriging analysis using the krigR R package (Davy and Kusch 2021; Kusch and Davy 2022). This approach first models the relationship between temperature and elevation at a coarse resolution (30 arcsecs), which is then used to interpolate the temperature to the target resolution (3 arcsecs) using elevation data with the same resolution. We used the raw MERIT DEM (Yamazaki et al. 2017) dataset, provided at 3 arcsecs, as the elevation data for the interpolation and the aggregated version at 30 arcsecs to model the coarse-resolution relationship (Suppl. material 1: table S2).

Finally, we considered that northern pikas can utilize the rock-interstice microclimate in the sub-ground space, which are known to be buffered from ambient air temperature (Millar et al. 2016; Smith et al. 2016; Sakiyama and García Molinos 2024), in areas with rocky landforms. However, given the absence of sufficient spatial distribution information regarding rocky landforms throughout our study extent, we selected three representative regions where we had sufficient data, separated by elevational ranges: Biei (1,100–2,000 m), Shikaribetsu (600–1,250 m), and Oketo (400–800 m) (Fig. 3d). For Biei, we were able to map the rocky landforms directly using aerial images. For Shikaribetsu and Oketo, we referred to the mapped occurrence of rocky landforms in previous literature (Kojima 2004; Kurumada 2010). For these regions, we then used an empirical relationship between ambient and interstice temperatures previously derived from northern pika taluses in our study area (Sakiyama and García Molinos 2024) to convert downscaled temperature (i.e., ambient air) into rock-interstice temperature for each of the 3 arc second grid cells containing rocky landforms (Suppl. material 1: fig. S3, appendix S1). Since rocky landforms are patchily distributed within each region, our assumption was that northern pika in grid cells without rocky landforms experience only ambient air temperatures, while those in grid cells with rocky landforms experience dual thermal conditions of ambient air and rock interstice temperatures. We validated our downscaling and microclimate incorporation approach by comparing thermal conditions measured in northern pika habitats with values of each bioclimatic variable in the grid cells where those observations had been made (see Suppl. material 1: appendix S2). Results confirmed that the fine-resolution data were moderately accurate in reflecting observed ambient temperature and that the microclimate data were highly accurate in reflecting observed rock interstice temperatures (Suppl. material 1: fig. S4; Fig. 5).

We performed the distribution modeling analysis of northern pikas in Hokkaido using maximum entropy (MaxEnt; Phillips et al. 2006). While model averaging is often used in SDM studies to remove potential bias in the predictions emerging from the algorithm used, we decided to use a single algorithm considering the overall methodological complexity and computation time. We chose MaxEnt because it has been shown to consistently rank among the best performing algorithms in comparative studies (Valavi et al. 2022). Models were developed separately for coarse- and fine-resolution settings considering years 1981–2010 as the baseline period using presence points obtained during this period for model fitting. Prior to model fitting, we reduced the effect of spatial sampling bias by thinning the presence data using the spThin R package (Aiello-Lammens et al. 2015) at a distance of 500 m as presence locations could be viewed as independent at this distance considering the complex mountainous terrain and the limited mobility of the northern pika. Considering that the diagonal distance of coarse grid cells is larger than this thinning distance and duplicate sampling thus could occur, we randomly selected a single point for each cell after thinning. This pre-selection process resulted in a total of 91 presence points remaining. We generated 10,000 background points throughout the study area, which are used in MaxEnt for characterizing the environmental conditions within the study area. To avoid multicollinearity issues, we computed the pairwise Pearson’s correlation coefficients among the continuous predictors for both the coarse- and fine-resolution settings and confirmed that all pairs did not have a strong correlation (|r| > 0.7) (Dormann et al. 2013) (Suppl. material 1: fig. S6).

Since MaxEnt models with default settings can lead to overfitting (Anderson and Gonzalez 2011), we initially tuned the model hyperparameters to control the complexity while fitting well to the data. To do this, we used various hyperparameter combinations of feature classes, which determines the shape of the response curves, and regularization multipliers, which controls the model complexity (Phillips et al. 2006; Phillips and Dudík 2008), and selected the combination that resulted in the highest predicative accuracy. For the feature class, we considered linear (L), quadratic (Q), product (P), and/or hinge (H) characteristics, and for the regularization multipliers, we used values ranging from 0.2 to 2.0 with intervals of 0.2. The predictive accuracy was assessed for each model by separating the presence data into training (80%) and test (20%) data and computing the area under the receiver operating characteristic curve (AUC), which indicates the accuracy condition on a range of 0 to 1 (Fielding and Bell 1997). All tuning processes were conducted using the optimizeModel function in the SDMtune R package (Vignali et al. 2020). After hyperparameter tuning, we evaluated the optimized model’s predictive accuracy based on 5-fold cross validation using Continuous Boyce Index (CBI), which indicates whether the predicted habitat suitability is consistent with the distribution of presence points (Hirzel et al. 2006). CBI values range from -1 to 1, with positive values indicating an accurate prediction, values close to zero indicating inaccurate prediction (i.e., prediction is close to random), and negative values indicating counter predictions.

For both coarse- and fine-resolution models, we predicted the habitat suitability for the baseline period across the study area and compared the predictions by computing the difference (Coarse – Fine). We also generated binarized maps of predicted species occurrence (i.e. presence or absence) based on the habitat suitability prediction across the whole study area. Since discriminant thresholds were likely to differ between the optimized models for coarse- and fine-resolutions and could be problematic for comparison, we instead used a sensitivity analysis approach by applying successive thresholds from 0 to 1 at 0.01 increments to discriminate presence or absence. Then, binary predictions were compared between coarse- and fine-resolution models by assessing the proportion of predicted species’ range in the study area calculated as the total amount of cells predicted as present divided by the number of all cells. We used proportion instead of simply calculating the total area of presence because the total study area was subtly different between the scales due to processing procedures of environmental data.

We further predicted the habitat suitability for the sub-ground space for the three selected regions by replacing local ambient air temperatures with rock interstice microclimate using the fine-resolution model. The resultant prediction based on rock-interstice microclimate was mapped and compared with those based on the coarse-resolution model and fine-resolution model using ambient air temperature by considering the average habitat suitability at the region-wide scale and in areas with rocky landforms. In each region and prediction setting, the region-wide habitat suitability was assessed by sampling the predicted values using the cell centroids of the coarse-resolution data, while the rocky landform habitat suitability was evaluated by sampling the predicted values at random points generated within areas with rocky landforms. The habitat suitability values were then compared among the prediction settings by conducting a beta regression analysis using the betareg R package (Cribari-Neto and Zeileis 2010) followed by a pairwise comparison between prediction settings using the Tukey method for p-value adjustments.

We use a hindcasting approach to evaluate the transferability of the models to project distributions over time (see Nogués‐Bravo 2009) by projecting the distribution of the species to the past and comparing it to known past presence points from the same past period for validation (Hodel et al. 2022). Although this approach is still rarely used in SDM studies (Yates et al. 2018), it allowed us to assess how confident we can be using the models to predict distributions for a time period different to that the model was calibrated for. We projected the distribution of northern pikas for the years 1961–1980, a period for which early presence records are available (Suppl. material 1: table S1, fig. S1). The same procedures as the baseline climate data were used for processing past data. We evaluated the model’s predicative accuracy using the CBI (Hirzel et al. 2006).

To project the future distribution of the northern pika, we obtained climate data for two future time periods, 2041–2070 and 2071–2100, considering two emission scenarios, SSP1-2.6 and SSP5-8.5, representing a globally sustainable and a fossil-fuel based development scenario, respectively. We predicted the future distribution for the whole study extent based on the coarse-resolution model and the fine-resolution model using ambient air temperature by applying future bioclimatic data, which were processed based on the same procedures for the baseline climate data (i.e., we assume that future topographical and geological variables will remain the same as present). We then binarized the suitability predictions to presence or absence by applying the same approach used for the baseline period to calculate the proportion of predicted species’ range.

Moreover, we also projected the future distribution of the northern pika for the three selected regions to assess how rock-interstice microclimate contributes to retaining suitable conditions in the future. When processing the future microclimate data, we assumed that the buffering effect of the rock interstices (i.e., the empirical relationship between local ambient and rock-interstice temperatures) will remain invariable over time. We believe that this assumption is robust because the field thermal measurements from which this relationship was derived contained a large elevational gradient of 350–2200 m and the model residual was stable across the ambient temperature gradient (Sakiyama and García Molinos 2024), although we note that some extrapolation did occur for some lower elevation areas where future conditions exceeded the observed thermal gradient. We binarized the suitability predictions for both fine-resolution and microclimate using the maximum true skill statistic computed for the fine-resolution model (Santini et al. 2021). The consensus among the resulting binarized predictions from each of the four GCMs was then computed by considering a cell as present whenever that cell was predicted as present by two or more individual GCMs (i.e., ≥50% consensus). The resultant prediction was compared between the fine-resolution prediction using ambient air temperature for each period and emission scenario.

The optimized models for both the coarse- and fine-resolution yielded CBI values of 0.73 and 0.80, respectively (Table 1). Among the multiple variables used, mean summer temperature was the most important predictor explaining the northern pika distribution in both models with a permutation importance of 74.5% and 76.2% in the coarse- and fine-resolution models, respectively (Fig. 2). Mean summer temperature showed a negative effect on occurrence probability in both models, while the slope of the response curve was slightly steeper in the fine-resolution model (Fig. 2). In both models, surface geology was the second most important predictor with 20.5% and 17.8% of permutation importance, respectively, which exhibited a positive relationship with occurrence probability (Fig. 2). The contribution of daily thermal range and slope was less than 5% in both models (Fig. 2).

Figure 2. 

Response plots of the model predictors for the best-fit coarse-resolution (top) and fine-resolution (bottom) models. Permutation importance of the model predictors is provided below each variable name on the horizontal axis.

Mountainous areas were predicted as suitable areas for the northern pika, with higher elevations predicted as highly suitable areas while low elevations were predicted as areas with low to moderate suitability in both the coarse- and fine-resolution models (Fig. 3a, b). Comparing the predictions, the habitat suitability was generally higher in the coarse-resolution model than the fine-resolution model (Fig. 3c). Difference in habitat suitability was relatively large in mountainous areas, particularly in the slopes rather than the summits (Fig. 3c). The proportion of predicted species’ range for both models exhibited similar inverse curves along the binarization threshold. However, the proportion was always higher in the coarse-resolution model than the fine-resolution model across the whole range of binarization threshold (Fig. 4a).

In the three selected regions, the regional habitat suitability based on both the coarse-resolution model and fine-resolution model using ambient air temperature was highest in Biei, followed by Shikaribetsu and Oketo, reflecting the regional elevational difference (Fig. 5). Within each region, the habitat suitability reflected the local elevational variation, with higher elevation areas predicted as more suitable (Fig. 5). In the fine-resolution model prediction, accounting for rock-interstice microclimate resulted in boosting the habitat suitability from that based on ambient air temperature in areas with rocky landforms in all regions (Fig. 5) Consequently, the habitat suitability in areas with rocky landforms was always significantly higher in the fine-resolution prediction based on rock-interstice microclimate than the coarse-resolution prediction and the fine-resolution prediction based on ambient air temperature (Fig. 5). However, the region-wide habitat suitability was not significantly different among the three prediction settings in two regions (Biei and Oketo), and it was only found to be different between the coarse-resolution prediction and the fine-resolution prediction based on ambient air temperature in Shikaribetsu (Fig. 5).

Both coarse- and fine-resolution models were supported by high transferability to predict distributional changes over time, as indicated by the positive CBI values in the hindcasting analysis (Table 1). This enabled projecting future distributions with confidence, in which both models indicated inverse curves in the proportion of predicted species’ range (Fig. 4b). Compared to the baseline period, future range was predicted to decrease by more than 60% in both time periods and emission scenarios considered (Suppl. material 1: fig. S7). Between the two future time periods, the difference in the habitat decrease rate was relatively small under the SSP1-2.6 scenario, suggesting that the northern pika’s range will shrink until 2041–2070 but will then remain relatively stable until 2071–2100. However, the habitat decrease was apparent under the SSP5-8.5 scenario, suggesting that suitable habitats will continue to decrease until 2071–2100. (Fig. 4b; Suppl. material 1: fig. S7).

Patterns of distributional change based on the fine-resolution model largely differed among the three selected regions between predictions based on ambient air temperature and rock-interstice microclimate. Based on ambient air temperature, the northern pika presence was predicted in Biei in both future time periods and scenarios and in Shikaribetsu under the SSP1-2.6 scenario, while the whole region was predicted to be absent in Oketo under all future settings (Fig. 6). Conversely, the northern pika presence was predicted in almost all regions in all future settings when accounting for rock-interstice microclimate, except in Oketo for the far future under the SSP5-8.5 scenario (Fig. 6). Discrepancies between predictions based on ambient air temperature and rock-interstice microclimate were high towards the lower elevations (Shikaribetsu and Oketo) and under the higher carbon emission scenario (SSP5-8.5).

Figure 3. 

Predicted habitat suitability of the whole distribution of the northern pika in Hokkaido, Japan, for the baseline period (1981–2010), based on (a) coarse- and (b) fine-resolution models. In (c), the habitat suitability difference between these two predictions (coarse – fine) is provided. Black points represent presence locations used for the analysis. In (d), the elevation of central Hokkaido is shown in the background with locations of three regions used for the subsequent analysis depicted in squares. The black border line denotes the model’s calibration extent with light blue areas in (d) representing areas excluded from the analysis (see Methods for details).

Figure 4. 

Proportion of predicted species’ range along the successive binarization threshold applied at 0.01 increments (see Methods for details) for the (a) baseline and (b) future periods. The lines and shade colors represent model settings (orange = coarse-resolution, green = fine-resolution). The future panels are divided by time period (2041–2070 and 2071–2100) and emission scenario (SSP1-2.6 and SSP5-8.5), while the filled envelopes represent the range of proportions based on four Global circulation models.

Figure 5. 

Predicted habitat suitability of the northern pika in (starting from top row) the Biei, Shikaribetsu, and Oketo regions for the baseline period (1981–2010), based on (starting from left column) the coarse-resolution model, fine-resolution model using ambient air (downscaled temperature), and fine-resolution model using rock-interstice microclimate. Areas with rocky landforms are depicted by black lines in the microclimate panels. The elevation for each region is shown in the fourth column. The last column shows boxplots comparing the predicted habitat suitability among the settings at the regionwide-scale (left) and within rocky landforms (right). The abbreviations on the x-axis correspond to: C = coarse-resolution model; F-aa = fine-resolution model using ambient air; and F-m = fine-resolution model using rock-interstice microclimate.

Figure 6. 

Predicted future distribution of the northern pika considering microclimates in (starting from top row) the Biei, Shikaribetsu, and Oketo regions based on the consensus approach. The 1^st and 2^nd columns are predictions for the 2041–2070 and 2071–2100 under the SSP1-2.6 scenario, and the 3^rd and 4^th those under the SSP5-8.5 scenario. Dark gray areas represent where suitable habitat is lost in the future, while light gray areas represent unsuitable areas at the baseline period. Areas with rocky landforms are depicted by black lines. The light blue areas represent suitable habitats supported by ambient air outside rocky landforms. Within rocky landforms, the dark blue areas represent suitable habitats supported by both ambient air and rock-interstice microclimate, whereas the pink areas represent those supported only by microclimate.