Irruptive dynamics of introduced caribou on Adak Island, Alaska: an evaluation of Riney‐Caughley model predictions

Abstract

Large mammalian herbivores introduced to islands without predators are predicted to undergo irruptive population and spatial dynamics, but only a few well‐documented case studies support this paradigm. We used the Riney‐Caughley model as a framework to test predictions of irruptive population growth and spatial expansion of caribou (Rangifer tarandus granti) introduced to Adak Island in the Aleutian archipelago of Alaska in 1958 and 1959. We utilized a time series of spatially explicit counts conducted on this population intermittently over a 54‐year period. Population size increased from 23 released animals to approximately 2900 animals in 2012. Population dynamics were characterized by two distinct periods of irruptive growth separated by a long time period of relative stability, and the catalyst for the initial irruption was more likely related to annual variation in hunting pressure than weather conditions. An unexpected pattern resembling logistic population growth occurred between the peak of the second irruption in 2005 and the next survey conducted seven years later in 2012. Model simulations indicated that an increase in reported harvest alone could not explain the deceleration in population growth, yet high levels of unreported harvest combined with increasing density‐dependent feedbacks on fecundity and survival were the most plausible explanation for the observed population trend. No studies of introduced island Rangifer have measured a time series of spatial use to the extent described in this study. Spatial use patterns during the post‐calving season strongly supported Riney‐Caughley model predictions, whereby high‐density core areas expanded outwardly as population size increased. During the calving season, caribou displayed marked site fidelity across the full range of population densities despite availability of other suitable habitats for calving. Finally, dispersal and reproduction on neighboring Kagalaska Island represented a new dispersal front for irruptive dynamics and a new challenge for resource managers. The future demography of caribou on both islands is far from certain, yet sustained and significant hunting pressure should be a vital management tool.

Introduction

Irruptive population dynamics often result when large mammalian herbivores are introduced to spatially confined habitats that lack top‐down control by predation or hunting (Leopold 1943, Gross et al. 2010), and these conditions often occur on islands. Riney (1964) and Caughley (1970) developed a four‐stage model of population growth and expansion of ungulates introduced to islands that predicts an irruptive rate of increase following release (stage 1), with population size eventually surpassing carrying capacity and growth rate slowing as time‐lagged density‐dependent feedbacks strengthen (stage 2). A steep decline in population size ensues (stage 3), and the subsequent population re‐equilibrates to a lower carrying capacity (stage 4). If islands are large enough, the model predicts a reset to stage 1 as animals move to new areas across a dispersal front or as barriers to dispersal are overcome. Importantly, the Riney‐Caughley model incorporates both a numerical and a spatial response, yet the shape and timing of irruptive dynamics can be influenced by different density‐dependent feedbacks with carrying capacity, species life history traits, and environmental conditions (McCullough 1997, Saether 1997, Clutton‐Brock and Coulson 2002, Forsyth and Caley 2006, Mysterud 2006, Gross et al. 2010, Tyler 2010). Moreover, the inevitability and consequences of irruptive behavior have been questioned (Caughley 1970, McCullough 1997, Gunn et al. 2003), which stems in large part from a paucity of time series data for populations that typically inhabit remote islands that are difficult to survey repeatedly.

Introductions of caribou and reindeer (Rangifer tarandus) on predator‐free arctic and sub‐arctic islands serve as useful natural experiments for examining complex interactions between intrinsic processes and extrinsic factors influencing ungulate population dynamics. These introductions were typically motivated by increased opportunity for sport and subsistence hunting or reindeer herding, but harvest levels often proved insufficient to limit population growth after humans abandoned remote islands or exerted nominal harvest pressure long‐term (Swanson and Barker 1992). At least 31 Rangifer introductions to islands have occurred world‐wide, with most displaying some degree of irruptive behavior numerically similar to the Riney‐Caughley model (Leader‐Williams 1988, Swanson and Barker 1992, Tyler 2010). However, mechanisms driving irruptions can vary widely among populations and habitats owing to four primary reasons. First, because islands typically lack a history of herbivory, herbivores released to islands encounter an accumulation of forage biomass that fuels rapid population growth, and overexploitation of this “surplus forage” in winter can be a major catalyst driving population crashes (Leopold 1943, Mysterud 2006, Gross et al. 2010). However, functional traits and herbivory tolerance of island vegetation can alter the magnitude of the irruption and subsequent stage 3 crash and stage 4 recovery. Although fat and somatic reserves acquired by ungulates during summer influences winter survival (Bardsen et al. 2010, Monteith et al. 2013), availability of winter forage largely sets carrying capacity for Rangifer and other northern ungulates (Klein 1991). Hence, populations introduced to islands with a winter forage base comprised of herbivory‐tolerant graminoids should be less prone to dramatic density‐dependent reductions in fecundity and calf survival, whereas populations subsisting on herbivory‐intolerant lichens, which often are abundant initially but then become rapidly depleted, are more prone to winter starvation, and are ultimately less persistent (Klein 1968, Leader‐Williams 1988, Mysterud 2006, Hansen et al. 2007, Gross et al. 2010). Some of the best known studies involve the latter case, where reindeer irrupted and then crashed to extinction on the Bering Sea islands of St. Paul, St. George (Hanna 1922, Scheffer 1951), and St. Matthew (Klein 1968, 1987, Klein and Shulski 2011) subsequent to the depletion of preferred winter forage‐lichens (e.g., Cladonia sp.) at high grazing densities. However, density‐independent weather events such as heavy snow‐fall or freezing rain over snow may also have contributed to the demise of these herds by creating “locked pastures” that rendered winter forage inaccessible (Klein 1968, Miller et al. 2005). Second, the amount of time that a population has persisted can influence the magnitude of irruptive growth and how that population responds to different sources of mortality (McCullough 1997). In a review of population dynamics of Rangifer on islands, including both introduced and long‐established natural populations, Tyler (2010) concluded that introduced populations display a consistent pattern of decline following an initial irruptive phase, whereas long‐established populations undergo frequent and unstable oscillations in population size. In both cases, however, density‐independent climatic conditions interacted with density‐dependent plant‐herbivore feedbacks to produce observed population dynamics, and evidence for locked‐pastures driving declines was largely anecdotal. Third, most introduced Rangifer are only hunted for a short time before human settlements are abandoned (Leader‐Williams 1988), yet populations subjected to multiple and discontinuous pulses of harvest intensity may exhibit growth patterns that diverge from the typical four‐stage shape. While long‐established ungulate populations released from hunting mortality can undergo weaker irruptions compared to newly introduced populations (McCullough 1997, Forsyth and Caley 2006, White et al. 2007, Tyler 2010), empirical studies of populations subjected to multiple periods of harvest are sparse. Fourth, the potential for prolonged delays in initial population growth and spatial expansion (i.e., latency or lag times) because of life history traits of the introduced species or novel environmental conditions is a fundamental concept of invasive species ecology (Crooks et al. 1999) that may influence the timing of irruptive dynamics (Forsyth and Caley 2006). Rangifer rarely produce twins and yearling females have low fecundity, so the species' reproductive potential is low relative to other cervids.

The Riney‐Caughley model is also spatially implicit. It predicts a dynamic spatial use pattern similar to a “rolling wave” (Caughley 1970) where herbivore density should be highest at the leading edge of the dispersal front and lowest at the initial release site where forage resources become depleted. Thus, if population growth and spatial expansion are not restricted by lag times (Crooks et al. 1999), island Rangifer should distribute themselves along a gradient of forage abundance in an ideal free fashion (e.g., Fretwell and Lucas 1970) where animal density matches per‐capita resource availability. This should then result in an isometric relationship between population size and density across areas with different probabilities of use (e.g., population core areas and home ranges), which provides further evidence of the Riney‐Caughley rolling wave of spatial expansion. However, Rangifer are highly vagile (Skoog 1968, Bergerud 2000), so introduced animals are capable of rapid and random expansion across unoccupied island habitats provided that barriers to intra‐island dispersal are lacking.

Ungulate spatial use patterns are also influenced strongly by sex‐ and age‐specific life history demands (McCullough 1999, Bowyer 2004) that could influence the rate of dispersal or spatial positioning of high density areas. In particular, mainland Rangifer can display strong inter‐annual site fidelity to calving areas (Gunn and Miller 1986) across a range of population densities. Site fidelity to calving areas has been proposed as an adaptation for parturition and calf rearing to coincide with low predation pressure (Bergerud 1996, 2000) or improving quantity and quality of forage during spring when energetic demands associated with lactation are high (White 1983, Post et al. 2008). On islands, however, whether or not Rangifer display site fidelity to calving areas becomes especially important during the irruptive periods of stages 1 and 2, which in turn could alter spatial predictions of the Riney‐Caughley model. Island habitat‐use should be unaffected by predation risk, so animals during the calving season should move freely to take advantage of the best available forage as a means of minimizing scramble competition, particularly as population size increases. This would lead to little spatial distinction between calving and post‐calving areas. In contrast, calving‐site fidelity would be reflected by calving areas that are spatially distinct from post‐calving areas. The center of the high‐density portion of the wave would also appear stationary across the range of population densities when spatial distributions during calving are examined on their own, even though the dispersal front is actually moving forward outside the calving season and the absolute size of calving areas increases with increasing abundance. This pattern of long‐term fidelity to calving areas across a range of population densities (from post‐release to the irruptive peak) could indicate that evolved fidelity to either better synchronize parturition with high resource availability or minimize predation risk (Bergerud 2000) does not breakdown when animals are placed in unfamiliar environments.

The introduction of caribou to Adak Island in the central Aleutian archipelago of Alaska set the stage for a type of natural experiment different from existing case studies that evaluated factors influencing population and spatial dynamics of island ungulates, and presents a unique case to test numerical and spatial predictions of the Riney‐Caughley model. Adak is remote and devoid of native land mammals, but it also supported a major U.S. military installation during World War II though the end of the Cold War. Hence, a large portion of island does not necessarily qualify as pristine wilderness. During 1958–1959, 23 barren‐ground caribou (R. tarandus granti) from the Nelchina herd on mainland south‐central Alaska were transplanted to the west side of Adak to provide recreational hunting and an emergency food supply for military personnel (Jones 1966). Adak caribou are similar to other insular Rangifer populations in that they are spatially restricted and have no predators other than humans, but differ in many key aspects. First, most studies (82%) of introduced Rangifer have occurred on islands above 60° latitude characterized by arctic‐tundra habitats and climates influenced by winter sea ice (see Tyler 2010). Less is known about populations introduced to lower‐latitude islands, such as the Aleutians (ca. 51° latitude), where mesic tundra habitats and maritime climates favor greater availability of graminoids that can be utilized as winter forage once lichens are overgrazed. Second, unlike most other island Rangifer, Adak caribou have been managed by sustained sport and subsistence hunting, and variation in harvest intensity could influence the timing and strength of irruptive dynamics (Forsyth and Caley 2006, White et al. 2007). In particular, we are aware of no studies of introduced Rangifer that have examined if an increase in harvest mortality can help decelerate population growth after the onset of a strong irruptive phase. Third, the Riney‐Caughley model describes the dispersal front as spreading into to new habitats within a single large island, whereas barriers to dispersal (i.e., water) typically restrict between‐island movements. However, Ricca et al. (2012) documented the very early stages of caribou invasion of Kagalaska, an island directly adjacent to Adak, but did not confirm the presence of reproductive females or their young. Evidence of an evolving dispersal front on Kagalaska that includes a reproductive component would represent a new twist in the manifestation of the Riney‐Caughley model. Last, the interaction between spatial use and irruptive population dynamics is an integral part of the Riney‐Caughley model, yet most studies of island Rangifer focus on the numerical response. Survey data available for Adak caribou allow examination of relations between population growth and spatial use.

We used a series of spatially‐explicit population counts spanning the 54‐year history of the Adak caribou herd to test numerical and spatial predictions based on the Riney‐Caughley model described heretofore. Numerically, we used population counts to test if population growth was characterized by irruptive dynamics, if post‐release growth rate was influenced by a latency period, and how periods of population stability and growth were influenced by harvest intensity and density‐independent weather conditions indexed by total snow and rainfall. We also conducted a simulation exercise to test if increases in harvest intensity following a period of relaxed harvest combined with strengthening density‐ dependent or independent mortality could explain the pattern of population growth during the last eight year of time series when a deceleration was detected. Spatially, we tested if high‐density areas moved outward from the point of release as population size increased (i.e., the rolling wave). We also measured the rate of expansion of core areas and home ranges (i.e., probabilities of use) at the population level to test if caribou densities within these probabilities of use increased isometrically with population size. On a seasonal basis, we tested if calving areas were spatially distinct from post‐calving areas. If spatially distinct, we then tested whether the location of calving areas remained stable across the entire range of caribou densities (i.e., post‐release period to irruption peak).

Methods

Study area

Adak is a large (725 km²) mountainous island located in the Andreanof group of the central Aleutian archipelago (51.872 N, 176.636 W). Kagalaska (116 km²) is the closest island directly to the east and is separated from Adak by a narrow (∼1 km) strait (Fig. 1). Adak has a mixed land‐ownership pattern that adds complexity to caribou management. The northern end comprises a former naval base now held in part by a private Alaska Native Corporation, and the remainder of the island is largely wilderness lying within the U.S. Fish and Wildlife Service (USFWS) Alaska Maritime National Wildlife Refuge (AMNWR). All of Kagalaska is managed by the AMNWR and has no human habitation. Human population size for the city of Adak dropped from >6000 in the early 1990s to <330 in 2010. Predation by Bald Eagles (Haliaeetus leucocephalus) on neonatal calves is a documented but insignificant source of mortality (Jones 1966).

Topographic map of Adak and Kagalaska Islands, located within the central Aleutian archipelago of Alaska. Area to the north of the space use boundary (North Adak) was excluded from analyses of spatial use, and the star indicates the initial release site for 23 barren‐ground caribou calves in 1958–1959 (Jones 1966).

The maritime climate and plant community of Adak and Kagalaska is similar to that of the entire Aleutian archipelago. Summers are characterized by cool (5–10°C) temperatures and persistent fog and rain, while winter temperatures fluctuate near 0°C and cyclonic storms occur frequently. Plant communities are compositionally similar to the maritime tundra described by Talbot et al. (2010) and can be broadly characterized by graminoid meadows (e.g., Calamagrostis nutkaensis, Carex sp., Erigeron peregrinus, Anemone narcissiflora), and heaths with evergreen (e.g., Empetrum nigrum, Phyllodoce aleutica, Loiseleuria procumbens) and deciduous (Vaccinium uliginosum, Salix sp.) dwarf‐shrubs. Caribou have significantly depleted the abundance of Cladonia lichens on Adak, whereas robust stands are ubiquitous on Kagalaska (Ricca 2013). The plant growing season is relatively short and occurs primarily from June to August.

Caribou surveys

Ground‐based counts were conducted near the original release site on west Adak where the population remained through 1964 (Jones 1966). Aerial surveys were conducted primarily from available military aircraft (fixed wing or helicopter) by AMNWR and other biologists at periodic intervals from 1967 to 1989 (n = 24) and from a Bell‐Jet Ranger helicopter from 1993 to 2012 (n = 4). Parallel transects 1.0–1.5 km wide were flown typically in areas with flat topography, whereas landscape contours were followed in mountainous areas. Specific flight speeds and paths varied among surveys according to aircraft used and weather conditions. Two or three observers counted caribou within groups and recorded group locations on 1:100,000 topographic maps. Animals were classified according to age class (adult or calf) during most surveys conducted during the calving season. During 2012, group locations were recorded by GPS, and ground‐based calf‐adult counts were conducted in observed calving areas (Teardrop Basin, Fig. 1) immediately after completing the aerial survey. Ratios of calves to adults obtained from aerial and ground‐based counts did not differ significantly (ANOVA: F = 0.2, P = 0.89), so we consider our aerial estimates accurate. The northern section of Adak was not surveyed consistently across the time series. However, few animals likely inhabited the north side of the island during early spring and summer due to its substantial road network and concentration of human activity. Less than 3% of the total population was detected on the north side during 2005 and 2012 when the entire island was surveyed. Kagalaska was surveyed along with Adak in 2012 following ground‐based detections of caribou sign and a single caribou group on Kagalaska during the prior two summers (Ricca et al. 2012). Surveys occurring before and after 1994 were categorized as military and post‐military, respectively, to help partition changes associated with the closure of the naval base and dramatic reduction in human population size.

Population dynamics

All population counts reflected summer population size after calving occurred within a given year. The few counts obtained during winter and early spring before calving (10% of the time series, all during the military period) were standardized to estimate population size during the preceding summer counts by adding the number of animals harvested to the observed count. Hunting season did not begin until the fall during the military period, so harvested animals were added to the observed count assuming they would have been alive and counted during the preceding summer (USFWS, unpublished manuscripts). Collection of annual harvest data was managed intensively by AMNWR during the military period (USFWS, unpublished manuscripts) when a permit system was strictly enforced and civilian access to the naval base was highly restricted. Post‐military harvest data based on harvest‐ticket returns were obtained from the Alaska Department of Fish Game's harvest website (ADFG 2012). Mail‐in harvest reporting, the lack of on‐site license vendors, and increased subsistence hunting during the post‐military period may have led to under‐reported harvest levels compared to levels during the military period.

We used multiple linear regression to test effects of harvest size and precipitation (snowfall and rainfall) on population growth rate (r) for the time series prior to the naval base closure. Explanatory variables in the regression model were log transformed and lagged one year to better synchronize effects on reproduction and survival influencing that year's count. Annual rainfall and snowfall data for Adak (Station Number 500026) from 1958 to 1993 were downloaded from the National Climate Data Center (NCDC; http://www.ncdc.noaa.gov/cdo‐web). Continuous weather monitoring on Adak ceased after the naval base closed. Precipitation measures from the closest weather stations, >300 km away at St. Paul and Dutch Harbor, correlated poorly (R² ≤ 0.46) with those from Adak. Therefore, we could not extend the analysis across the entire time series.

We evaluated the relationship between log population size and year using segmented regression (Muggeo 2003) to test for period‐related growth rates (e.g., pre‐ and post‐military) characterized by significantly abrupt changes in slope across the entire time series. Reported harvest levels from both periods were included as a covariate. We compared linear, two‐, and three‐slope models using Akaike Information Criterion (AIC_c; Burnham and Anderson 2002) and evaluated the significance of each slope (r value) with 95% confidence intervals. Regressions were run using the segmented package (Muggeo 2008) in Program R (R Development Core Team 2012). Lag times were evaluated by comparing measured population growth during the initial years of the time series with the theoretical maximum rate for Rangifer (r_max = 0.30–0.31; Leader‐Williams 1988, Heard 1990).

Using known attributes of caribou life history and realistic mechanisms associated with irruptive population dynamics, we ran a set of eight simple‐simulation models (e.g., White 2000) to help explain the population trajectory between the 2005 and 2012 surveys (see Results). Single factor models consisted of: (1) the known increase in reported harvest (H) after 2005; (2) high levels of unreported harvest leading to higher total harvest (high H); (3) intensifying density‐dependence (DD), and (4) episodic density‐independent mortality (DID). Multi‐factor models included: (1) H + DD; (2) high H + DD; (3) H + DD + DID; and (4) high H + DD + DID. We did not evaluate two‐factor models only consisting of harvest and density‐independence because density‐independent limiting factors often interact with density‐dependent regulating factors (Saether 1997, Sinclair 2003, Tyler 2010). Plausible model outcomes would yield ranges overlapping the observed 2012 population size and calf:adult ratio. Vital rates and herd composition ratios were unknown, so we used published values from other growing Rangifer herds in predator‐free, sub‐arctic habitats as input parameter estimates for our simulations. For all years, we set the percentage of breeding adult (≥2 year) females at 37.5 % (±7.5%) and adult males at 30% (±5%) (Skoog 1968, Leader‐Williams 1988). Vital rates for the first two years were set as follows for all models: pregnancy rates = 90% (±5%), neonatal (0–3 month) calf survival = 97.5% (±2.5%), and adult and overwinter calf survival = 95% (± 5%) (Skoog 1968, Skogland 1986, Leader‐Williams 1988, Bergerud 2000). Mortality was 1 − survival. Beginning in spring 2005, we set N_{pre‐calving(t)} at 2100 caribou (± 200). We calculated N_calves(t) as N_{pre‐calving breeding females(t)} × pregnancy rate_(t) × neonatal calf mortality_(t), and N_recruits(t) as N_calves()(t) − overwinter calf mortality_(t). Total mortality_(t) equaled the sum of adult mortality_(t), total harvest_(t), and calf mortality_(t). Net recruitment equaled N_recruits(t) − total adult and harvest mortality_(t). Pre‐calving N_(t+1) equaled N_{post‐calving(t)} + net recruitment_(t). We incorporated demographic stochasticity to each unmeasured parameter by taking 1000 draws from a uniform distribution spanning the range of input values, and ran each simulation 1000 times. We nominally increased H by 15% (range: ±5) to account for the lack of hunting license vendors on Adak that likely contributed to under‐reported harvest (Adak residents, personal communication), and to provide a more realistic estimate since harvest reporting rates are rarely 100% (Kilpatrick et al. 2005). High H values were calculated by adding a 55% (±5%) increase to reported harvest to account for potentially high levels of unreported local subsistence take (Adak residents, personal communication). Because additive harvest mortality by definition cannot occur if density‐dependence is operative in a population (Sinclair and Pech 1996), we treated harvest mortality as ‘semi‐additive' in all eight models by estimating the number of adults dying from natural death in winter independently from those harvested. Moreover, density‐dependence in Rangifer (Adams and Dale 1998, Bergerud 2000) and most cervids (Gaillard et al. 2000) manifests primarily through reductions in fecundity and recruitment, so we allowed vital rates to remain unchanged for two years and then reduced pregnancy rates annually by 3.3% (±10%), neonatal calf survival by 2% (±5%), and overwinter calf survival by 5% (±5%) from 2007 to 2012. In contrast, adult caribou survival is dramatically less sensitive to density‐dependent reductions, yet significant over‐winter mortality can occur when deep snow or freezing rain over snow precludes access to forage (Scheffer 1951, Klein 1968, Gates et al. 1986, Gunn et al. 2003). Although we could not model winter weather effects due to the lack of weather data for Adak after the base closure, local knowledge indicated that the winters of 2006 and 2011 had higher than normal and more persistent snow (L. Spitler, USFWS, personal communication). Significant reindeer winter kill on St. Paul Island (approximately 700 km NE of Adak) was also reported by local media in the winter of 2011–2012. Thus, we conservatively modeled density‐independent mortality by lowering adult (and calf) survival rates in 2006 to 15–20% and 20–25% in 2012. These reductions are comparable to those experienced by other predator‐free herds (Reimers 1982, Leader‐Williams 1988, Reimers 2012).

We also evaluated if the population trajectory could be explained by more compensatory harvest mortality in two different sets of model runs. First, we re‐ran all eight models with the identical suite of parameter inputs except for adult mortality, which was re‐calculated by subtracting the total harvest from the number of adults at risk of over‐winter mortality. Hence, fewer adults were at risk of dying in winter as harvest increased. Second, density‐dependence should weaken as harvest increases, particularly if exact compensation in recruitment occurs whereby equilibrium population size remains unchanged as harvest increases (Sinclair and Pech 1996). Accordingly, we re‐ran the 4 models that combined harvest and density‐dependence (i.e., H + DD; high H + DD; H + DD + DID; high H + DD + DID), and only kept density‐dependent vital rates at the original set levels for the first year (2005) to match the known population size. We then reversed the strength of density‐dependence applied over time for the original (and semi‐additive) model runs, which simulated an immediate and strong reduction in vital rates associated with a K‐overshoot followed by an increase in recruitment associated with lagged harvest effects. Specifically, pregnancy rates, neonatal calf survival, and over‐winter calf survival were lowered to 65% (+10%), 20% (±5%), and 70% (±5%) for years two and three, respectively, and then allowed increase incrementally over the remaining years to the original starting levels. Compensatory adult mortality was modelled in the same manner described above.

Spatial use estimation

Group locations and sizes recorded on topographic maps during aerial surveys were digitized into ArcGIS version 10.0 (Environmental Systems Research Institute, Redlands, California, USA). We used fixed‐kernel probability density functions to calculate utilization distributions (UDs) that approximated different intensities of caribou use. We calculated fixed kernel UDs using the kde function in Geospatial Modeling Environment (Beyer 2012), with likelihood cross‐validation (CVh) smoothing parameters calculated by the ‘ks' package (Duong 2012) in Program R (R Development Core Team 2012). Likelihood based cross‐validation performs better than least‐squares methods (e.g., LSCV) for data sets containing <50 observations (Horne and Garton 2006), which was a common feature of all of our military‐period data sets. Locations were weighted by group size so that larger groups contributed more local density to the kernel estimated UD. We calculated 95 and 99% UDs as estimates of total range size, and 50% UDs as estimates of intensive ‘core area' use. We only used surveys occurring from late spring to early fall and estimated UDs for surveys containing >17 groups, which resulted in a total of 13 annual surveys suitable for spatial analyses. Calving for the parent Nelchina herd (Skoog 1968) and most other Rangifer herds (Bergerud 2000) peaks in late May, so surveys conducted from late May through early June (1967, 1984, 1993, 2012) were classified as ‘calving', and surveys conducted from late July to early October (1977, 1978, 1982, 1983, 1985, 1988, 1989, 1998, 2005) were classified as ‘post‐calving'. Areas of high post‐partum caribou density often delineate calving grounds for management purposes (Taillon et al. 2012); hence, we considered 50% UDs from the calving season to be representative of calving areas. Estimating separate UDs for groups with and without calves would have resulted in insufficient sample size and inflated smoothing for at least two of the four calving season surveys. Hence, calving season UDs were estimating using all group locations and groups with and without calves were plotted with UD estimates for visual representation. We excluded the infrequently‐surveyed north side of Adak (Fig. 1) from all spatial use analyses to maintain comparability across the time series. Kagalaska was not surveyed systematically by air in conjunction with Adak and was excluded from spatial use analyses prior to 2012. Notably, caribou were not detected in an aerial survey of Kagalaska in 2003 (J. C. Williams, unpublished manuscript).

We evaluated spatial similarity and temporal stability within and among core areas estimated during calving and post‐calving with two approaches. The UD approach quantified spatial overlap and shifts in spatial distribution over time in a probabilistic manner, but was prone to sample‐size dependent variation in the estimation of smoothing parameters that can lead to over‐ or under‐estimation of UD size (Seaman et al. 1999). The point‐based approach tested for differences among distributions of caribou groups and centroids, which was less prone to sample size related variance but lacked estimates of areal extent.

For the UD approach, we calculated spatial overlap between all pairwise comparisons for 50 and 95% calving and post‐calving UDs as the average of UD_1,2 (areal proportion of UD₁ overlapped by UD₂) and UD_2,1 (areal proportion of UD₂ overlapped by UD₁) (Kernohan et al. 2001). We constructed linear contrasts within analyses of variance to test: (1) if intra‐seasonal spatial overlap within calving and post‐ calving seasons was greater than inter‐seasonal overlap; and (2) if overlap within just calving seasons was greater than inter‐seasonal overlap. We then spatially merged all 50% UDs by season to create an aggregate calving and post‐calving core area representative of the entire time series, and calculated the aggregate overlap between calving and post‐calving areas. We used a modification of Gunn and Miller's (1986) criteria for defining stable calving grounds, where <50% overlap in calving season core areas among years represented a lack of calving ground site fidelity.

For the point approach, we calculated the median centroid weighted by group size for each calving and post‐calving distribution, measured centroid distance from release site as a measure of dispersal, and calculated all pairwise centroid distance combinations within and among seasons. We constructed the same set of linear contrasts used for the UD overlap analysis to ask: (1) if all intra‐seasonal distances were less than inter‐seasonal distances; and (2) if distances within calving seasons only were less than inter‐seasonal distances. Multiple Response Permutation Procedure (MRPP) (Kernohan et al. 2001) was used to test if calving and post‐calving caribou group locations came from different probability distributions. Significant differences were evaluated by the ‘A' statistic, where A > 0 indicates greater within group homogeneity than expected by chance, along with the standard P value (McCune and Grace 2002). We also calculated standard distance (a measure of point dispersion around a geometric mean center) weighted by group size as measure of spatial variation. Centroids and standard distances were calculated using the spatial statistics toolbox in ArcGis 10, and t‐tests were used to evaluate differences between seasons.

We quantified the rate of range expansion by calculating the cumulative area added to aggregate 50 and 95% UDs across sequential surveys divided by the number of years elapsed between surveys. We then calculated caribou density within 50, 95, and 99% UDs and regressed log_e caribou density against log_e population size. A hypoallometric (slope < 1) relation would indicate fast spreading and widely dispersed animals at low population sizes, a hyperallometric (slope > 1) relation would indicate that animals remained more tightly clumped in particular areas as population size increased, and a isometric relation (slope = 1) would indicate that caribou congregated at an even rate with population size.

Results

Population dynamics

The introduced population (N = 23) grew to approximately 200 animals by the late 1960s, and then fluctuated between 200 and 600 animals through the remainder of the military period (Fig. 2). Population size irrupted to 2750 animals roughly 10 years after the closure of the naval base, but then appeared to level off at nearly 2900 animals in 2012. Hunting began in 1964 with >20% of the herd harvested annually by the early 1970s (Fig. 2). Annual harvest averaged 32% (range: 19–48%) with females accounting for 48% (range: 32–62%) of the total through the remainder of the military period. Hunting was allowed in all years with the exception of the 1983–1984 hunting season, which was closed following two successive years of population decline. Reported harvest plummeted to approximately 3% (1–11%) annually during the first half of the post‐military period before increasing during the latter half (2005–2011) to numerical levels similar to those reported during the military period. However, the harvested percentage of the herd remained low at 9% (8–11%) owing to the much larger population size. Females comprised 22% (4–64%) of the total reported harvest from 1994 to 2006 before increasing to 58% (43–67%) after new regulations limiting male harvest and allowing unlimited female harvest were enacted in 2007.

Time series of population size (grey circles) and reported hunter harvest (open circles) for Adak caribou, 1958–2012.

The rate of population growth during the military period was more strongly related to harvest than weather (F_{3, 18} = 7.4, P = 0.002). Harvest in the previous year correlated with lower population growth rate (t = −3.3, P = 0.004), but total rainfall or snowfall in the previous year had no effect (t ≤ 0.88, P ≥ 0.39). Across the entire time series, the relationship between harvest‐adjusted population size and year was best described by a three‐slope model that was 14 and 38 ΔAIC_c units lower than the two‐ and linear‐slope models, respectively (Fig. 3). Extremely rapid population growth approaching the theoretical maximum for the species occurred from 1958 to 1969 (r = 0.29, 95% CI = 0.11–0.47). Growth then slowed markedly from 1970 to 1985 (r = 0.01, 95% CI = −0.02–0.04) before undergoing a second irruptive period from 1986 to 2012 characterized by a slower overall growth rate across the 24‐year time span (r = 0.09, 95% CI = 0.07–0.10) but a steeper eight‐fold increase on a numerical basis.

Segmented regression of Adak caribou population size (partial residuals accounting for one‐year lag harvest effects) against year illustrating three significantly different slopes from 1958 to 2012. R values (slopes) are indicated next to each significant segment.

Our first set of simulated models indicated that semi‐additive harvest alone likely did not explain the similarity between the 2005 and 2012 population estimates (Fig. 4). Harvest (H) and high harvest due to unreported harvest (high H) models were not realistic since they yielded continuous irruptive population growth along with calf:adult ratios (0.47, range: 0.33–0.62) that clearly exceeded the 2012 observed value (0.28). All four density‐dependent and density‐independent models combined with either level of harvest yielded calf:adult ratios overlapping the 2012 value. However, only models incorporating high harvest and density‐dependence (high H + DD) or harvest, density‐dependence, and density‐independence (H + DD + DID) yielded a range of population estimates overlapping the 2012 count. High H + DD yielded decelerated growth through 2009 and declining growth thereafter, while H + DD + DID yielded a similar trend with a more dramatic decline from 2011 to 2012. In contrast, H + DD yielded logistic growth stabilizing at approximately 4000 animals, and high H + DD + DID yielded a relatively stable population of approximately 3000 animals through 2009 but a precipitous decline thereafter. Excluding harvest from density‐dependent and density‐independent models resulted in continued irruptive growth (not shown in Fig. 4).

Time series of simulated post‐calving population size estimates (box‐plots) for Adak caribou, 2005–2012 modeled with reported (H) and high‐unreported (High H) adult harvest, intensifying density‐dependent fecundity and mortality (DD), and episodic density‐independent mortality (DID). Population size estimated from the 2012 survey is denoted by an asterisk; note that scales differ across y‐axes.

None of the models with harvest treated as compensatory in relation to adult over‐winter mortality yielded a range of population estimates that overlapped the 2012 count (Appendix: Table A1). While the high H + DD and H + DD + DID models best approximated the observed trend and yielded matching calf:adult ratios, their minimum simulated counts overestimated the 2012 count by 125 and 195 animals, respectively. Similarly, none of the models with harvest treated as compensatory in relation to recruitment and adult‐overwinter mortality yielded a range of population estimates that overlapped the 2012 count, although the maximum simulated count from the H + DD + DID model only underestimated the 2012 count by 76 animals (Appendix: Table A1). Irruptive growth resumed after 2010 with the H + DD model, while decreasing population trends that resulted from the high H + DD and high H + DD + DID models indicated modelled levels of compensatory recruitment and over‐winter adult mortality were insufficient to compensate for high rates of unreported harvest. Calf:adult ratios from all models in both sets (range: 0.34–0.63) overestimated the observed 2012 calf‐adult ratio (0.28), which provided more evidence against compensatory harvest.

Table A1. Simulated post‐calving population size estimates for Adak caribou from 2005 to 2012 calculated from two sets of models where harvest mortality was compensatory relative to: (1) adult overwinter mortality; and (2) recruitment and overwinter mortality.

Spatial use

Total ranges during the post‐calving season represented by 95 and 99% UDs expanded across the entire study area by the early 1980s, whereas high‐density core areas represented by 50% UDs shifted to the central and east side of the island during the height of the second irruption in 2005 (Fig. 5). Total ranges during the calving season also expanded across the study area during the 1980s, but core areas were located consistently on the south‐central portion of Adak across the range of population densities. Importantly, a single caribou group comprised of three adults and one calf (<3 weeks old) was observed on the southeast corner of Kagalaska during the 2012 survey (Fig. 6).

Distribution of groups detected during aerial surveys and fixed kernel 50% (core area), 95%, and 99% (population range) UD estimates for Adak caribou during the post‐calving period (late July–October), 1977–2005.

Distribution of groups detected during aerial surveys and fixed kernel 50% (core area), 95%, and 99% (population range) UD estimates for Adak caribou during the calving period (May–June), 1967–2012. Note that groups with and without calves were not categorized in 1984.

Utilization distribution and point‐based approaches indicated spatially distinct and temporally stable calving areas, and more widely dispersed and temporally variable post‐calving areas. Average pairwise overlap between individual core areas differed significantly among all seasonal comparisons (F_2,75 = 12.1, P < 0.0001). Intra‐seasonal overlap ([IMAGE OMMITED. SEE PDF.] = 35%) and overlap within only the calving season ([IMAGE OMMITED. SEE PDF.] = 43%) was greater than inter‐seasonal overlap ([IMAGE OMMITED. SEE PDF.] = 16%) (t ≤ −2.5, P ≤ 0.01). Calving season core areas remained stable in the south‐central side of Adak since confidence intervals for spatial overlap (29–57%) included 50%, whereas confidence intervals for overlap within the post‐calving season (27–39%) and between seasons (11–22%) were less than 50%. Furthermore, only 24% of the aggregate post‐ calving core area shared space with the aggregate calving core area (Fig. 7A). In contrast, calving and post‐calving 95% UDs were not spatially distinct since average pairwise overlap did not differ among any intra‐ or inter‐ seasonal comparisons ([IMAGE OMMITED. SEE PDF.]_calving = 71%,[IMAGE OMMITED. SEE PDF.]_{post‐calving} = 72%,[IMAGE OMMITED. SEE PDF.]_{calving‐post‐calving} = 67%; F_2,74 = 1.7, P = 0.19). Aggregate post‐calving and calving 95% UDs overlapped by >88%.

Overlap between aggregate post‐calving and calving season 50% UD range estimates (A) and median centroids (B) for post‐calving (pink to red) and calving (light blue to dark blue) point distributions for Adak caribou, 1967–2012.

All seasonal comparisons for centroid distances were marginally significant (F_2,75 = 2.5, P = 0.08), yet intra‐seasonal ([IMAGE OMMITED. SEE PDF.] = 6.1 km) and within calving only ([IMAGE OMMITED. SEE PDF.] = 3.7 km) distances between centroids were significantly less than inter‐seasonal distances ([IMAGE OMMITED. SEE PDF.] = 7.2 km) (t ≥ 2.2, P ≤ 0.03). Calving area centroids were also farther away ([IMAGE OMMITED. SEE PDF.] = 9.3 km) from the release site than post‐calving centroids ([IMAGE OMMITED. SEE PDF.] = 6.2 km) (Fig. 7B). MRPP indicated that calving group distributions were spatially distinct from post‐calving group distributions (P = 0.02, A = 0.16). Post‐calving groups were also more widely dispersed (standard distance = 9.2 km) than calving groups (standard distance = 6.8 km) (t₁₁ = −2.6, P = 0.04).

Among sequential surveys, the aggregate post‐calving core area expanded more rapidly ([IMAGE OMMITED. SEE PDF.] = 17.8 km²/yr) than the aggregate calving core area ([IMAGE OMMITED. SEE PDF.] = 1.5 km²/yr) (t₇ = −2.4, P = 0.02). However, the cumulative percent area added became asymptotic in the late 1980s, and then increased again during the height of the second irruption for both seasonal distributions (Figs. 5, 6, and 8). Similarly, the aggregate post‐calving range expanded more rapidly (30.3 km²/yr) than the aggregate calving population range (4.1 km²/yr) (t₇ = −1.9, P = 0.05). The cumulative percent area added to the total range became asymptotic near 100% by the early 1980s as most of the island became occupied (Figs. 5, 6, and 8). Significant (t₁₁ ≥ 7.5, P ≤ 0.0001) isometric relationships occurred between log population size and log density in 50, 95, and 99% UDs with 95% confidence intervals for all slopes overlapping 1.0 (50%: 0.8–1.4, 95%: 0.7–1.1, 99%: 0.8–1.1) (Fig. 9). All relationships remained significantly (t₇ ≥ 4.8, P ≤ 0.002) isometric when only post‐calving surveys were included (CIs: 50% UD = 0.5–1.5, 95% UD = 0.6–1.2, 99% UD = 0.6–1.2).

Cumulative area added to core areas (50% UD, circles) and range area (95% UD, inverted triangles) per year during calving (open symbols) and post‐calving (grey symbols) seasons for Adak caribou, 1967–2012.

Relationship between total population size and density within core areas (50% UD, open circles) and range area (95% UD, grey circles) for Adak caribou, 1967–2012.

Discussion

Population dynamics

Irruptive population dynamics occurred for Adak caribou throughout the 54‐year history of the herd, but the shape of the pattern did not completely fit the Riney‐Caughley model. We detected two distinct periods of irruptive population growth separated by a long period of relative slow growth, and this pattern was linked to variation in hunting pressure during the military period. The first irruption occurred immediately after release whereby the population grew at a rate near the maximum for the species until managed harvest exceeded approximately 20% of the population in the mid 1960s. Hence, numerical latency was essentially non‐existent post‐introduction. Although Rangifer are notably less fecund compared to other cervids, documented rapid growth (r = 0.23–0.26) and irruptive peaks occurring ≤20 years post‐release for several populations inhabiting sub‐arctic and high arctic islands (Klein 1968, Heard 1990, Swanson and Barker 1992) combined with our results indicate an overall lack of numerical latency for the species. Like other newly introduced island herds (Klein 1968, Leader‐Williams 1988), Adak caribou attained body sizes that were larger than those from mainland herds (Jones 1966, Klein 1968). Small changes in forage quality and abundance can enhance Rangifer body condition, which in turn yields disproportionately large increases in fecundity and population growth (i.e., mutiplier effect; White 1983, Cebrian et al. 2008). Similar to previous studies of Rangifer, access to ample and high quality forage resources likely enabled large female body size that drove successful reproduction year after year during the first period of irruptive growth on Adak when exceptionally high diet quality was reported (Klein 1968).

The hunting closure of 1983–1984 appeared to act as a strong catalyst for the second irruption that began in 1986. The doubling in population size shortly after the one‐year hunting closure may have enabled population growth to withstand subsequent increased harvest levels, and the irruption was further exacerbated by the dramatic decreases in harvest following the naval base closure. For established Rangifer populations on islands, harvests exceeding 15% often correlate with a decline in growth rates or exacerbate stage 3 crashes (Gates et al. 1986, Ouellet et al. 1994, Tyler 2010). In the case of Adak, high (ca. 30% total, 50% female) reported harvest rates were necessary to slow but could not stop population growth during the military period. Established ungulate populations released from harvest should also experience weaker irruptive growth relative to newly introduced populations because herbivory from established populations has already reduced surplus forage resources (Caughley 1979, McCullough 1997, Forsyth and Caley 2006). Adak caribou do not appear to fit this paradigm. Although the growth rate during the second irruption (r = 0.09) was much lower than the post‐release growth rate (r = 0.29), population size increased ∼800% after the onset of the second irruption and nearly 300% after the naval base closure and concomitant decline in hunting pressure. This increase was much higher than those displayed by other ungulates released from long‐term hunting pressure, including an approximate 150% increase for bighorn sheep (Festa‐Bianchet et al. 2003, as cited by Forsyth and Caley 2006) and a 60% increase for red deer (Clutton‐Brock et al. 1997, as cited by Forsyth and Caley 2006).

Perhaps most intriguing was the marked deceleration of population growth between 2005 and 2012 that gave the appearance of theta‐logistic growth for a population approaching carrying capacity. However, unmeasured processes such as delayed density‐dependence typical of many populations that display irruptive dynamics (Forsyth and Caley 2006) may have been working during the seven years between surveys, which led us to conduct our ‘what if' modeling exercise. Caution needs to be exercised when interpreting model output (White 2000) since input variables were comprised largely of demographic parameters measured for other herds and local knowledge of unreported harvest rates and weather severity. Still, our demographic parameters were based on average published values for introduced Rangifer herds inhabiting similar environments, and we feel confident in local knowledge that matches our years of experience on Adak. In addition, the limited empirical data on recruitment and survival did not allow a priori identification of whether harvest was additive or compensatory. Un‐modelled variation in vital rates in relation to harvest and carrying‐capacity could have also yielded different outcomes, yet results from two sets of simulations for compensatory harvest mortality were not supported since they did not overlap with the observed 2012 count and calf:adult ratios. Moreover, the mere presence of density‐dependence does not necessitate compensatory mortality, and exact compensation from harvest is difficult to prove without specific data on density‐dependent mortality and changes in population size in situations with and without hunting (Sinclair and Pech 1996). Thus, harvest likely acted in a more additive than compensatory manner on Adak where predators were absent and caribou experienced a milder maritime climate.

Our modeled results point to the importance of continued hunting pressure in limiting population growth of Adak caribou, but they also suggest that reported or unreported harvest alone was insufficient to halt the second period of irruptive growth without additional density‐ dependent or independent mechanisms. The observed fit of the semi‐additive harvest + density‐ dependence and independence (H + DD + DID) model is plausible given that the prevailing paradigm of large herbivore population dynamics recognizes the interaction of density‐ dependent and independent factors (Saether 1997, Tyler 2010). However, we failed to detect large numbers of winter‐kill carcasses during our 2012 aerial and ground surveys that followed a winter with heavy snowfall. Hard ice and deep snow conditions may stimulate irruptions and crashes on high arctic islands and polar deserts with low standing crops of winter forage biomass (Gates et al. 1986, Ouellet et al. 1994, Aanes et al. 2000, Gunn et al. 2003, Miller et al. 2005), but these conditions are less likely in the more benign climate of the Aleutians. Time‐lagged weather effects also did not strongly influence population growth during the military period, which lends less support for the H + DD + DID model. Thus, the increase in known plus unreported harvest combined with strengthening density‐dependence presents a more likely set of conditions that helped slow the second irruption. High levels of harvest can replace density‐dependent processes regulating population growth (White et al. 2007), but our high H model that culled 15% of the population resulted in continued growth, so additional density‐dependent mechanisms were needed to slow the irruption. Notably, our modeled density‐dependence was relatively moderate when compared to the extreme lower end of vital rates measured in other predator‐free Rangifer herds (e.g., calf survival = 50%, pregnancy rate = 7% (Bergerud 2000)). Observed calf:adult ratios in 2012 that were lower than those associated with rapidly growing populations (>0.40) but higher than those associated with rapidly declining herds (≤0.17) (Valkenburg et al. 1996) suggest that our modelled moderate‐level of density‐dependence was reasonable.

Habitat productivity can further influence the strength and timing of density‐dependence (Leader‐Williams 1988, Mysterud 2006). The concept of surplus forage is a key component of the Riney‐Caughley and other irruptive models (Gross et al. 2010), and herds that forage within productive, graminoid‐dominated communities such as those on Adak may be less prone to extreme irruptions and crashes compared to herds exploiting surplus biomasses of herbivory‐intolerant lichens in less productive habitats as has been demonstrated for island (Klein 1968, Leader‐Williams 1988) and mainland populations (Olofsson et al. 2004, van der Wal 2006). In addition, a trend towards a warmer and drier climate in the Bering Sea region may further facilitate graminoid production and inhibit lichen regeneration (Klein and Shulski 2011). While we stress strongly that the post‐irruption shape of the times series is still unfolding, it is plausible that Adak caribou underwent a delayed density‐dependent response mediated by hunting pressure between 2005 and 2012.

Spatial use

To the best of our knowledge, no studies of introduced island Rangifer have measured a time series of spatial use to the extent described in this study. The pattern of spatial expansion for Adak caribou, particularly core areas during the post‐calving season, aligns well with the predictions of the Riney‐Caughley model, and fidelity to distinct calving areas across the range of population densities was supported by UD and point‐based methods. The Riney‐Caughley model predicts that the highest densities will appear at the leading edge of the dispersal front, and post‐calving 50% UDs followed that pattern as population size increased. Core areas were initially located near the release site on the west and southwest side of the island, expanded and stayed towards the center of the island during the military period when population size fluctuated between roughly 200 and 600 animals, and then expanded into the east side of the island during the height of the second irruption. In contrast, the total population range expanded across the island area by the early 1980s. It then expanded again into Kagalaska during 2012 after the height of the second irruption, forming the basis for a new stage 1. Spatial expansion of caribou across Adak also appeared to follow an ideal free distribution (Fretwell and Lucas 1970) owing to the isometric relation between density (within all utilization distributions) and total population size, as well as the overall pattern of eastward core‐area expansion. In particular, the patterns for post‐calving core areas indicate that high numbers of caribou likely concentrated within spatially limited areas that met per‐capita resource needs without extensive exploration, and then dispersed into new areas to the east as population size and competition potential increased. Other studies have suggested that caribou distribution and range expansion is driven by density thresholds. For mainland herds, Skoog (1968) predicted that caribou spatial distributions would shift once densities surpassed 2.0 animals km⁻², where the density threshold is a function of competition for space instead of available forage, but Hinkes et al. (2005) found variable support for this hypothesis among three caribou herds from southwest Alaska. Range expansion of Adak caribou does not appear to be linked to a set density threshold since periods of expansion during the 1980s, 1990s, and 2000s correspond to densities of approximately 0.5, 1.6, and 5.0 caribou km⁻², respectively. This is not necessarily surprising owing to the spatial constraints imposed by insularity compared to much larger movement opportunities for caribou on the mainland.

Adak caribou are quite capable of traversing the minimum 30‐km width of the island within a day's travel (Skoog 1968, Bergerud 2000), yet their range was largely restricted to the west side of the island until the early 1980s. Why did they stay on the west side for so long? Unlike the pattern observed for numerical population dynamics, the initially slow range expansion may represent a type of latency. According to early reports, the original population presumably stayed within 6 km² of the release site for four years before venturing toward the center of the island (Jones 1966; R. Thabes, unpublished manuscript). The original animals were pen‐raised as calves and released after two months of captivity (Jones 1966). Without established adult groups to follow, these animals likely did not begin extensive exploration of their new habitat before becoming mature adults. It is also important to note that while sporadic sightings of caribou in the Kagalaska Strait region of Adak (Fig. 1) occurred during the 1970s (the time frame when our range estimates do not extend eastwardly), biologists at the time noted that “from most reports, very few animals are presently in the eastern half of the island” (T. Early, unpublished manuscript).

While consistent inter‐annual use of calving areas has been documented for a number of mainland herds with varying densities over time (Skoog 1968, Gunn and Miller 1986), shifts to a new calving area have occurred for some mainland herds in association with increasing population size (Taillon et al. 2012), interchange of individuals between herds (Hinkes et al. 2005), and anthropogenic disturbance (Wolfe et al. 2000, Cameron et al. 2005). Island caribou, however, inhabit a unique environment free from predation pressure, inter‐specific competition, and new gene flow. Harassment from mosquitos (Culicidae), black flies (Simuliidae), and bot flies (Oestridae) can also drive seasonal movement patterns for mainland herds (Skoog 1968), but biting insects are of negligible importance on Adak where only white sox (a Simuliid species) are active on rare calm days. In addition, the original released population successfully conceived as yearlings and produced calves as two‐year olds in 1960 when animals were reported to stay close to the release site on the west side of the island (Jones 1966). Thus, a founder effect characterized by animals repeatedly returning to the natal sites of the original population on the west side of Adak is not apparent. However, caribou began to calve on the south‐central side of Adak within no more than seven years after the first non‐captive calves were born, and have shown fidelity to that site thereafter.

We propose that calving on the south‐central side of Adak likely became established once females located an area with high quality resources that coincided with parturition. The calving area has an open bowl shape comprised of mountainous terrain to the west, north, and east draining into a large basin (Teardrop Basin) with a south facing aspect (Fig. 1) where snow melts more quickly. In turn, the topographic relief and southern exposure may promote spatial variation in plant phenology during the early growing season that correlate positively with calf production (Post et al. 2008), and create conditions conducive for large congregations of calf‐cow groups typical of those encountered during aerial surveys after the peak of calving. As the summer growing season progresses, animals may disperse more widely to track phenological changes in other habitats, which can explain the rapidly expanding and non‐spatially restricted post‐calving areas and isometric relation between density with utilization distributions and total population size. Observed calving site fidelity could also be influenced by evolved anti‐predation mechanisms that have not had adequate time to break down (Berger 1999, 2007). Anti‐predation mechanisms can relax substantially after 10–20 generations of caribou have been removed from predation risk (Berger 2007), but at least 10 generations have already elapsed on Adak since the introduction. Breakdown of site‐fidelity could be expected in future years, but completely disentangling outcomes based on the resource availability or predation hypotheses is untenable without knowledge of changing resource availability or adding pressure from predation.

Management implications

The removal of the Adak herd was identified as a preferred management action at the end of the military period (USFWS and ADFG 1994), but the challenges of mixed land ownership largely prevented plan implementation. Notably, introduced reindeer have persisted without intensive harvest for over 100 years on the Aleutian Islands of Atka and Umnak, and these herds may provide an analog for an unmanaged Adak herd. While informative time series counts do not exist for Atka and Umnak, these populations are suspected to have undergone several irruptive oscillations fitting the Riney‐Caughley model (Swanson and Barker 1992, USFWS and ADFG 1994), and at least in the case of Atka, may represent a herbivory‐induced alternative ecosystem state (Ricca 2013). It is also clear that sustained harvest pressure is a critical tool for managing the magnitude of irruptive population and spatial dynamics of island caribou. A basic model where λ = 1 + (recruitment rate − mortality rate) supports our contention. If we let recruitment rate simply equal the proportion of calves in the total 2012 count (0.11) and the mortality rate equal the average reported harvest from 2005 to 2011 (0.09), then λ = 1.02. So, assuming only hunting related mortality and continued density‐dependent reductions in productivity (as our simulated high H + DD model suggests), current harvest rates may be able to dampen future oscillations in population growth. Moreover, our discovery of calving (albeit small) on Kagalaska provides further evidence of Adak caribou evolving into a meta‐population (Ricca et al. 2012). Kagalaska now represents a new stage 1 dispersal front, which also presents new challenges for managers charged with controlling invasive species.

Given the rather unique combination of hunter harvest, mild climate, and human habitation and infrastructure on Adak, we caution that our results may not apply directly towards the management of all Rangifer populations introduced to islands at mid and high latitudes. Nevertheless, our study illustrates how the Riney‐Caughley model remains a useful paradigm to broadly describe the population and spatial dynamics of introduced island ungulates.

Acknowledgments

We recognize the many biologists, island residents, and pilots whose collective efforts over 50 years produced much of the time series used in this manuscript, including (but not limited to) D. Boone, D. Brown, V. Byrd, G. Brogdon, C. Dean, F. Dienes, T. Early, S. Ebbert, J. Fischer, J. Fuller, K. Hall, R. Jones, E. Klett, K. Livezey, J. Martin, M. Masteller, J. Meehan, J. Mueller, B. Reiswig, M. Smith, L. Slater, L. Spitler, T. Spitler, B. Springer, G. Tarbet, R. Thabes, and V. Tutiakoff Jr. This projected was supported in part by funding from the USFWS Region 7 Invasives with Volunteers Program, USGS Western Ecological Research Center, USFWS‐AMNWR, Friends of Alaska National Wildlife Refuges, and Texas State University Department of Biology. Helicopter use for the 2012 and other surveys was graciously facilitated through cooperative agreements between the USGS Alaska Volcano Observatory (J. Paskievitch) and USFWS‐AMNWR. M. Ricca, F. Weckerly, and J. Williams (co‐authors) along with V. Tutiakoff Jr. conducted the 2012 survey. We thank A. Duarte, V. Eviner, G. Herring, D. Klein, T. Kimball, and an anonymous reviewer for helpful comments on prior manuscript drafts. J. Estes and T. Tinker kindly allowed use of their skiff to access remote areas for ground counts, and conversations with Adak residents and hunters were very informative. Use of product names does not imply endorsement by the U.S. Government.

Supplemental Material

Appendix