The frequency and severity of extreme weather events such as heat waves are increasing globally, revealing ecological responses that provide valuable insights towards the conservation of species in a changing climate. In this study, we utilized data from two populations of GPS-collared female wood bison (Bison bison athabascae) in the boreal forest of northwestern Canada to investigate their movement behaviours in response to the 2021 Western North American Heat Wave. Using generalized additive mixed effect models and a model selection framework, we identified a behavioural temperature threshold for wood bison at 21°C. Above this threshold movement rates decreased from ~100 m/hr at 21°C to a low of ~25 m/hr at 39°C (150% decrease; -9%/°C). Extreme heat also contributed to changes in diurnal movement patterns, reducing wood bison movement rates and shifting the timing of peak activity from midday to early morning. These findings highlight the behavioural adaptations of female wood bison and underscore the need to understand the behavioural and physiological responses of cold-adapted mammals to extreme weather events. Subsequent effects of thermoregulatory behaviour may impact individual fitness and population viability, particularly at high latitudes where cold-adapted species are increasingly exposed to severe weather resulting from anthropogenic climate change.

Study system

We used movement data from females of two bison populations in northwestern Canada: the Aishihik population (n = 1,951; 95% CI = 1,688 - 2,295; Jung et al. 2023a) found in southwestern Yukon (~61.4°, -137.3°) and the Ronald Lake population (n = 272; T. Hegel, Alberta Environment and Parks, pers. comm. 2021) located in northeastern Alberta (~57.9°, -111.7°). The populations are approximately 1,500 km apart. The Aishihik population occupied an 8,000 km² area east of Kluane National Park and was primarily within the Traditional Territories of the Champagne and Aishihik First Nations and Little Salmon/Carmacks First Nation (Fig. 1; Jung et al. 2015a; Clark et al. 2016). The range of the Aishihik population had a cold and semi-arid climate, with snow cover extending from October to May (Jung 2020). Their range was mainly above treeline (~950 m above sea level [ASL]) and characterized by a mountainous landscape with several peaks ≥1,600 m ASL and alpine plateaus bisected by large river valleys and lakes (Jung 2020). Lowlands consisted of open-canopy boreal forest, wet sedge meadows, and relict boreal grasslands (Jung 2020). The Aishihik population was free ranging with unrestricted movements (Jung 2017), and experienced ecological and evolutionary processes such as competition (Jung et al. 2015a, 2015b, 2018) and predation (Jung 2011; Jung et al. 2023b). Anecdotal observations from field surveys indicate that calving for the Aishihik bison population typically occurs between early May and late June. Rarely, newborn calves have been observed as early as April and as late as December (Jung et al. 2019). The rut occurs in late summer (Jung 2020). Major summer diet components of the Aishihik bison population included sedges (Carex spp.), rushes (Juncus spp., and Eriophorum spp.), and grasses (Calamagrostis purpurea and Poa spp.) (Jung et al. 2015b). The population is subject to predation by wolves (Canis lupus) and brown bears (Ursus arctos); however, these events are rare (Jung et al. 2023b).

The Ronald Lake wood bison population was located in northeastern Alberta, immediately south of Wood Buffalo National Park, and occupied a range of ~4,500 km² within the Traditional Territory of Treaty 8 First Nations (Fig. 1; Tan et al. 2014; DeMars et al. 2020). The Ronald Lake population experienced short, warm summers with a mean daily temperature above 15°C and long, cold winters with a mean daily temperature below 10°C (Downing and Pettapiece 2006). Their home range was characterized by a mixture of upland deciduous, coniferous, and mixedwood forests and a network of lowland marshes, bogs, and other peatlands across undulating terrain (240 to 300 m ASL; Downing and Pettapiece 2006). Calving in this population typically occurred between 3 May and 28 June, and the rut occurred in late summer as is typical of wood bison (Komers et al. 1993; Hecker et al. 2020). Major summer diet components for this population included prickly rose (Rosa acicularis), fireweed (Chamerion angustifolium), currants (Ribes spp.) and willows (Salix spp.; Hecker et al. 2021b). Scat analyses indicated that bison constituted a relatively small portion of the summer diet of wolves (Dewart et al. 2020).

Temperature data

The 2021 Western North American Heat Wave occurred between late June and early July throughout western Canada (Overland 2021; Cotlier and Jimenez 2022). In northern Alberta, record highs were recorded between 29 June and 2 July 2021, including a high of 40.3°C in Fort McMurray (Environment and Climate Change Canada 2021a), at the southern edge of the range of the Ronald Lake population. In Whitehorse, southeast of the Aishihik population’s range, the heatwave peaked at 30.3°C on 28 June (Environment and Climate Change Canada 2021b). Given the short intensity of the heat wave experienced by both populations, we were interested in observing and comparing potential variation in bison movement rates before, during, and after the heat wave. Thus, we defined our study period as one month starting 15 June and ending 15 July 2021.

We acquired temperature data for the Aishihik population from Braeburn, Carmacks, Champagne, and Haines Junction weather stations (Fig 1; Government of Yukon 2023). We acquired temperature data for the Ronald Lake population from the Birch Mountain and the Mildred Lake weather stations (Fig 1; Alberta Agriculture and Forestry 2020). For all stations, we downloaded hourly average, maximum, and minimum temperatures (℃). Because fine-resolution temperature data were not available in our remote study regions, we averaged temperatures across stations for each population to capture potential regional variability.

Animal location data and daily movement rates

We used location data obtained from 34 adult female wood bison wearing GPS collars during summer 2021 to calculate movement metrics, including 21 and 13 adult females from the Aishihik and Ronald Lake populations, respectively. Adult female bison were chosen due to their importance in population dynamics and tendency to aggregate in large groups representative of a large portion of the population, particularly at this time of year (Gaillard et al. 2000). Moreover, we did not have data from male bison because deployment success rates for this age-sex class were dismal (Jung and Kuba 2015; Jung et al. 2018) and collaring operations were challenging (Jung and Larter 2021). However, we acknowledge that movement metrics for bison differ among sexes. Aishihik bison wore Lotek LiteTrack Iridium GPS collars (Lotek, Newmarket, Ontario, Canada), while those from the Ronald Lake population wore Vectronic Vertex Plus GPS collars (Vectronic Aerospace, Berlin, Germany). Bison were handled in accordance with approved protocols and procedures of the Alberta Wildlife Animal Care Committee and in accordance with the Yukon Wildlife Act. GPS collars for the Aishihik and Ronald Lake populations were programmed to provide a location at 60 min and 120 min fix intervals, respectively. Thus, we rarified locations from the Aishihik population to match the fix rate of the Ronald Lake population due to the influence of fix interval on the calculation of movement metrics (e.g., movement rate) for highly mobile animals (Prichard et al. 2014).

Our analysis included 5,193 locations from 21 GPS-collared bison from the Aishihik population and 1,407 locations from 13 GPS-collared bison from the Ronald Lake population. We filtered location data to align with the heat wave period, which conveniently corresponded to the time between calving and rut for bison, helping to minimize potential confounding effects of behavior (Melton et al. 1989; Komers et al. 1993; Jung et al. 2019; Hecker et al. 2020). Prior to analyses, we removed records with no coordinate information or low fix accuracy with a dilution of precision (DOP) value greater than 10 m (Bjørneraas et al. 2010), and manually reviewed and removed erroneous movements (i.e., movements far exceeding the maximum speed of bison). No bison were collared <2.5 months prior to our analyses so there was no need to censor data for capture effects on movement rates (Jung et al. 2019b). Fix rate success and location precision for GPS-collared bison in our region are typically >90% (Jung and Kuba 2015; Jung et al. 2018). For each bison, we calculated movement rate (v_i) in m/hr as:

v_i = l_i/t_i

where l_i is the step length between location i and location i + 1 and t_i is the time between location i and location i + 1 (Johnson et al. 2002; Sheppard et al. 2021).

 Modelling variation in movement rates across temporal scales

We created three sets of generalized additive mixed effect models (GAMMs) to assess non-linear associations between bison movement rates and temperature during the 2021 Western North American Heat Wave. We created a model set for each population to compare and assess heat wave responses between populations. Additionally, we made models combining both populations to generalize associations between bison movement and temperature. We used the R package “mgcv” (Wood 2022) in R version 1.3.1093 (R Core Team 2022) to assess the predictive power of temperature, time of day, and their interaction, related to movements.

We included a temperature covariate to investigate the biological response of bison to heat, and time of day and to quantify diurnal variation in movement rates. We fit cubic splines to avoid discontinuous movement rates throughout a 24-hour period (Wood 2006; Schmidt et al. 2016), and the interaction between temperature and time of day was fit as a tensor product interaction (time of day fit as a cubic spline) to investigate temperature-mediated diurnal variation in movement rates throughout the heat wave (Wood 2006). We considered Bison ID and Julian day as random effects in all models to account for variation among individuals and between days. In models including both populations, we also included population as a random effect to account for variation between populations. We fit all splines with low basis dimensions to avoid overfitting data (Sebastian-Azcona 2019; Sheppard et al. 2021). Concurvity is a common issue in generalized additive mixed models (GAMMs) that include time and time-varying covariates (e.g., temperature) (Wood 2006; Wood 2023). While there is no criterion for identifying an acceptable level of pairwise concurvity, a threshold of 0.5 or lower has been recommended for semiparametric models (Ramsay et al. 2003; Chien, 2009). Additionally, including or removing concurved variables can lead to changes in the sign or shape of smooth terms which may or may not be desired (Simpson 2013). To address this, we calculated and reported pairwise concurvity indices for each model and evaluated changes in the sign and shape of smooth terms to guide decisions on variable retention. Specifically, we retained all variables of ecological interest unless the inclusion of concurved variables (e.g., temperature or time of day) substantially influenced the sign or shape of smooth terms and/or reduced model performance.

We log-transformed movement rates (common logarithm) due to the expected high frequency of slow movements (Schmidt et al. 2016; Sheppard et al. 2021), specifying a Gaussian family identity link function. We fit models using maximum-likelihood estimation, which is necessary when comparing mixed-effect models (Pinheiro and Bates 2000). We visually confirmed assumptions and ranked candidate models using Akaike’s information criterion corrected for small sample sizes (AICc), selecting the model with the lowest AICc as the most parsimonious (Burnham and Anderson 2002). While GAMMs are powerful predictive tools, we wanted to make explicit statistical inferences about model predictions. Thus, we fit breakpoint regression models using the “segmented” package (Muggeo 2022) to approximate breakpoints and coefficients of individual splines from the most supported GAMM. We reported back-transformed effect sizes as the percent change in movement rates (m/hr) between predictor variable breakpoints. For interaction terms, we reported predicted movement rates at the 15^th and 85^th temperature percentiles to provide simplified ‘low’ and ‘high’ temperature scenarios.