Moose at their bioclimatic edge alter their behavior based on weather, landscape, and predators.

Populations inhabiting the bioclimatic edges of a species' geographic range face an increasing amount of stress from alterations to their environment associated with climate change. Moose Alces alces are large-bodied ungulates that are sensitive to heat stress and have exhibited population declines and range contractions along their southern geographic extent. Using a hidden Markov model to analyze movement and accelerometer data, we assigned behaviors (rest, forage, or travel) to all locations of global positioning system-collared moose (n = 13, moose-years = 19) living near the southern edge of the species' range in and around Voyageurs National Park, MN, USA. We assessed how moose behavior changed relative to weather, landscape, and the presence of predators. Moose significantly reduced travel and increased resting behaviors at ambient temperatures as low as 15 °C and 24 °C during the spring and summer, respectively. In general, moose behavior changed seasonally in association with distance to lakes and ponds. Moose used wetlands for travel throughout the year, rested in conifer forests, and foraged in shrublands. The influence of wolves Canis lupus varied among individual moose and season, but the largest influence was a reduction in travel during spring when near a wolf home range core, primarily by pregnant females. Our analysis goes beyond habitat selection to capture how moose alter their activities based on their environment. Our findings, along with climate change forecasts, suggest that moose in this area will be required to further alter their activity patterns and space use in order to find sufficient forage and avoid heat stress.

Climate change is forecast to drive dramatic changes to the distribution (Parmesan 2006; Chen et al. 2011) and abundance of species worldwide (Dirzo et al. 2014; Ripple et al. 2016). Species will face temperatures that more frequently exceed their thermal thresholds (Pörtner 2001), and vegetative communities will shift in potentially rapid and unpredictable ways, thus altering ecological communities (Parmesan and Yohe 2003). Ecosystem shifts can lead to greater abundance of competitors and predators (Huey et al. 2009; Gilman et al. 2010), or changes in predator behavior that modify competition or predator–prey dynamics (Post et al. 1999). Animals inhabiting areas at the bioclimatic edges of their range will encounter the most dramatic and earliest effects of climate change. Understanding how individual animals respond to existing variation in those environmental conditions predicted to sustain the greatest changes will provide the best indication of the future persistence (Franco et al. 2006), abundance trends (Forchhammer et al. 2001), and range shifts (Hampe and Petit 2005; Hickling et al. 2006) of a population.

For mammals, it may take years for the effects of climate change to result in detectable population declines (Forchhammer et al. 2001; Parmesan 2006; Mason et al. 2014). Studies that attempt to understand changes in mammal behavior using traditional approaches that lack direct observation, such as habitat selection studies, may not fully capture the altered behaviors of individuals coping with novel climatic and ecological conditions. For example, animals may select for multiple habitats, but this simple association may not capture how an animal uses each (e.g., one habitat may provide thermal refuge and the other foraging opportunities; Street et al. 2016). Likewise, these models may not identify whether an animal is only using a habitat at certain times to reduce its risk of predation (Latombe et al. 2014). Fortunately, improved methodologies of analyzing animal movement can help us discern what types of behaviors occur in different habitats, and how these behaviors may change based on current ecological conditions (Edelhoff et al. 2016; Gurarie et al. 2016).

Hidden Markov models (HMMs) and similar approaches (e.g., state–space models) can classify animal behavior based on movement characteristics derived from the locations of individual animals (Morales et al. 2004; Forester et al. 2007; Patterson et al. 2008; McClintock et al. 2012a; Beyer et al. 2013). These models assign an unobservable behavioral state, assumed to be the true behavior of the individual, to each location of the animal. Aside from movement characteristics, HMMs can incorporate ancillary data, such as physiology measurements from biologgers and activity counts from accelerometers (McClintock et al. 2012b; Fehlmann et al. 2017; Leos-Barajas et al. 2017) to more accurately assign behavioral states or define additional states (Nams 2014). Behaviors linked with each relocation can be associated with outside biotic and abiotic factors to gain inference on how animals respond or activity levels change when encountering different conditions and habitats (e.g., Russell et al. 2015).

Moose Alces alces are keystone herbivores and are a particularly good candidate species to examine how a changing environment may affect the behavior of a mammal at the edge of its bioclimatic threshold. Moose are physiologically sensitive to heat (Renecker and Hudson 1986, 1989; McCann et al. 2013) and dramatically alter their habitat selection when ambient temperatures rise (Schwab and Pitt 1991; van Beest et al. 2012; Street et al. 2015, 2016) by selecting for habitats that act as thermal refuges (Dussault et al. 2004; McCann et al. 2016). Energetic costs are tied to higher ambient temperatures primarily in two ways: 1) direct costs of dissipating heat through increased metabolic rates, including in extreme cases, panting to reduce excess heat (Renecker and Hudson 1986; Renecker and Schwartz 1998), and 2) indirect costs of forgoing foraging opportunities while resting to avoid overheating (Street et al. 2016). These costs may lead to a reduction in physical condition and thus an increased risk of mortality (Renecker and Hudson 1992; Joly and Messier 2004; van Beest and Milner 2013).

Throughout the southern extent of their range, moose show signs of poorer health (Ruprecht et al. 2016), reduced fecundity (Monteith et al. 2015; Ruprecht et al. 2016), lower calf survival (Grøtan et al. 2009; Severud et al. 2015), or range contractions (Dou et al. 2013). Summer heat stress was found to be one of several potential factors in nearly extirpating a moose population in northwestern Minnesota, USA (Murray et al. 2006). Similarly, moose populations in the northeastern part of Minnesota have experienced population declines of more than 50% between 2005 and 2016 (DelGiudice 2016). Factors such as warmer year-round temperatures increased parasite loads and diseases from the northward advance of white-tailed deer Odocoileus virginianus, and an increasing wolf Canis lupus population are thought to all play an increasing role in declining moose populations (Murray et al. 2006; Lenarz et al. 2009, 2010, DelGiudice et al. 2011; Mech and Fieberg 2014). In addition to the current threats facing moose, the boreal forests of northern Minnesota that moose inhabit are predicted to shift dramatically to drier and more open cover types in the coming decades (Galatowitsch et al. 2009), thus reducing the ability of moose to find habitats typically used for foraging (Street et al. 2015) and thermal refuges (McCann et al. 2016).

Here, we utilize HMMs to examine how moose living in a protected area along the southern extent of the species’ range alter their behaviors when responding to variability in weather (including ambient temperature), and the presence of wolves, along with encountering multiple habitats and landscape features. Our HMM characterized all moose global positioning system (GPS) locations as one of three behavioral states: traveling, resting, or foraging. We incorporated activity levels from accelerometers in the GPS collars to better differentiate foraging versus resting behaviors that may be indistinguishable based on movement data alone. Our approach enabled us to ask three important questions about moose behavior across seasons: 1) can we capture the thermal thresholds at which moose alter their behavior, 2) does moose behavior change in areas regularly used by wolves, and 3) how do moose change their behavior in different primary habitats and relative to landscape features such as lakes? Many studies have attempted to determine the thermal thresholds of moose, but they were often conducted in laboratory or captive settings. Here we develop a fuller understanding of moose behavior by analyzing behavioral changes at a fine scale as moose encounter different stressors and habitats in the wild. We expand on our findings in the context of long-term climate forecasts of the region and discuss how the southern extent of the geographic distribution of moose is expected to change in the future.

Materials and Methods

Study area

All GPS-collared moose inhabited the Kabetogama Peninsula of Voyageurs National Park (VNP; 48°30′N, 92°50′W; n = 11; moose-years = 17) and the adjacent Rat Root Lake area (n = 2; moose-years = 2). The combined study area is located in north-central Minnesota, USA along the southern edge of the geographic distribution of moose. The Kabetogama Peninsula (305 km2) is a roadless protected area, but maintains hiking trails in the summer and snowmobile trails during the winter months. The Rat Root Lake area (∼12 km2), which is located near the western edge of VNP, is primarily state forest land with a low road density. Area summers, while relatively short in duration, are humid and reach average daily July temperatures of 18.8 °C (10-year average: 2003−2012). Winters are typically dry and cold with an average daily temperature of −13.7 °C (10-year average: 2003−2012) during January. The 10-year averages are based on the National Oceanic Atmosphere Administration’s National Climatic Data Center [cited August, 10th 2016 (https://www.ncdc.noaa.gov/)] Global Climate Station Summary for International Falls, MN, USA. Lakes and ponds throughout the park are generally ice-covered from late November until late April or early May (Kallemeyn et al. 2003).

The landscape of the VNP ecosystem contains numerous lakes and ponds (23.3% open water, based on % areal coverage) with islands containing rocky outcrops and shoreline bluffs and a mosaic of beaver-influenced wetlands (Johnston and Windels 2015a). Vegetation is primarily forest (36.4%), dominated by aspen (Populus spp.), white birch Betula papyrifera, balsam fir Abies balsamea, spruce (Picea spp.), pine (Pinus spp.), and red maple (Acer rubrum; Faber-Langendoen et al. 2007). Shrublands (1.2%), typically alder (Alnus spp.) and willow (Salix spp.), along with both woody (30.5%; largely tamarack Larix laricina and black ash Fraxinus nigra) and herbaceous wetlands (7.8%) are pervasive throughout the VNP ecosystem. These habitat classes represent significant variability in thermal properties (i.e., operative temperatures) at various times of day and year (Olson et al. 2014).

Wolves, the main predator of moose, are abundant inside and outside of VNP with ∼16–22 individuals in two to three packs inhabiting the centrally located Kabetogama Peninsula of the park alone (Olson and Windels 2014). However, wolves regularly make use of other prey (Gogan et al. 2004; Chenaux-Ibrahim 2015; Gable et al. 2016) such as white-tailed deer (∼3.8 deer/km2) and beaver (∼5.0 beaver/km2) that are present at much higher densities relative to moose (0.13 moose/km2; Windels and Olson 2016), albeit with much smaller body mass per individual.

Capture and handling of moose

From 2010 to 2012, we captured adult moose (males: n = 4, moose-years = 5; females: n = 9, moose-years = 14) during February and March using helicopters to dart individuals (Quicksilver Air, Inc., Fairbanks, AK, USA). We immobilized and anesthetized moose with 1.2 mL (4.0 mg/mL) carfentanil citrate and 1.2 mL (100 mg/mL) xylazine HCl, and used 7.2 mL (50 mg/mL) naltrexone HCl. We used 3 mL (5 mg/mL) yohimbine HCl as antagonist. During each moose capture, we extracted a blood sample for progesterone analysis. We classified individuals with progesterone values larger than 2 ng/mL as pregnant in a given year (see Schwartz et al. 1995; Schwartz 1998). We outfitted each moose with a GPS collar (Sirtrack Limited, Hawkes Bay, New Zealand). Fix attempts were scheduled at 15-min intervals during 2010 and 20-min intervals during 2011–2012. We estimated the average GPS error of locations from stationary Sirtrack collars at ∼7 m for a 50% circular error probable (McCann et al. 2016). Each GPS collar contained an accelerometer that provided an activity count on a 1-min average basis. All animal capture and handling protocols were approved by the University of Minnesota and National Park Service Animal Care and Use committees.

We removed any GPS locations associated with a mortality event, any fix with a horizontal dilution of precision (HDOP) >15, and checked that the final rate of movement was biologically feasible (fastest rate = 52.6 m per min). Additionally, we removed moose that did not have an accelerometer in the collar, and GPS locations collected after an accelerometer had become faulty or failed to record activity entirely.

HMMs require temporal regularization of the data (i.e., no missed GPS fixes), so we used the package waddle (Gurarie and Bracis 2013) in program R (R Core Team 2016) to linearly interpolate missing locations in the moose data for both the 15- and 20-min interval data independently. For our moose locations, we interpolated ∼6% of all locations for an average of 4.9% per moose-year. The waddle package also calculates the movement rate and relative turning angle for all movement steps. We refer to this as the full movement dataset.

We associated activity level with each movement step by binning accelerometry data based on the beginning and end timestamps from each movement step in the full movement dataset and averaging the corresponding activity values.

Spatial covariates

For covariates relying on the locations of moose, we used the coordinates at the end of each step from the full movement dataset. We determined habitat type by overlaying moose locations on a raster of the 2011 National Land Cover Database (Homer et al. 2015) and extracting the habitat classification of each raster cell using the package raster (Hijmans 2015) in program R (R Core Team 2016). Lakes and beaver-influenced wetlands within VNP were provided as ESRI Shapefiles of polygons (Johnston and Windels 2015b). We created a raster in ArcMap 10.3.1 (Redwoods, CA), which expanded beyond the boundaries of all moose locations, where each 30 m × 30 m cell contained the distance from the center of each cell to the nearest edge of a lake or pond. We overlaid all moose locations onto the distance-to-water raster in program R and extracted the distance value.

Weather covariates

We used weather conditions from a weather station located at the International Falls airport located ∼18 km from the eastern edge of VNP (weather station: 48.561, −93.398). The data, collected approximately once per hour, were available through the National Oceanic Atmosphere Administration’s National Centers for Environmental Information [March 1, 2016 (https://www.ncdc.noaa.gov/isd/data-access)]. For our weather variables of interest—ambient temperature (°C), wind speed (km per hour), and liquid precipitation (previous hour in centimeter)—we interpolated the values, using package zoo (Zeileis and Grothendieck 2005) in program R (R Core Team 2016), to create a timestamp match with our moose location data.

Temporal covariates

Using the timestamps associated with each moose location in the full dataset, we calculated the altitude of the sun in the sky at VNP using the package maptools (Bivand and Lewin-Koh 2017) in program R (R Core Team 2016). Values less than zero represent times of day where the sun was below the horizon, and positive values represent times where the sun was above the horizon.

The presence of snow cover, activities associated with calving, and the availability and quality of forage may strongly influence movement. Accordingly, our analyses considered three seasons that corresponded with declining snow cover and calving (spring, 1 April–30 June), the snow free period (summer/fall, 1 July–31 October; hereafter, “summer”), and full snow cover (winter, 1 November–31 March). Within season variation was considered by assigning a unique day number to each date within a season, starting with 1 on the first day of each season then increasing by one every day until the last day of the season.

Capture and handling of wolves

We captured adult wolves in VNP using padded-foothold traps (Livestock Protection Company, Alpine, TX, USA) from 2012 to 2014 during June–October. We fit each individual with either an Argos GPS (Telonics, Inc., Mesa, AZ, USA) or an Iridium GPS collar (Lotek Wireless, Inc., Newmarket, Ontario, Canada; Vectronic Aerospace GmbH, Berlin, Germany). The collars attempted a fix once every 20 min–6 h, but most were for every 4–6 h, for up to 2 years. We removed any location with a HDOP value >15.

We collected GPS data from 26 individual wolves between 9 October 2012 and 17 November 2015. We removed any locations that occurred more frequently than ∼4 h (3 h, 57 min) of the last included location to reduce the bias in home range delineation among packs. We pooled data among individuals that we visually determined to be in the same pack based on proximity in time and space of GPS locations. We assumed the GPS locations of the collars represented the pack’s location for the duration of the study. We split all pack data into seasons (same seasonal delineations as the moose) and only included pack-seasons that had a minimum of 30 fixes that met all of our criteria. We retained 18,367 total fixes from 20 distinct packs or individuals inhabiting unique areas of study for at least one full season. Seasonal home ranges contained 6,093 locations from 17 pack-seasons in winter, 1,704 locations from 6 pack-seasons in spring, and 7,087 locations from 17 pack-seasons during summer. Areal coverage of combined annual wolf home ranges contained most of the VNP ecosystem and overlapped the home ranges of GPS-collared moose.

We calculated seasonal home ranges for wolf packs with Program R (R Core Team 2016) using package adehabitat (Calenge 2006). We used the function kernelUD to create a bivariate normal kernel for each seasonal home range and used the ad hoc method for the smoothing parameter. We created a utilization distribution raster for each wolf pack’s seasonal home range where each raster cell contained the value of the smallest isopleth of the home range covering the cell. For each moose location in the full dataset, we extracted the values of all wolf pack utilization distributions calculated at that location and season. We binned the extracted values into 0–19, 20–39, 40–59, 60–79, 80–99, and >99 categories of the isopleth values of wolf home ranges (i.e., a value of 100 indicates that the moose location was not overlapped by a wolf home range, a value between 0 and 19 is within the very core of a pack’s home range). We assigned each bin a value of 0–5 with larger values corresponding with smaller utilization distribution values from the wolf pack home ranges. The largest bin value for each moose location across all wolf packs was retained for modeling. We refer to this value as “wolf home range” going forward.

Calibration of collar activity

We combined direct human observations of moose behavior with data from SirTrack GPS collars outfitted with accelerometers (the same make and model as the wild moose in VNP area) to determine what metric of activity count best corresponded to an observed moose’s behavior. From 30 July to 2 August 2009 an observer followed a single human-habituated moose inside a ∼2.6 km2 pen at the Kenai Moose Research Center (KMRC) near Soldotna, Alaska. The KMRC research facility is owned by the Alaska Department of Game and Fish and houses three to six moose per 2.6 km2 holding pen. Within each pen, the moose forage on naturally occurring vegetation and habitat primarily consists of a mix of mature white (Picea glauca) and black (Picea mariana) spruce along with deciduous trees such as aspen Populus tremuloides, willow (Salix sp.), and cottonwood Populus trichocarpa. Moose also had access to ponds, small lakes, and sedge meadows.

The observer classified the behaviors of the focal moose from a distance of 3–10 m. The time and a classification of behavior (walking, feeding, drinking, standing, standing and ruminating, bedded, or bedded and ruminating) that occurred for a minimum of 3 s were recorded. The moose was observed for a total of 2,315 observation minutes over four 7–11-h observation periods (see Ness [2010] for additional details).

The mean, maximum, and standard deviation of activity levels from the accelerometer, and the mode of the behavior classification were compared based on 10,000 sets of sub-sampled data. Activity data were sub-sampled in both 15- and 20-min time spans to correspond with the length of GPS fix intervals from wild moose in VNP using random starting points with replacement. All consecutive data following the randomized starting point for the given time span (i.e., the next 15 or 20 min of data for both activity counts and direct observational data) were then included in the analysis. We performed an ANOVA with a post-hoc Tukey’s Honest Significant Difference (HSD) test to determine significant differences in the accelerometer activity levels between the three behaviors resulting from the mode calculation of the direct observational data: feeding, ruminating, and moving. We considered these three mode behaviors to correspond to our three behavioral states in our HMM model. We tested which metric of activity levels (mean, median, max, or sd) provided the most significant differences among the three behaviors.

Activity levels from the accelerometer data of the captive moose were significantly different among the four summary metrics (mean, median, max, and sd) for the observed mode behavioral state of either feeding, ruminating, or moving (adjusted P < 0.001 and 95% confidence intervals of the Tukey’s HSD did not overlap zero for all metrics for both 15 and 20 minute intervals). Based on the ANOVAs, mean activity had the largest differences among the observed behaviors for both 15-min intervals (mean: F2,9998 = 4387, SD: F2,9998 = 2096, median: F2,9998 = 1676, max: F2,9998 = 552) and 20-min intervals (mean: F2,9998 = 6390, SD: F2,9998 = 2894, median: F2,9998 = 1978, max: F2,9998 = 765). Mean activity counts for the 15-min intervals were 13.9, 4.8, 0.67, and 13.1, 4.7, 0.75 for the 20-min intervals corresponding to movement, ruminating, and rest, respectively.

We used the mean activity counts from the captive moose behaviors as informative priors for the true mean activity level in our HMM model for the 20-min interval data. For the data corresponding to the 15-min location moose, we re-ran the same procedure on the same captive moose activity data using a minimum activity level of 0.01, and adding a random variance of 0.001 for each observation as was done with the VNP data to better fit the assumption of a normal distribution (there were a large number of low activity counts from the VNP moose). The resulting values were used for the informative priors in the HMM models for moose with 15-min GPS interval data.

Statistical analysis

Our approach to understanding moose behavior involved two distinct steps. First, we assigned each moose location from the full dataset a behavioral state, using a Bayesian three-state switching HMM. This HMM utilized step lengths and turning angles derived from sequential GPS locations, and activity counts from accelerometers within the same GPS collars to classify each location into one of three behavioral states. We did not use our HMM to test how weather conditions, spatial attributes of the landscape and habitat, and wolf home range areas affected moose behaviors because the large increase in number of parameters would have made processing times and synthesis prohibitive. Instead, we utilized a frequentist approach, generalized additive models (GAMs) with a multinomial distribution (categorical behavior classification was the response), to model the influence of these covariates on moose behavior (e.g., Russell et al. 2015).

HMM: state assignment

We expanded on code provided in the Supplementary Material in Gurarie et al. (2016) for an HMM three-state switching model. Our aim was to classify each moose location as one of three categorical latent behavioral states: 1) traveling, 2) foraging, or 3) resting. Following the assumptions of Gurarie et al. (2016), we used a wrapped Cauchy distribution to describe turning angles and a Weibull distribution to describe step lengths (Kareiva and Shigesada 1983). We assumed that traveling would be associated with larger travel distances and directional persistence (i.e., turn angles close to 0, straight-line movements). We assumed resting would be associated with the smallest travel distances and turning angles that were more evenly distributed compared with those of the travel state, but with a higher probability of being toward 180° (this is because all locations have some GPS error and stationary collars will yield “steps” with a tendency to have a tight turning angle back toward the actual location of the collar). The foraging state should be characterized by step-length distances commonly larger than the resting state, but smaller than the traveling state, and a generally uniform distribution of turning angles as the moose meanders through its habitat searching for and consuming forage or standing still while browsing during the 15- or 20-min period. To increase our ability to differentiate among the three latent behavioral states, we incorporated a third data source: activity data from the GPS accelerometer. We used a log-normal distribution for the activity parameter, and assumed that activity levels would increase from resting to foraging to traveling.

We used vague priors for all parameters except those of movement scale and activity. Specifying slightly informative priors (based on an exploratory analysis of a few individual moose and data from captive animals) helped to keep the parameters aligned for all of our MCMC chains (without specifying the mean prior values, the order of the behavior vectors would flip activity and travel). Other priors, such as those for movement shape, turn angles, and transition probabilities between states were unchanged from those of Gurarie et al. (2016). For the zero-truncated Gaussian priors of the movement scale parameters (one for each behavioral state: travel, forage, and rest), we provided mean values of 100, 40, and 10 for the 20-min data and 75% of these values for the 15-min data (75, 30, and 7.5); all of these priors had precisions set to 0.001. The informative priors for the mean of activity were derived from the captive moose data at both the 15- and 20-min intervals with a large variance term (variance = 1,000, JAGs precision value of 0.001). The prior for the variance of activity was uniform from 0 to 30 for all three states. We added a maximum step length to be considered in the rest state of 50 m. This value is larger than our GPS error estimates of ∼7 m for the 50% error, but given a time step of 15–20 min, we believed this would allow for a small amount of movement (e.g., moving to a nearby bedding site) along with GPS error. All step lengths and activity measurements less than 0.001 were set to 0.001. Because large gaps in the GPS data can cause model fitting issues when using the interpolated values to fill in missing locations (e.g., all turn angles are 0 and the step lengths are constant), we removed all data interpolated over large gaps. This removal was necessary for three moose-years of data where the collar had stopped acquiring fixes and re-started later in the year.

We fit the models using a Bayesian MCMC framework with JAGs through the RJAGS (Plummer 2016) and R2jags (Yu-Sung and Yajima 2015) packages. We ran three chains of 5,000 iterations with a burn in of 1,000 and we thinned the chains by eight resulting in 500 samples from each chain. We checked the Gelman–Rubin statistic, making sure it was close to a value of 1, and visually assessed the three chains of each parameter to judge convergence. If a model did not converge, we then ran additional iterations in intervals of 1,000. Typically, the initial 5,000 iterations were enough for convergence when proper initial values were provided for the initial movement and activity values. We used the mode latent state associated with the samples as the behavioral state assigned for each movement step. We report the mean and its associated variation for our parameters by randomly sampling, with replacement, 1,000 posterior point estimates from each moose-year’s HMM posterior distribution, calculating the mean of these samples (i.e., 1 random sample from each moose-year) and estimating the 95% credible interval of the distribution of means by parameter, behavior, and fix rate (i.e., separate for 15 and 20 min fix data).

To ensure the assigned behavioral states from the HMM-provided reasonable classifications, we utilized a suite of diagnostic plots to visually inspect that behaviors were assigned with expectations based on movement distance (i.e., the largest distances for travel and the smallest for rest), turning angles (i.e., most directional persistence for travel and the least for rest), and activity counts (i.e., largest activity counts for travel and the least for rest). A moose may exhibit multiple behaviors in the 15 or 20 min time of our GPS fix intervals, but the HMM approach allows us to assign a single behavior associated with each GPS fix interval that is the most likely primary behavior based on the distributions for movement rate, turning angle, and activity level. We determined the mean proportions of moose locations assigned each behavior by season and used a bootstrap to estimate the uncertainty among moose-years. The bootstrap estimated 95% confidence intervals (using package boot; Canty and Ripley 2016) based on the adjusted bootstrap percentile method (Davison and Hinkley 1997) using 10,000 bootstrap samples of the proportions of behaviors assigned to individual moose-years by season. Additionally, we calculated the average bout length (in minutes; a bout is defined by consecutive moose locations assigned the same behavior) per individual by behavior type, and generated bootstrap estimates of the predictive distribution.

Modeling behavioral states and covariates

Moose behavioral states were analyzed using a GAM with a multinomial distribution using package mgcv (Wood 2011) in program R (R Core Team 2016). The multinomial distribution allowed us to estimate how each covariate influenced the likelihood of moose being in either the resting or traveling states relative to the foraging state (i.e., changes in the likelihood of foraging are contingent upon changes in non-foraging behavior [traveling or resting]; e.g., a reduced likelihood of traveling with a constant likelihood of resting must result in more foraging activity). This estimated value is called the log odds ratio and it measures the change in the log odds of the probability of being found in one state against the baseline. Because the coefficient values are relative to the baseline state of foraging, a negative coefficient value for the rest state and a positive for the travel state suggest that as the covariate increases, moose are more likely to be found in the travel state and less likely to be found in the rest state relative to foraging. Additionally, GAMs enabled us to incorporate smoothers for covariates we believed to have non-linear relationships within the behavioral state of the moose. Smoothers offer a greater flexibility relative to regression splines, allowing us to estimate at what temperature moose reduce travel and increase rest, and estimate a general daily activity budget by including a smoother for the altitude of the sun.

We modeled each moose’s behavior on a seasonal level using GAMs. We standardized and centered covariates to compare the effect size of each covariate to one another and interpret how strong of an effect each covariate has on the likelihood of being in either of the rest or travel states. We modeled the behavioral state of each moose location reported from the HMMs as a function of smoothers for (1) ambient temperature and (2) altitude of the sun, and linear predictors based on the (3) categorical habitat classification (deciduous forest, mixed forest, conifer forest, shrub/scrub, emergent herbaceous wetland, woody wetland, and open water), (4) distance to lake or ponds, (5, 6) weather conditions (precipitation and wind speed), (7) day number of the season, and (8) wolf home range. In a few instances we removed the wolf home range covariate or a given habitat class from a seasonal individual moose GAM model if an individual did not have enough variation in either overlap in a known wolf pack’s range (e.g., a moose could not be found exclusively within the 40–60% isopleth level for a pack or the model could not run; locations were required within multiple isopleth values), or habitat type (e.g., for moose-years with relatively few locations in a season, a less common habitat type such as shrub, may have three locations that were all classified the same) because of the multinomial distribution of the model.

To summarize the results of the individual moose GAMs, as a way of providing a “population-level” result and to provide a measure of variability among individuals, we report a mean, based on the estimates of individual moose-year model coefficients and bootstrapped 95% confidence interval using the same bootstrapping method previously described. Relative effect sizes can be determined by comparing coefficient values among categorical and continuous variables because they are scaled and centered. These effects highlight trends in the associations between covariates values and behaviors from the HMM. We considered P-values at or below an alpha level of 0.05 to be significant for linear predictors and smoothers. To summarize the effects of the smoothed covariates (ambient temperature and sun altitude) across all individual moose models, we reported on the significance of the smoothed terms in the same manner as the linear predictors (i.e., based on the number of individual moose models out of the total). The predicted mean response and 95% confidence interval were then plotted by behavioral state. We calculated the mean and confidence interval with the same bootstrapping method previously described using the predicted values of the response at a given grouped level of the independent variable (e.g., increments of 3 °C).

Because we hypothesized that wolf home range might have a larger influence on the behavior of pregnant moose during the spring, when moose in VNP give birth, we re-ran our GAMs based on moose-year instead of individual moose. This enabled us to link the coefficient values for wolf home range of the resulting moose-year GAMs with classifications for: 1) pregnant moose, and 2) non-pregnant females and males. We performed an ANOVA with a post-hoc Tukey’s HSD test to determine differences between the groups.

Results

We classified 348,226 observations of moose movement and activity from 13 individuals (males = 4; females = 9) and 19 moose-years (males = 5; females = 14; parameter point estimates for HMM models in Online Appendix Table 1). Our spring (observations = 120,603, n = 11, moose-years = 17), and summer HMM models (observations = 87,847, n = 8, moose-years = 9) resulted in similar proportions of assigned behaviors; proportions of each behavior had overlapping confidence intervals, and within seasons, travel was assigned to significantly fewer moose locations than forage and rest ( [95% CI]; spring: forage = 0.44 [0.39–0.47], rest = 0.38 [0.33–0.42], travel = 0.19 [0.13–0.27]; summer: forage = 0.40 [0.39–0.42], rest = 0.42 [0.32–0.48], travel = 0.18 [0.12–0.29]). During winter (observations = 139,756, n = 13 individuals, moose-years = 19), moose were less active (0.45 forage or travel) relative to spring and summer (0.62 and 0.58 forage or travel, respectively) due to decreased travel (  = 0.07, 95% CI = 0.05–0.11) and increased rest ( = 0.55, 95% CI = 0.51–0.57). Across seasons, length of continuous behaviors (minutes per bout) was least for traveling ( = 66.8, 95% CI = 58.6–75.0), intermediate for foraging ( = 75.8, 95% CI = 62.9–88.8), and greatest for resting ( = 98.8, 95% CI = 79–118.6).

Weather conditions and temporal effects

Nearly all moose altered their behavior based on the altitude of the sun (i.e., amount of daylight; smoother P < 0.05) for both the travel (96.9% of moose) and rest states (96.9% of moose) across all seasons (32 moose-seasons total). During spring, but especially during summer, moose were more likely to rest near-midday and the middle of the night, and were more likely to travel and forage during crepuscular times (Figure 1). Moose activity patterns were nearly reversed during winter when rest occurred much more frequently in the dark, and moose foraged and traveled during daylight hours (Figure 1). Within seasons, moose exhibited changes in behavior resulting in large effect sizes, by foraging more (less rest and travel) as winter turned into spring, traveling and resting more (less foraging) as spring progressed, and reducing travel between the start and end of summer (see Day # of season; Table 1; summaries of coefficient direction [±] and significance in Online Appendix Table 2).

Figure 1.

Predicted values of the relative-risk ratios from smoothers modeling the influence of solar elevation (i.e., amount of daylight; −2 = dark/nightime and 2 = brightest/noon) on the resting behavior of GPS-collared moose in VNP, Minnesota, USA. Each moose’s behavioral state, classified using a HMM, was modeled individually on a seasonal basis using a GAM with a smoother for solar elevation for the behavior of rest relative to foraging. Here, the predicted log odds ratio represents the likelihood of choosing one outcome category (rest) relative to the baseline category (forage). Standardized and centered values of solar elevation were pooled by increments of 0.2 to predict the mean and 95% confidence interval for all pooled individual moose estimates using bootstrapping.

Table 1.

Mean (±95% bootstrapped confidence intervals) of estimated coefficients from generalized additive model (GAM) results of individual moose-year data by covariate and season

Variable	Season	n	Rest	Travel
Day number of season	Spr	11	0.16 (0.06–0.39)*	0.25 (0.07–0.46)*
	Sum	7	0.03 (−0.15–0.12)	−0.24 (−0.37–0.08)*
	Win	13	−0.11 (−0.31–0.02)*	−0.39 (−0.91–0.08)*
Lake/pond distance	Spr	11	0.03 (−0.01–0.07)	−0.14 (−0.18–0.09)*
	Sum	7	0.09 (0.02–0.22)*	0.03 (−0.09–0.16)
	Win	13	−0.01 (−0.07–0.03)	−0.23 (−0.81–0.06)*
Precipitation	Spr	11	−0.02 (−0.06–0.02)	0.05 (0.01–0.1)*
	Sum	7	0.01 (−0.07–0.22)	0.06 (−0.02–0.28)
	Win	13	−0.03 (−0.06–0.01)*	0.03 (0–0.06)*
Wind speed	Spr	10	0.07 (0.04–0.14)*	−0.02 (−0.07–0.07)
	Sum	7	−0.03 (−0.06–0.01)*	−0.08 (−0.17–0.01)*
	Win	13	0.02 (0.00–0.04)*	−0.02 (−0.15–0.05)
Wolf home range	Spr	8	0.04 (−0.02–0.11)	−0.12 (−0.22–0.03)*
	Sum	7	−0.03 (−0.05–0.01)*	−0.03 (−0.1–0.03)
	Win	12	−0.01 (−0.04–0.02)	−0.10 (−0.68–0.11)
Emergent herb. wetland	Spr	11	0.09 (−0.16–0.21)	0.35 (0.06–0.61)*
	Sum	7	0.32 (0.14–0.55)*	0.13 (−0.45–0.38)
	Win	13	−0.14 (−0.62–0.07)	0.65 (0.19–1)*
Open water	Spr	9	−0.25 (−0.94–0.03)	0.74 (0.33–1.38)*
	Sum	7	−0.24 (−1.07–0.1)	0.12 (−0.19–0.42)
	Win	7	−0.2 (−0.63–0.14)	1.61 (0.97–2.28)*
Woody wetland	Spr	11	0.12 (0.04–0.21)*	0.04 (−0.22–0.22)
	Sum	7	0.25 (−0.01–0.63)	−0.07 (−0.58–0.18)
	Win	13	0.13 (0.05–0.2)*	0.47 (0.3–0.67)*
Shrub/scrub	Spr	9	0.10 (−0.06–0.3)	0.15 (−0.14–0.59)
	Sum	6	−0.18 (−0.33–0.03)*	−0.15 (−0.33–0.11)
	Win	9	−0.22 (−0.30–0.08)*	−0.36 (−0.85–0.00)*
Mixed forest	Spr	11	0.03 (−0.03–0.1)	−0.09 (−0.25–0.07)
	Sum	7	0.07 (−0.05–0.19)	0.08 (−0.1–0.36)
	Win	12	0.07 (0.02–0.13)*	0.12 (−0.05–0.26)
Evergreen forest	Spr	11	−0.11 (−0.6–0.16)	−0.14 (−0.36–0.09)
	Sum	7	0.06 (−0.49–0.49)	0.19 (−0.04–0.48)
	Win	11	0.13 (0.06–0.21)*	0.00 (−0.46–0.46)

Notes: Behaviors (travel, forage, and rest) of GPS-collared moose in VNP, Minnesota, USA, were classified using a three-state HMM with switching. These behaviors were analyzed using GAMs with a multinomial distribution by season. Habitat cover types (last six covariates listed) were categorical and relative to the cover type for deciduous forest. Results for the behavioral states of rest and travel are relative to the foraging state (log odds ratio). Sample size was based on unique moose-year.

95% bootstrapped confidence interval did not overlap zero.

Beyond the influence of the sun’s altitude, moose behavior changed in association with variation in ambient temperature (Figure 2). Nearly all moose altered their behavior in association with changing ambient temperature throughout the seasons (smoother based on temperature: % of moose for which P < 0.05; winter: travel = 84.6%, rest = 92.3%; spring: travel = 100%, rest = 100%; summer: travel = 85.7%, rest = 85.7%). The one individual who did not significantly respond to temperature during summer had considerably less data relative to other moose (1,741 fixes or ∼24.2 days of fixes; mean fixes per moose during summer for all moose= 10,980). In response to warm temperatures, moose reduced their travel in both spring and summer (Figure 2) and correspondingly increased resting behavior. The change from a positive association between travel and ambient temperature to a significantly negative one occurred at a higher ambient temperature during summer (24–27 °C) relative to spring (12–15 °C; Figure 2).

Figure 2.

Predicted values of the relative-risk ratios of smoothers modeling the influence of ambient temperature (°C) on the travel behavior of GPS collared in VNP, Minnesota, USA, during spring and summer. Each moose’s behavioral state (rest, forage, or travel) was classified using a HMM. The behavioral states were modeled individually on a seasonal basis using GAM with a smoother for ambient temperature. The log odds ratio represents the likelihood of choosing one outcome category (travel) relative to the probability of choosing the baseline category (forage). The predicted log odds ratio predictions were pooled by increments of 3 °C, and the mean and 95% confidence interval were calculated for all pooled individual moose estimates using bootstrapping.

Precipitation and wind speed did not have strong effects on moose behavior, but moose consistently responded to changes in both (Table 1; Figure 3). Moose generally increased resting behavior with increased wind speed, especially during spring, and increased foraging when wind speed increased during summer (Table 1; Figure 3). In contrast, moose were more likely to be active (i.e., traveling or foraging) during precipitation events especially during spring and summer (Table 1; Figure 3).

Figure 3.

Coefficient values from GAMs showing the influence of variables on the behavior of GPS-collared moose in VNP, MN, USA. Moose behaviors, as classified by a HMMs, were modeled by individual moose and season using GAMs with multinomial distributions that provide the log odds ratio estimate. The log odds ratio estimates the likelihood of choosing one outcome category (rest or travel) relative to the probability of choosing the baseline category (forage). Each point represents an estimate from a single moose GAM for each given behavior. The model coefficients for the resting behavioral state response during the A) spring, B) summer, and C) winter, and the traveling behavioral state response during the D) spring, E) summer, and F) winter are shown. Covariate coefficients of log odds ratio estimates were classified as significant if the P-value was ≤0.05.

Habitat and landscape effects

Based on the effect sizes of the cover types, moose behavior was strongly influenced by habitat cover types and landscape features throughout the seasons, albeit with a large degree of variability among individuals (Figure 4). Moose used wetland habitats for travel and rest more than for foraging throughout the year (Table 1; Figure 4). Moose were rarely located in open water (% of locations per moose; spring: = 0.1, range ≤ 0.1–2.6; summer: = 1.7, range = 0–5.5; winter: ≤ 0.1, range = 0–0.1). But when they did use these areas—which likely consist of open areas along the frozen shoreline in winter and spring—it was consistently for travel (Table 1, Figure 4). Distance to the nearest lake or pond had relatively large effects on moose, but the behaviors associated with these areas changed throughout the year (Figure 3). In the spring, most moose reduced travel and reduced rest (Table 1, Figure 3A,D) in areas further from lakes and ponds, which suggests these areas are primarily used for foraging. The opposite was true for summer as moose significantly increased resting and traveling further from shore (travel not significant [Table 1] but most moose had positive coefficients [71%; Figure 3B,E;Online Appendix Table 2]). During winter, moose decreased travel when closer to shore (Table 1; Figure 3F).

Figure 4.

Coefficient values from GAMs showing the influence of habitat types on the behavior of GPS-collared moose in VNP, MN, USA. Moose behaviors, as classified by a HMMs, were modeled by individual moose and season using GAMs with multinomial distributions that provide the log odds ratio estimate. The log odds ratio estimates the likelihood of choosing one outcome category (rest or travel) relative to the probability of choosing the baseline category (forage). Each point represents an estimate from a single moose GAM for each given behavior. The model coefficients for the resting behavioral state response during the A) spring, B) summer, and C) winter, and the traveling behavioral state response during the D) spring, E) summer, and F) winter are shown. Covariate coefficients of log odds ratio estimates were classified as significant if the P-value was ≤0.05.

Moose used evergreen/conifer forests primarily for rest during winter, and while there was a large degree of variability among individual moose-years, 43% of moose-years had a significant positive association between rest and evergreen/conifer forests during summer (Table 1, Figure 4; Online Appendix Table 2). Mixed forests were also primarily associated with resting behavior of moose during winter and summer and generally with reduced travel during spring (Table 1; Figure 4). However, while the influence of mixed forests was consistent across seasons, the influence relative to deciduous forests was not large (Figure 4). Moose used shrub/scrub habitats primarily for foraging during summer and winter (i.e., moose consistently reduced travel and rest in scrub/shrub; Table 1; Figure 4).

Wolf home ranges and their effects on moose

Wolves used the smallest area during the spring (n = 6; 50% Isopleth = 25.3 km2, 95% CI: 11.9–38.7 km2; 95% Isopleth = 130.4 km2, 95% CI: 87.3–173.5 km2) resulting in the lowest percentages of moose locations within known wolf home ranges (25.3% of all locations, 72.7% of moose with any spatial overlap). A high percentage of moose locations that did overlap with wolf home ranges was primarily in the outermost area of wolf home ranges (80–99% isopleth; 47.5% of locations with overlap) and no moose was located within the 40% isopleth value of any wolf pack’s home range. A majority of moose (n = 8, moose years; 63%) located in these outer areas of wolf pack home ranges increased resting behavior closer to wolf pack cores (Online Appendix Table 2) and the same percentage of moose reduced traveling in these areas substantially (Table 1; Figure 3D). During the spring, females that were pregnant at the time of capture (n = 5 moose-years) reduced their travel more than non-pregnant females and males (n = 7 moose years) in areas closer to wolf pack cores (pregnant vs. other coefficient for wolf home range = −0.16). The largest reductions in travel were associated with three of the five moose-years for which the female was pregnant that year (β = −0.25, −0.26, −0.35). However, the small sample size did not result in a significant difference between the two groups (Padj. = 0.08, 95% CI pregnant vs. other = −0.35–0.02). Females that tested pregnant were only slightly less likely to rest when located closer to wolf home range cores (difference pregnant vs. other = −0.07, Padj. = 0.26, 95% CI pregnant vs. other = −0.19–0.06) suggesting an increase in foraging or another behavior with similar behavioral characteristics.

During summer, wolves expanded their home ranges (n = 17; 50% Isopleth = 42.1 km2, 95% CI: 28.2–56.0 km2; 95% Isopleth = 170.6 km2, 95% CI: 112.7–228.5 km2) resulting in all moose locations overlapping at least one wolf pack’s home range, and nearly half of locations falling in the inner <40% isopleth core (48.2% of locations). Moose in close proximity to the core of wolf home ranges reduced rest (Table 1), but the effect was small (Figure 3B).

Wolf home ranges were largest on average during winter months (n = 17; 50% Isopleth = 83.6 km2, 95% CI: 83.6–145.6 km2; 95% Isopleth = 345.1 km2, 95% CI: 71.1–619.1 km2). Wolf packs overlapped all moose locations again, and the majority of moose (76.9% of individuals) were located in areas closest to the core of the wolf home ranges (<40% isopleth) during at least some portion of the winter. Most moose, with adequate variation in wolf home range, responded to the closer proximity to wolf home range cores by altering their amount of travel (58% of individuals with significant coefficients; Online Appendix Table 2). While some moose had exhibited large effects in travel behavior from wolf home range cores (Figure 3C,F) the response was not consistent (50% increased travel and 50% decreased travel).

Discussion

Our approach, utilizing HMMs to classify the behavior of moose, provides additional insights about moose living along the southern extent of their bioclimatic range beyond that of traditional habitat selection studies. By incorporating movement characteristics and activity counts from accelerometers (validated by direct observations on captive moose), we were able to distinguish between three behaviors of moose remotely, at fine temporal and spatial scales, and associate the changes in behavior with different environmental variables. While individual moose exhibited variable behavior associated with changes in weather, landscape, and the presence of predators, we were able to identify several consistent behavioral responses, especially regarding time of day and ambient temperature. Our findings highlight the difficult situation facing thermally sensitive species coping with climate change; rising temperatures will require individuals to reduce heat production and heat exposure by resting more frequently or shifting activity to nighttime when temperatures abate, thus reducing opportunities to seek out new food sources and forage, especially when predators are present.

Assessing thermal thresholds for moose can be difficult outside of laboratory or captive settings, but our approach, which accounted for the influence of solar elevation, time of year, and habitat, characterized moose behavior at a fine scale and captured the ambient temperatures when moose became less active (i.e., travel) and more likely to rest. The thermal thresholds of moose found by Renecker and Hudson (1986), who studied captive moose living in enclosed pens, were accepted as the best estimates of heat stress for many years. They found that moose increased respiration rates at 14 °C and started open-mouth panting at 20 °C during the summer. Recently, McCann et al. (2013) found that heat stress thresholds vary for moose housed in an outdoor enclosure and that shade and wind speed influenced the thresholds. With no wind, moose increased respiration rates for evaporative cooling at 17 °C, while with wind respiration rates did not increase until 24 °C. Similarly, free-ranging moose in VNP began resting more frequently in summer when temperature reached 21 °C and reduced travel when they reached 24 °C. While we were able to capture strong behavioral changes at these ambient temperatures, without the ability to directly observe the moose or take physiological measurements, we can only state that we estimated behavioral-related changes in activity and not necessarily heat stress. A 24 °C threshold was also corroborated by Broders et al. (2012), who found that moose in Nova Scotia sought thermal shelter when daytime temperatures reached 24 °C. However, it is also important to recognize that the thermal thresholds for moose are likely to vary from region to region, and individuals with poorer health, as a result of disease or parasite loads, may exhibit a lower threshold of thermal tolerance, as was indicated for populations near VNP (Murray et al. 2006; Lenarz et al. 2009) and for captive moose in Minnesota (McCann et al. 2013).

Spring thermal behavioral thresholds affecting moose were much lower than summer thresholds. Despite the rapid increases in spring temperatures observed globally, especially in the Midwestern United States (Schwartz et al. 2006), they often do not receive as much attention as the more extreme values of ambient temperature forecasted for summer months. We expected moose to exhibit a lower spring behavioral thresholds because seasonally shifting thermal tolerance is well established across many taxa of animals (Pörtner 2002). Our findings, which suggest a spring thermal behavioral threshold around 15 °C, may be a result of natural acclimation to cold winter temperatures and the rapid transitions in weather and phenology that occur during spring, exacerbated by remnant winter coats. It could also reflect a lack of canopy and horizontal cover within deciduous vegetation at this time. Research by Lenarz et al. (2009) highlighted the importance of spring temperatures by associating warmer spring ambient temperatures with lower adult survival in Minnesota (although see Mech and Fieberg [2014] which questions this association). However, earlier and warmer spring temperatures can also benefit moose by making spring forage more abundant and available earlier in the year. Grøtan et al. (2009) linked these potential benefits of warmer spring time temperatures experienced by moose in Norway to subsequent increases in calf survival, while Monteith et al. (2015) found a negative effect on recruitment associated with warmer spring temperatures due to a mismatch in phenology of early season forage and nutrition. While the overall impacts of a warmer and earlier spring season are still debated, other environmental conditions, such as more frequent precipitation events, may help counteract times of heat stress more so during spring than summer.

Temperatures are forecasted to increase further in the coming decades, which means that moose will likely need to rest more frequently and thereby forego travel and foraging to avoid heat stress. Galatowitsch et al. (2009) forecasts that the VNP area will experience an increase in the average daily minimum and maximum summer (June–August) temperatures of around 1.6 °C by 2030 and 3.3 °C by the year 2060. Studies by Dussault et al. (2004) and Street et al. (2015) found that moose coped with warm summer daytime temperatures by reducing daytime activity and switching to nocturnal activity. Our findings show that moose in VNP follow a similar seasonal pattern of resting during the heat of the mid-day in summer months, and decreasing rest at night. Moose responded in a similar, but weaker pattern during the spring, and reversed this pattern during the winter by becoming active during daylight hours. This behavioral plasticity may help moose cope with rising temperatures for years to come; however, it is clear that moose will be faced with a difficult tradeoff between energy acquisition and resting to deal with higher ambient temperatures. Street et al. (2016) examined this inherent tradeoff and found that a more northerly population in Ontario, Canada strongly selected for areas of better forage while the more southerly population in Minnesota, USA selected for a balance between habitats associated with better thermal cooling properties and those with better forage. Average wind speeds are also forecasted to be reduced under some climate change scenarios (e.g., Pryor et al. 2009), thus also impacting another primary source of cooling for moose (McCann et al. 2013).

Climate change in the region is already shifting vegetative communities in boreal ecosystems (Soja et al. 2007) and the pace of this change is expected to accelerate (Galatowitsch et al. 2009). This shift may eliminate the ability of moose to cope with heat by selecting for habitats with the best thermal properties. Numerous studies have found that moose select for habitats that provide better shade and cooling when ambient temperatures increase (Renecker and Hudson 1989; Schwab and Pitt 1991; Demarchi and Bunnell 1995; van Beest et al. 2012; Street et al. 2015, 2016; McCann et al. 2016). These studies generally agree that lowland forests with dense canopies are the most preferred by moose during times of high ambient temperatures.

Although some studies have found little difference in the abilities of different habitats to provide thermal shelter (Lowe et al. 2010), Olson et al. (2014) found large differences in thermal properties among habitats in the VNP ecosystem that drive the operative temperatures experienced by moose in the various habitat types. Open habitat types, such as scrub/shrub and wetlands had an average difference in afternoon temperatures during the summer of 3.38 °C (and a maximum of 8.10 °C) relative to forested habitats. These findings support the patterns found in our results that moose forage in mixed forests and shrublands, and travel through wetlands—both active behaviors—but likely only when temperatures fall below the maximum thermal thresholds of a given season. Moose select for better thermal cover in conifer and deciduous forests with dense canopies when the thresholds are exceeded—a finding that was supported by Street et al. (2016) using step-selection functions to analyze fine-scale habitat selection elsewhere in Minnesota, and an analysis of moose bed site selection by McCann et al (2016). We did not include an interaction between ambient temperature and habitat type in our models due to the large increase in the size and complexity of the results of our models, especially while incorporating the effects of smoothers.

The overall lack of a strong, consistent response of moose to wolves in VNP is not surprising for both methodological and biological reasons. However, while the overall effect size of wolves on the behavior of all moose was relatively small, we did find some evidence that pregnant females close to the expected calving season (i.e., before and after presumed parturition) were less likely to be traveling in areas more frequented by wolves relative to adult females without calves and males. This response is likely a way to reduce the exposure of calves to predation risk, but our results also suggest that pregnant females increased foraging behavior in these same areas. Increased foraging in areas with a higher likelihood of attack seems counterintuitive, but our HMM classification scheme was limited to three basic behaviors. By relying solely on movement and activity data, we potentially misclassified some vigilance behavior as foraging. These behaviors may appear relatively similar in terms of movement and activity characteristics. Our inability to identify vigilance behavior may be exacerbated because wolf GPS data only temporally overlapped with the collection of moose GPS locations for 1 year. Even with better temporal and spatial overlap, it can be difficult to identify responses of prey species to predators when relying on GPS data alone (e.g., Eriksen et al. 2011). Biologically, studies in other parts of Minnesota have found the recently recovered wolf population to play a major role in moose population abundance (Mech and Fieberg 2014). However, studies in VNP, where wolf and moose population levels have been fairly constant over the last decade (Olson and Windels 2014), have found that moose comprise a relatively small percentage of wolf diet (∼2–3%; Gogan et al. 2004; Chenaux-Ibrahim 2015; Gable 2016). Wolves in VNP can avoid the task of predating moose because less-risky prey alternatives such as white-tailed deer and beaver are present at relatively higher densities than moose (Windels and Olson 2016) and consistently make up large portions of wolf diets within VNP (Gogan et al. 2004; Chenaux-Ibrahim 2015; Gable 2016).

While the direct influences of wolves may have been muted due to methodological hurdles, moose patterns of forage and rest in relation to lakes and ponds may be an indirect function of predator avoidance and thermal cooling. Lakes in many regions of moose range are typically frozen well into spring, and thus do not provide forage opportunities until summer. Moose may utilize the cooling properties of aquatic habitats during warmer months, and use the often forage-rich shorelines, wetlands, or shallow lakes and ponds in the region, while quickly dissipating heat (Schwab and Pitt 1991). However, these same areas along the shore may be the riskiest areas of the landscape during winter and early spring when wolves often select for frozen surfaces and shorelines (Kuzyk et al. 2004) to travel and hunt (Mech 1991). In a nearby study on Isle Royale, Michigan, USA, Montgomery et al. (2014) found a strong negative association between the distance to shore and successful predation events of wolves on moose. Wolves in VNP forage near water bodies and wetland habitats to prey heavily on beavers in summer (Gable et al. 2016). Our findings of VNP moose potentially support the idea of a tradeoff between forage intake, thermal cooling, and predator avoidance whereby VNP moose may risk being around shorelines where wolves frequent in order to access forage-rich habitat patches, such as wetlands created by beavers (Johnston and Windels 2015a), that also allows them to stay cool during the summer months, but utilize areas further from lakes in the winter and spring when water bodies are frozen. Morris (2014) reviewed studies that reported on how moose used aquatic habitats and found little support for the hypothesis of minimization of direct predation risk and heat stress amelioration. Instead, the analysis found greater support for maximizing nutrition and avoidance of biting insects. More research, possibly incorporating novel technologies and methodologies, may help researchers better determine how and why moose use aquatic habitats, beyond travel, forage, or rest, so that we may better understand how or if moose can cope with future drier environments.

Methods that enable researchers to more fully utilize the data from GPS units can offer new insights beyond habitat selection by discerning how an animal responds behaviorally to their environment (Edelhoff et al. 2016; Gurarie et al. 2016). Activity data from collar accelerometers are often unused or are analyzed in relative isolation from animal movement. Our HMM enabled us to combine activity and movement data to more accurately separate moose behavior among resting and foraging. Anecdotally, the HMM model struggled to assign behavioral states to large numbers of moose locations in a way that made biological sense when utilizing only step lengths and turn angles. Additionally, by utilizing a captive study of moose wearing identical GPS collars (including accelerometers) to most of the moose in our study, we were able to provide biologically appropriate priors that helped the model converge more quickly and give more sensible results. However, we caution that there are several factors to consider when utilizing a three-state HMM model with activity data: 1) several accelerometers failed or provided unreliable and infrequent measurements later in the lifespan of the collar and thus reduced our sample size, 2) Moen et al. (1996) found the association between directly observed behaviors and activity data were no longer reliable at GPS fix rates of >1 h due to multiple behaviors occurring within that time frame (our fix attempts were 15 or 20 min), and 3) the relationship we observed in wild moose between movement rate and activity had a large degree of variability among moose and even between moose-years of the same individual. Our behavioral classification might have been improved if were able to directly observe moose behavior in different habitats, and create unique priors for the mean of activity dependent on the habitat associated with the location of the moose. For instance, movements through thicker vegetation associated with certain habitats may result in higher activity counts than the same movement rate through a generally open habitat.

To further understand how moose along the southern edge of their range will respond to climate change, additional research may be needed to understand the effects of climate change on the quality of browse available to moose. Moose make habitat selection decisions by balancing the tradeoff between thermal cooling properties and forage (Street et al. 2016), but finding this equilibrium may become extremely difficult if climate change also reduces nutrition of available forage. Findings by Dearing (2012) and Kurnath et al. (2016) show that plant secondary compounds in herbivore forage may have temperature-dependent toxicity properties; as temperatures rise, these toxic properties may reduce the nutritional content of forage currently assumed to be important for moose diets. As an additional potential avenue of research, biologger technologies that allow for the remote measurement of physiological changes (e.g., abdominal and vaginal implant transmitters and cardiac biologgers) offer the next step in understanding how moose respond to stressors. When biologger technologies are combined with traditional GPS units and accelerometers, HMMs could identify other unique behavioral and physiological states that may not have been apparent with GPS and activity data alone (e.g., Ditmer et al. 2015). Additional states such as “resting and heat stressed” or “traveling and stressed” potentially from a predator would greatly enhance our understanding of both animal behavior and physiology.

Understanding how populations inhabiting areas along the geographic edges of their distribution mechanistically respond to changes in their environment can help predict future range shifts, and inform conservation or management needs. Analytical approaches that provide insights into how and why species alter their behavior in response to environmental change at fine spatial and temporal scales are a critical first step to understanding the potential for a population to persist. Our findings suggest that in the coming 50−100 years, the adaptability of moose inhabiting the VNP ecosystem will be tested in the face of a shifting bioclimatic range (Galatowitsch et al. 2009), whereby moose will either need to find more or higher quality forage in habitats used for thermal refuge, or risk the physiological impacts of heat stress. Our study highlights this predicament for moose in VNP and other southerly moose populations by estimating the seasonal maximum behavioral thermal thresholds and determining how moose altered their behavior based on habitat, landscape features, and risk of predation. However, moose are an adaptable species and have shown persistence in other southern extents of their range with even less vegetative cover (e.g., southern Manitoba, Canada and North Dakota, USA). While our study provides a closer look at the behavioral response to many aspects of their changing environment, future research should both expand upon new technologies and methodologies to better discern animal behavioral responses while more holistically incorporating other critical factors that may influence population persistence such as disease, parasite loads, and interactions with new competitors and predators.

Supplementary Material

Supplementary material can be found at https://academic.oup.com/cz.

Click here for additional data file.