The WAH utilizes over 350,000 km² of northwestern Alaska, typically migrating from wintering areas in the south, which vary by year and individual, to the calving ground and summer range in the north (Figure 2; Lent, 1966; Dau, 2015, Joly and Cameron, 2019). Calving generally occurs May 31–June 13 (Cameron et al., 2018). Beginning in 2009, GPS collars (Telonics, Mesa, AZ) were deployed annually on adult female caribou (≥2years old) as they swam across the Kobuk River during fall migration (Dau, 1997; Joly et al., 2012). Captures were conducted using procedures approved by the State of Alaska Institutional Animal Care and Use Committee (IACUC; 0040-2017-40). Collars were programmed to record locations every 8 h and by 2017, 203 collars had been deployed. From 2003 to 2016, the herd decreased from a high of 490,000 to 201,000 caribou and then increased to 244,000 in 2019 (Alaska Department of Fish and Game, 2019). Caribou populations are known to fluctuate at decadal time scales, and this is generally linked with large-scale climate patterns (Gunn, 2003; Joly et al., 2011). The northeast extent of the WAH range overlaps with the neighboring Teshekpuk Herd, and individuals between the two herds have been known to mix (Mager et al., 2013; Prichard et al., 2020).

We applied two different approaches to infer calving events from the 2010–2017 GPS data: an individual-based method and a population-based method (DeMars et al., 2013; Cameron et al., 2018). The former fit two a priori movement models (parturient and non-parturient) to movement rate data and model fit was evaluated using information criteria. The second method established a herd-specific movement rate threshold for calving from known events and then analyzed movement rates for individuals that dropped below the threshold using a 3-day smoothing parameter (DeMars et al., 2013). Using instances when the two methods agreed resulted in accurately classifying calving events 89% of the time (n = 119) when compared to aerial observation data (Cameron et al., 2018).

For data spanning 2010–2015, we used the calving events as reported in Cameron et al. (2018), in which aerial data were used to validate identified calving events from the movement-based approaches. For the data spanning the calving period of 2016–2017, we followed the procedures outlined in Cameron et al. (2018) to identify calving events without relying on supporting aerial data. However, because there was a record number of non-migratory individuals during the winter of 2016–2017 (Joly and Cameron, 2019) and the individual-based method is ill-suited for individuals not exhibiting migration movements prior to calving (Cameron et al., 2018), we incorporated a designation of migratory and non-migratory for each individual and adjusted the analysis as follows. For individuals that migrated (identified as crossing at least one of the three major rivers separating summer and wintering areas), we used the calving events from instances of method agreement. For individuals that did not migrate to a southern wintering area that year and that the two model results disagreed, we used calving events identified by the population-based method. For calving sites, we used the GPS location that corresponded with the identified calving event from the population-based method, because the individual-based method appeared to label events one GPS interval early (Cameron et al., 2018).

To address our first goal of spatial trends in calving areas, we defined an annual calving area as the area used by the majority (>80%) of individuals for calving in the herd in a given year (Gunn and Miller, 1986). We calculated a kernel utilization distribution (Worton, 1989) based on the calving sites for each year using the package “adehabitatHR” version 0.4.14 (Calenge, 2006) in the R statistical program version 3.4.3 (R Core Team, 2017). In this approach, a bivariate normal probability distribution is centered over each calving site in a given year and averaged together, resulting in one distribution (the kernel) for each year. Kernels were generated using a 500 × 500 m grid in an Albers equal area projection, which minimizes distortion along the latitudinal gradient given the relatively high latitude of our study area and ensured valid comparisons between years (Snyder, 1987). All kernels were generated using the same bandwidth smoothing parameter (h = 25,000) and we used the 95% contour as they resulted in unbroken range delineations for all years (Hooten et al., 2017). This approach, which is based on the explicit calving sites, minimizes potential bias in delineating calving areas that can be introduced by mismatches between calving timing and aerial observation timing during traditional surveys (Gunn and Miller, 1986).

To test the null hypothesis that the spatial distribution of annual calving areas did not vary by year, we employed a kernel randomization analysis outlined by Breed et al. (2006). For comparisons between 2 years, we randomly assigned (without replacement) a year designation to each calving site. Then, kernels for both years were generated using the same grid and smoothing parameters as outlined above. The area of both randomized kernels was then computed, as well as the area of overlap between the two kernels. Last, we calculated the test statistic as the area of the kernel overlap divided by the largest area of the two kernel regions. We repeated this process 250 times without duplicating any random year assignments. The p-value was calculated as the proportion of random overlaps smaller than the observed overlap for the 2 years being considered, so that if the observed overlap was smaller than all observed overlap values, the p-value was <0.004 (see Supplementary Figure 1 for illustration). We performed this analysis for all combinations of annual comparisons, ranging from sequential up to 7-year intervals, and considered our alpha level as 0.05 for a one-tailed test.

Our other goals were to understand the biotic and abiotic factors driving caribou calving site selection at the landscape level and the navigational mechanisms caribou employ to arrive there. We performed a resource selection function analysis (RSF; Manly et al., 2002) using the calving sites each year and compared them with random locations from the herd’s range, representing the third-order of selection (Johnson, 1980). To define range-wide availability, we drew a 100% minimum-convex polygon, constrained to the coastal boundary, around all GPS locations during the study period. Defining availability is a particular challenge for resource selection studies, with the implicit assumption that available points are unused and available to all individuals (Keating and Cherry, 2004; Aarts et al., 2008). We focused on a range-wide scale for this analysis because individuals in the herd used the polygon area throughout the 8 years of study and we detected calving events at the extreme southwestern and northeastern extent of the range, far outside of the traditionally defined calving area. For each of the 8 years from 2010 to 2017, we created 10,000 random locations within the polygon, for a total of 80,000 available points.

We attributed both used and available points with a combination of physiographic attributes and annually varying environmental indices. We attributed elevation values from a 5 m resolution digital terrain model derived from the Interferometric Synthetic Aperture Radar (U.S Geological Survey, 2017) and calculated a solar radiation index (Keating et al., 2007) for each point using slope and aspect derived from the terrain model. This index ranges from −1 to 1, with low index values corresponding to north-facing steep slopes, high values south-facing steep slopes, and flatter slopes around 0.35. We calculated a vector ruggedness measure (VRM; Sappington et al., 2007), which is a measure of the ruggedness of the terrain, for each point using the digital terrain model and a 15 × 15 m swath. We used a land cover classification map (Boggs et al., 2016) to attribute all points with land cover type and reduced the classifications into four categories from the original 20 based on diet categories of the predominant vegetation (forest, shrub, herbaceous, and lichen/sparse; Supplementary Table 1). We filtered points that occurred in pixels originally categorized as bare ground, fire scar, ice/snow, and water.

For environmental indices, we attributed the annual snow off date (day of year) specific to that year for each point as determined from Moderate Resolution Imaging Spectroradiometer (MODIS) data (Macander et al., 2015). We included two measures of primary productivity at multiple time intervals using the normalized difference vegetative index (NDVI, for review see Pettorelli et al., 2005a) acquired from the MODIS V6 and compiled into 7-day composites with 250 m resolution (data available from the U.S. Geological Survey; Jenkerson et al., 2010). For an index of forage quantity, we used the raw NDVI value at a weekly temporal resolution and as an index of forage quality we calculated the change in NDVI values between sequential NDVI composites (NDVI rate) for the same time period by calculating the difference between sequential coverages (denoted “ΔNDVI”). ΔNDVI has been used in prior studies as an index for forage quality, including in Africa (Boone et al., 2006) and Alaska (Griffith et al., 2002), and also used to calculate a similar measure, the Instantaneous Rate of Green-up (Bischof et al., 2012). For arctic vegetation, a positive change in NDVI during spring corresponds with phenological periods of high nutrient concentrations and rapid vegetation growth (Finstad, 2008; Gustine et al., 2017).

We included five temporal windows (1 week before peak calving, the week of the peak, and the following 3 weeks after peak calving) for both NDVI metrics to assess at what temporal scale caribou may be responding to vegetation signals. To evaluate support for perception-based selection, we assigned the five temporal windows for both NDVI metrics relative to peak calving for that specific year, with the effect that the week of peak calving NDVI values differed between years and corresponded to the timing of calving observed the given year (perception of current conditions). To evaluate the potential for memory-based selection, we assigned these temporal windows relative to the herd’s average peak calving across all 8 years (June 3, see section “Results”), such that regardless of when peak calving was in a given year, both NDVI metrics represented consistent weeks across years (average conditions). This framework is similar to work assessing the influence of perception and memory in zebra (Equus burchelli) migration (Bracis and Mueller, 2017).

We tested the influence of these biotic and abiotic factors on caribou calving site selection using mixed-effects logistic regression, with use of a calving site as the response. We log-transformed elevation and VRM to approximate a normal distribution and standardized continuous covariates (mean centered and divided by the standard deviation) for model fitting. Correlation coefficients among physiographic attributes were under 0.5, and they all were under 0.2 when compared with environmental variables. We included a random intercept term for year to account for sampling across time and considered random slope terms for the environmental variables to account for stochastic annual variability (Gillies et al., 2006). We performed model selection at two stages – the first to select a random effect structure and the second to select fixed effect variables and structures (Bolker et al., 2009) using model selection based on Akaike’s Information Criterion corrected for small sample sizes (AICc; Hurvich and Tsai, 1989; Burnham and Anderson, 2002). For all NDVI, ΔNDVI, and snow-free variables, we fitted full fixed-effects models with a random slope term for each environmental covariate (including a random intercept for year) and compared performance with an intercept-only random effects model. In the second stage of model selection, we proceeded with fixed-effects selection using the top-performing random effect structure from the previous stage and included all biologically justifiable interactions and combinations. All analyses were performed in the R statistical program using the package “lme4” (Bates et al., 2015). We used our top model to generate a predictive map for calving sites by averaging the selected environmental covariate raster across the 8 years, as well as generated year-specific predictive maps with the corresponding environmental data for that year. We calculated the average of the year-specific predictive values within each annual calving area and compared these with the calving site values for the given year.

From 2010 to 2017, we detected 214 total calving events, ranging from 15 to 52 in a given year, and the average calving date was June 3 (Table 1). We identified calving in one non-migratory individual in 2016 and 14 in 2017 for which we used results from only the population-based model. We detected four calving events outside of the historical calving grounds: one in 2012 for an individual that remained on the winter range of the Seward Peninsula and three in 2017 for WAH individuals that calved in the Teshekpuk Herd calving area to the east. For the subsequent analyses of calving area trends and selection, we excluded these four events because they greatly skewed calving distribution estimates, leaving us with 210 total calving events across 8 years.

Across the 8 years we analyzed, the WAH calving areas exhibited variation at the annual scale, but the general area was characterized by remarkable fidelity. The average extent of the calving area for the herd in a given year was 24,772 km² (Table 1). Calving areas exhibited both latitudinal and longitudinal variation across years, with calving occurring in the Brooks Range and as far south as the Noatak River in some years (Figure 3). On an annual basis, calving areas had significantly less overlap than expected by chance three out of seven times (p < 0.05; Table 2). This trend of significant differentiation among years was evident at all further levels of comparison: at 2-year (p < 0.05 for five out of six), 3-year (p < 0.05 for four out of five), 4-year (p < 0.05 for three out of four), 5-year (p < 0.05 for one out of three), 6-year (p < 0.05 for one out of two), and 7-year intervals (p < 0.05 for the one comparison). When considered across years, the calving area of WAH females always shared a 7,281 km² core area of overlapping extent that was used every year of the study, with calving areas of less frequent use stretching as far away as the Noatak River (Figure 4) for a total calving ground extent of 53,330 km².

The selection of calving sites was characterized by mostly flat tundra within a band of elevation that was greening up at the time of average calving for the herd across all years of the study. The environmental covariate that explained the most variance in calving site selection was ΔNDVI at the average peak calving date for the study (“ΔNDVI.148”), from May 21–27 to May 28–June 3 every year (day of year 141–147 to 148–154; Supplementary Table 2) and substantially outperformed the next best model in random effect selection, which included a random slope for ΔNDVI at peak calving specific to each year (ΔNDVI.Calve; Δ AICc = 29.2). In model selection for fixed effects, the top performing model included terms for land cover, quadratic terms for elevation and solar radiation that indicate selection for intermediate values for both, an interaction between elevation and ΔNDVI.148, and terrain ruggedness (Supplementary Table 3). Females strongly selected for sites with high ΔNDVI at the time of peak calving (Table 3).

Calving sites were associated with a band of low elevation areas, indicating selection of elevations between approximately 50–600 m above sea level. Elevation and ΔNDVI exhibited an interactive effect, with females most strongly selecting for elevations approximately 100–175 m that were experiencing the fastest green-up at the time of peak calving. Of the four land cover classes we considered, we did not detect calving in any forested sites and found the strongest selection for herbaceous cover at calving (Table 3). The solar radiation index also exhibited a quadratic relationship for calving site selection (Table 3), indicating selection of sites with a positive index ranging from approximately 0.15 to 0.5, which correspond to lower angle slopes and encompass nearly all aspects. The negative linear coefficient for terrain ruggedness supported this result, indicating that females selected for less rugged terrain (Table 3). Our predictive map of calving habitat indicates that calving for the WAH occurs in the largest, continuous expanse of habitat characterized by these unique factors within their range, and that the attributes associated with calving sites extend to the east beyond documented calving areas (Figure 5). Importantly, calving occurred on sites with higher predicted value from the top model compared to the average of the calving area in the given year (Supplementary Figure 2).

Migratory ungulates rely on large expanses of range to maximize fitness (Hebblewhite et al., 2008; Joly et al., 2019) and migration routes of animals that rely on spatial memory are more susceptible to disturbance as they are likely more inflexible (Bracis and Mueller, 2017). Once lost, migratory patterns can take many generations for a population to learn and re-establish (Jesmer et al., 2018). Previous studies recommend that to fully conserve calving grounds for species such as caribou, managers should consider the full extent of calving at a decadal scale as the goal (Carroll et al., 2005; Taillon et al., 2012). Across 8 years of study, the WAH used an approximately 7,000 km² core area along the Utukok River for calving and a broader area of 53,330 km² to respond to environmental variability experienced each year on the calving ground. Comparing our findings with previous studies of the WAH up to six decades prior highlights the remarkable fidelity of this herd to its general calving ground (Supplementary Figure 3; Lent, 1966; Kelleyhouse, 2001) and local indigenous knowledge suggests this pattern extends before the 20th century (Lent, 1966; Burch, 2012). We recommend managers adopt the extent of the calving ground as the management goal for migratory caribou herds such as the WAH to ensure adequate space to respond to the annual environmental variability faced by caribou populations. We expect this recommendation has immediate utility for WAH management, for the area where the majority of calving occurs is on the National Petroleum Reserve – Alaska and specifically within the Utukok River Special Area. The Bureau of Land Management is currently revising the Integrated Activity Plan, which will designate conservation areas within the Reserve and stipulations on development in the greater area, and decisions made now have the potential to impact the WAH calving grounds for decades.