Understanding how animals learn in novel environments is crucial for predicting behavioral responses to rapid environmental change, yet we lack knowledge about how long different behaviors take to develop and refine, and how exploration facilitates learning. We tested the exploration-refinement hypothesis using movement and diet data from GPS-collared bison (Bison bison; n=10) monitored for five years following reintroduction to Banff National Park. We examined how exploration influenced movement efficiency, habitat selection behaviors, and home range establishment. Different behaviors showed distinct learning trajectories. Movement efficiency in familiar areas reached an inflection point at 384 days, but when facing unfamiliar terrain, efficiency decreased as animals needed time to learn the new area. Habitat selection behaviors showed rapid initial improvement followed by extended refinement periods lasting up to three years. Exploratory movements occurred primarily in the first year (57%) but continued throughout the study and were positively correlated with improvements in habitat selection behaviors and home range stabilization, but not with movement efficiency. Our results demonstrate that exploration facilitates spatial learning—the process of acquiring and using information about environmental structure, distances, and relative positions to navigate effectively, locate resources, and avoid threats. This suggests that allowing individuals to explore while preventing large dispersal movements may be required for successful reintroductions. The extended timeline required for spatial behaviors to stabilize (3-4 years) indicates that behavioral refinement takes substantially longer than typically assumed, highlighting the need to align monitoring periods with the biological timeframes required for animals to learn novel environments.

Population size, animal handling, and collaring

The initial population of 16 wild-born bison (10 adult females, 6 adult males) from Elk Island National Park increased from 34 individuals in 2018 to 90 by 2022. Throughout the study period, 30 – 100% (from 100% initially to ~30% by study end) of the adult female population was monitored via GPS collars(Vectronic Aerospace, Berlin, Germany). We only included individuals with > 1 year of data. All data were cleaned to only include high quality 2-hour fixes (3 dimensional fixes with a Horizontal dilution of precision < 2). While bison initially moved as a single group, they exhibited typical fusion-fission dynamics as the population grew, regularly splitting into multiple subgroups throughout the study period (Figure S2, see Supplementary Methods for details on group dynamics). The animal study and capture protocol was approved by Parks Canada’s Animal Care Task Force under a 2016 Parks Canada Research and Collection Permit.

Metrics of spatial learning

Foraging behavior

We first evaluated the pace at which naïve individuals develop foraging behaviors. Bison are bulk feeders and mainly eat grasses (Steuter et al. 1995; Knapp et al. 1999) which predominantly grow in meadows within the study area. Thus, we considering foraging behaviors as each individual’s exposure to percent herbaceous cover and calculating the amount of time spent in meadows. For each individual’s GPS location, we extracted percent herbaceous cover from the Vegetation Resources Inventory (British Columbia Landcover Classification Scheme level 4; https://catalogue.data.gov.bc.ca/group/vegetation-resource-inventory [VRI level 4]). To quantify time spent in meadows, we first delineated potential foraging meadows for bison in this system. We delineated potential foraging meadows using a combination of landcover type, percent herbaceous cover, and shrub crown closure data. Meadows were required to meet minimum size and vegetation quality thresholds to be considered viable foraging habitat (full delineation criteria in Supplementary Methods). We then intersected GPS locations with the meadow polygons and summed the number of locations falling within a meadow over the course of each month an individual was monitored.

Movement behavior

To evaluate the pace at which individuals develop movement behavior, we evaluated how closely individuals followed to predicted least-cost paths when moving between heavily used foraging/resting areas. First, we delineated the population-level core use areas (i.e., 50% utilization distribution polygons) on the landscape by fitting a kernel density estimator (adehabitatHR package) to all bison GPS locations (i.e., the entire study period July 2018 – September 2023) using a 30 m grid and delineating the cells that included the top 50% of volume of use (i.e., core area; Calenge 2006).

We analyzed movement efficiency by comparing observed travel paths between core areas to predicted least-cost paths based on topographic slope. Least-cost paths were calculated using elevation-based cost functions optimized for mountainous terrain (technical details in Supplementary Materials).

We compared observed traveling trajectories with predicted least-cost paths between core areas over time. To render GPS location data to such a framework, we first identified whether each bison location was inside or outside the core area and removed locations within core areas. We then identified consecutive steps between core areas in the movement data and connected those sequential locations to form a trajectory. We compared observed paths with least-cost paths using the normalized path deviation index. The index compares how far apart the observed movement path is from the predicted least-cost path along the entirety of the trajectory (the area between the two paths divided by the distance of the shortest path to normalize the metric so that indexes can be compared across individuals with different path lengths). The normalized path deviation index is reported as the percent deviation between the two paths (i.e., if normalized path deviation index is 30%, the average distance between the two paths is 30% the length of the distance of the shortest path; Lewis 2023). 

Emergent patterns: Home range settlement and diet quality

We evaluated the pace at which emergent patterns such as home range area and overlap and diet stabilize over time. First, we estimated individual monthly home ranges using an auto-correlated kernel density estimator (Fleming and Calabrese 2017, Signer et al. 2019). We removed any individual-month where the individual in question had <20 days of GPS locations. As a result of a few missing months, we only compared area and overlap for consecutive months. We compared monthly home range overlap using Bhattacharyya's affinity index which incorporates both the spatial extent and the intensity of space use, providing a comprehensive measure of the similarity between the utilization distributions of individuals (Fieberg and Kochanny 2005, Signer et al. 2019). Second, to estimate diet quality, we analyzed a number of diet quality metrics from fecal samples collected opportunistically over time (from December 2018 to November 2020) and across different areas of the study area. Fecal samples were analyzed for diet quality by the Grazing Animal Nutrition (GAN) Lab in Texas (USA). Samples were tested for fecal nitrogen, digestible organic matter, phosphorus, and nitrogen, all metrics that have been shown to be relevant to bison diet and foraging (Geremia et al. 2019).

Modelling approach

To evaluate the pace of learning, we plotted learning curves, which have often been used in behavioral ecology experiments to test for learning (e.g., Papachristos and Gallistel 2006, Dukas 2008) and more recently to quantify the development of movement behavior (Hertel et al. 2023). Similar to the methods of Hertel et al. (2023), we evaluated the pace of learning and whether foraging behavior or movement behaviors are learned more quickly, by comparing a set of 3 models for each foraging, movement, and home range variable. We evaluated how time influenced: 1) exposure to herbaceous cover, 2) time spent in meadows, 3) diet quality, 4) deviation from least-cost paths, and 5) home range overlap. If the distribution of our explanatory variable was right-skewed, we log transformed it prior to fitting the model. We fit each variable as a separate model where time was a) a linear predictor, b) a squared predictor, and c) a cubed, polynomial predictor and compared models using Akaike information criterion (AIC; Aho et al. 2014), allowing the general pace and shape of each learning curve to be determined. To account for repeated measures of individuals and seasonal effects, each model except for diet quality was fit with a random intercept term of animal ID and season. Season was operationally defined as spring (March – May), summer (June – August), fall (September – November) and winter (December – February). Diet quality samples were not tied to individual animals, therefore the random intercept term was only season.

We then determined changepoints in the predicted learning curves (i.e. inflection points) to identify when the pace in learning significantly changed (Muggeo et al. 2014). We fit segmented regression models to the predicted relationship of each variable (e.g., time in meadows, adherence to least-cost paths) to time based on the top learning curve model from above. By identifying changepoints, we can detect significant behavioral changes, in this case when learning peaked or began to slow down. Segmented models are regression models where the relationship between the response and predictor variables are piecewise linear (i.e., two or more straight lines connected at unknown values, the breakpoints or changepoints; Muggeo et al. 2014, Muggeo 2017). We set the number of change points to reflect the predicted curve from the best fitting model for each variable (e.g., for an asymptotic curve, we used one change point, but for a sigmoid or similarly complicated curve, we used two change points, etc.).

To evaluate whether exposure to novel places (i.e., exploration) leads to spatial learning, we fit models with each of our foraging and movement behavior variables as a function of the number of exploratory events undertaken by the bison throughout the study period. We fit each variable as a separate model where the number of exploratory events was fit with the same structure as the learning curve (i.e., if the top learning curve model was a polynomial, we fit exploratory events as a polynomial term). Similar to the learning curve models, each model included a random intercept term for animal ID and season.