The rise of an apex predator following deglaciation

Sea otters (Enhydra lutris) are an apex predator of the nearshore marine community and nearly went extinct at the turn of the 20th century. Reintroductions and legal protection allowed sea otters to re‐colonize much of their former range. Our objective was to chronicle the colonization of this apex predator in Glacier Bay, Alaska, to help understand the mechanisms that governed their successful colonization.


| INTRODUC TI ON
Apex predators have experienced substantial declines worldwide due to human impacts including harvest and habitat loss (Estes et al., 2011;Prugh et al., 2009;Ripple et al., 2014). Loss of apex predators has resulted in global change in the structure of communities, leading to ecological, economic and social impacts around the world (Prugh et al., 2009;Simenstad, Estes, & Kenyon, 1978). Recently, the generality of these impacts have been documented (Estes et al., 2011;Ripple et al., 2014;Ritchie & Johnson, 2009), and reintroductions of locally extirpated apex predators by translocation have become a management tool for restoring community assemblages (e.g., rewildling; Smith, Peterson, & Houston, 2003;Enserink & Vogel, 2006;Palomares, Rodriguez, Revilla, Lopez-Bao, & Calzada, 2011). Ecosystem transitions resulting from loss of apex predators can be difficult or impossible to reverse (Estes et al., 2011), and the success of reintroducing apex predators into ecological communities following their extirpation has been mixed and is often confounded with human-induced changes to the landscape (e.g., habitat fragmentation, climate change; Crooks & Soulé, 1999;Schadt et al., 2002;Steury & Murray, 2004;Ritchie et al., 2012). There are few documented cases that chronicle the colonization dynamics of an apex predator following extirpation, where direct human interactions have been limited. The sea otter provides one of the few exceptions, as remnant and reintroduced populations across the North Pacific have been the focus of long-term studies of population dynamics and community ecology (Bodkin, Ballachey, Cronin, & Scribner, 1999;Estes, 1990;Estes & Palmisano, 1974;Kenyon, 1969;Larson, Bodkin, & VanBlaricom, 2015).
Sea otters are an apex predator of the nearshore marine food web whose presence alters community dynamics through a trophic cascade of reduced prey (e.g., sea urchins), and resultingly, more macrophytic marine vegetation (e.g., kelp forest; Estes & Palmisano, 1974;Estes, Smith, & Palmisano, 1978;Duggins, 1980;Estes & Duggins, 1995). An apparent exception to the sea otter role as an apex predator is in the Aleutian Archipelago, Alaska, where sea otter abundance has recently been regulated by killer whales (Orcinus orca; Estes, Tinker, Williams, & Doak, 1998).
Sea otters were overharvested rangewide during the multi-national commercial fur trade in the 18th and 19th centuries (Bodkin, 2015;Kenyon, 1969). By the end of the 19th century, sea otters were extirpated from southeastern Alaska, and nearly all of the North Pacific, with only a few small isolated populations persisting (Bodkin, 2015;Kenyon, 1969). Legislation following the maritime fur trade, including the International Fur Seal Treaty (1911), the Marine Mammal Protection Act (1972) and the Endangered Species Act (1977), provided legal protection to sea otters from most harvest. These efforts, combined with translocation by humans, have resulted in the return of sea otters to much of their former distribution, after being absent for periods ranging from decades to centuries (Bigg & MacAskie, 1978;Bodkin, 2015;Doroff, Estes, Tinker, Burn, & Evans, 2003;Jameson, Kenyon, Johnson, & Wight, 1982;Kenyon, 1969;Lubina & Levin, 1988).
By the time sea otters were extirpated from southeastern Alaska, glaciers in Glacier Bay had already begun retreating. Glacier Bay is a tidewater glacier fjord located in southeastern Alaska, USA (58.67°N, 136.90°W). Since 1750, the glaciers in Glacier Bay have retreated >100 km inland at a mean speed of 0.4 km per year, representing the most rapid and extensive retreat in modern times, 15 times faster than any other recorded tidewater glacier (Figure 1; Lawrence, 1958;Chapin, Walker, Fastie, & Sharman, 1994;Fastie, 1995).
The ecological succession of the nearshore community in Glacier Bay occurred in the absence of sea otters as an apex predator for >250 years following deglaciation. By the time sea otters returned to southeastern Alaska, and subsequently arrived near Glacier Bay in the late 1980s (Figure 1), the nearshore ecosystem supported a diverse and abundant community of benthic invertebrates, including abundant populations of sea otter prey items (Weitzman, 2013).
The spread of sea otters into Glacier Bay provides a unique opportunity to examine the potential of an apex predator to colonize vacant habitat in an area (a) where ecological succession following deglaciation occurred in the absence of apex predation, (b) where the apex predator is legally protected from harvest, and (c) that has relatively few direct impacts from humans. Using novel statistical tools for estimating spatio-temporal colonization dynamics, we chronicled the change in distribution and abundance of sea otters in Glacier Bay. We were focused on the twenty-year period starting when sea otters first began consistently using Glacier Bay in 1993, through 2012, a time that we estimated most of the bay had been colonized by sea otters. Finally, we discuss potential mechanisms for the extraordinary success of sea otters in Glacier Bay. distribution and abundance, reshaping future marine communities in the wake of deglaciation and global loss of sea ice.

| Data collection
A variety of data sources were collected to chronicle the distribution, abundance and colonization dynamics of sea otters in Glacier Bay. We leveraged information from three data sets collected between 1993 and 2012. The data sets included (a) design-based aerial surveys, (b) intensive search units to estimate detection probability during design-based aerial surveys and (c) distribution aerial surveys, each described in the next three sections. Other opportunistic observations of sea otters in Glacier Bay are reported in Figure 1.

| Design-based aerial surveys
The first data set we used included probabilistic, design-based aerial surveys. The design-based aerial surveys had the most rigorous sampling design of the three data sets (Bodkin & Udevitz, 1999), but were concentrated between 1999 and 2006. Therefore, they did not cover the early colonization process that started in the early 1990s. The survey methods are described in detail in Bodkin and Udevitz (1999) and Esslinger, Esler, Howlin, and Starcevich (2015).
Briefly, this survey consisted of an observer flying in a single-engine fixed-wing airplane following pre-determined linear transects placed systematically across Glacier Bay, with a randomly selected location for the initial transect. The transects were stratified based on two criteria, ocean depth and distance from shore. Areas with depth < 40 m received higher sampling effort than areas with depth >40 m, and areas closer to shore received higher sampling effort. Transects were flown in years 1999-2004, 2006 and 2012. Transects were 400 m wide, indicated by strut marks on the aircraft, flown at a velocity of 29 m/s and a height of 91 m. Observers searched for, and located sea otters within transects, and subsequently counted them. Data collected during these surveys and the associated intensive search units (described in the next section) have been paired with design-based estimators of abundance and uncertainty that account for total area surveyed (Bodkin & Udevitz, 1999;Esslinger et al., 2015). Abundance estimates and 95% confidence intervals from these data and the associated design-based estimators are provided in Figure 2.

| Intensive search units
The design-based surveys undercounted sea otters due to imperfect detection (e.g., diving sea otters are often undetected; Williams, Hooten, Womble, & Bower, 2017a). To estimate detection probability, additional data were collected during the surveys in intensive search units (ISUs). ISUs were a randomly selected subset of 469 groups of sea otters from the design-based survey. Each group of ≥1 sea otter observed during the design-based surveys received a group identification number (group ID). The ISUs occurred on a pre-determined random subset of group ID, so ISUs would not occur over areas with no sea otters present. ISU data were collected every year the design-based surveys occurred, with effort distributed approximately equally among years. At these 469 random sites, after a group of ≥1 sea otters was detected and counted using the procedures from the design-based survey, the pilot deviated from the original design-based transect and circled the group of sea otters.
The pilot flew five concentric circles around the group of sea otters so observers could obtain accurate counts of abundance within the group. The five circles were flown in approximately 3.6 min, a time chosen based on the aerobic dive limit of sea otters (so diving sea otters could be included in the counts when they resurfaced).

| Distribution surveys
Sea otter distribution surveys were conducted in Glacier Bay by one or more observers in a single-engine high-wing airplane. Distribution surveys had the least rigorous sampling design, but were conducted during early years of colonization, and provide information on the initial colonization process. They were also conducted in years when the design-based surveys were not conducted, providing information on what occurred between design-based surveys. In an attempt to survey all shoreline habitat <40 m deep, swaths were flown parallel to the shoreline at an altitude of 152 m during calm sea conditions. Distribution surveys used in our analysis were flown in 1993, 1994-1998, 2009 and 2010. ISUs were not flown during distribution F I G U R E 2 Abundance estimates of sea otters in Glacier Bay National Park, southeastern Alaska. Model-based estimates are the posterior mean and 95% credible intervals from our spatiotemporal integrated population model. Design-based estimates are the estimated mean and 95% confidence intervals from the design-based surveys that occurred during 1999-2004, 2006 and 2012. Count data are the number of sea otters counted during the distribution surveys that occurred during 1993, 1995-1998, 2005, 2009 and 2010. Data from both the design-based surveys and the distributional surveys were used to fit the spatio-temporal integrated population model Year N Model−based estimates Design−based estimates Counts from distributional surveys surveys. Total counts of sea otters obtained during distribution surveys are provided in Figure 2 (blue points). All known occupied habitats were surveyed, but results were not corrected for areas not surveyed or for individuals not detected, and thus may not reflect actual abundance. Instead, distribution survey data provide information on minimum number of sea otters present, additional information on the spatial distribution of sea otters and relative abundance of sea otters.

| Statistical analysis
Previous analysis of the design-based data described above has provided robust and unbiased estimates of abundance of sea otters in Glacier Bay (e.g., Figure 2; Esslinger et al., 2015). However, inference from design-based estimators is often limited to annual abundance estimates (and their associated uncertainty) and provides limited information on spatial or spatio-temporal processes such as colonization dynamics. Therefore, we used a Bayesian hierarchical spatio-temporal integrated data model to simultaneously estimate occupancy, abundance and colonization dynamics of sea otters in Glacier Bay. The general model that we adapt was originally described in Williams et al. (2017b). We tailored this model to our analysis, including rigorous model selection and model-checking procedures, and additional derived quantities of interest for characterizing a colonizing population. We used R statistical software (R Core Team 2013) to conduct all analyses.

| Model specification
We represent counts of sea otters obtained during design-based and distribution aerial surveys as y i,t , at location i = 1, …, q during time t = 1993, …, 2012. Sea otter counts under-represent the true number of sea otters present (a latent parameter represented with n i,t ) due to incomplete detection of individuals. We modelled the relationship between counts and the true abundance using where p t is the detection probability during year t. Equation 1 is a special case of a more general data model that permits estimation of both group and individual level detection probability. Specifically, Equation 1 assumes p t is the probability of detecting an individual sea otter, and the probability of detecting a group of sea otters, in Supporting Information for additional information). We allowed detection probability to vary by year to account for differences in observers and survey conditions. The ISU data provided information for 469 values of n i,t , and allowed us to estimate p t for each year ISU data were collected. For all n i,t for which we did not have data, we modelled it as a negative binomial random variable with expected abundance λ i,t , and overdispersion parameter τ. We used a negative binomial distribution (vs. a Poisson distribution, for example) to account for possible overdispersion in the data (Ver Hoef & Boveng, 2007) and assessed this choice using a cross-validation model-checking procedure (Conn, Johnson, Williams, Hooten, & Melin, 2018).
We modelled the spatio-temporal dynamic process of sea otter spread using a continuous-space, continuous-time reaction-diffusion model. Reaction-diffusion models are commonly used in theoretical ecology (e.g., Cantrell & Cosner, 2003;Cressie & Wikle, 2011;Hooten & Wikle, 2008;Okubo, 1980;Wikle & Hooten, 2010) to model animal movement and spatio-temporal population processes (Holmes, Lewis, Banks, & Veit, 1994), but are often deterministic, and not statistical in nature. We used our data to estimate model parameters, and quantified parameter uncertainty in a Bayesian hierarchical framework, permitting us to couple mechanistic and mathematical models of animal movement and growth with statistical uncertainty (Wikle & Hooten, 2010).
We aligned the discrete nature of the data collection process, where data were collected in 400 × 400 m grid cells, with our continuous-space, continuous-time model, by integrating over grid cells where λ(s, t) is the intensity at any location s = (s 1 , s 2 )′ (e.g., s 1 = latitude, s 2 = longitude) during any time t, and λ i,t is the expected abundance of sea otters in a grid cell. We used a partial differential equation called Ecological Diffusion for our reaction-diffusion model (Garlick, Powell, Hooten, & McFarlane, 2011;Hefley, Hooten, Russell, Walsh, & Powell, 2017;Hooten, Garlick, & Powell, 2013;Turchin, 1998;Williams, Hooten, Womble, Esslinger, & Bower, 2018;Williams et al., 2017b). Ecological diffusion is a mechanistic diffusion model that emerges from individual random walks with heterogeneous movement probabilities based on local environmental information, resulting in variable motility rates (see Hefley et al., 2017, for a thorough discussion of ecological diffusion). The term (s,t) t in Equation 2 represents the instantaneous change in abundance intensity (i.e., change in mean abundance) as a function of intrinsic population growth and local immigration and emigration. The diffusion parameter δ(s) represents sea otter motility and depends on spatial exogenous covariates. Because sea otters are influenced by characteristics of the environment and move slowly through areas that provide necessary resources and quickly through areas that do not, allowing motility to vary based on environmental characteristics provides a mechanism for estimating resource selection for a colonizing species. Motility has physically interpretable units (km 2 year −1 ) that are inversely related to residence time (km −2 year; Hefley et al., 2017). We consider a log-linear relationship between motility δ(s) and four spatial covariates that we expected to be associated with sea otter resource selection, based on previous research.
We used a log-linear relationship because expected residence time (and therefore motility) must be positive (Hefley et al., 2017 We centred and scaled all covariates except depth, which was an indicator variable. We expected that sea otters would use shallow areas, close to shore, with steep bottoms, and relatively simple shorelines, and therefore, diffusion out of sites with these traits would occur slower than at sites not having some or all of these traits. We used the interaction slope(s) × depth(s) in Equation 3 because we hypothesized that slope of the ocean floor was only important in shallow areas (i.e., <40 m, depth covariate = 1) where the ocean floor was readily accessible to sea otters. We hypothesized that sea otters preferred simple shoreline complexity based on previous observations of large groups of sea otters in Glacier Bay around islands or linear shorelines. Based on these hypotheses, we expected β 1 to be negative, reducing diffusion out of shallow areas (where depth = 1), β 2 to be positive, increasing diffusion as distance from shore increased, β 3 to be negative, reducing diffusion as bottom slope increased in shallow areas and β 4 to be positive, increasing diffusion in areas with a large shoreline complexity index.
The reaction parameter γ in Equation 2 represents instantaneous growth rate of sea otters, excluding local immigration and emigration. Thus, the function γλ(s, t) in Equation 2 represents Malthusian growth. Malthusian growth for the Glacier Bay sea otter population during our time frame was supported by abundance estimates (Figure 2; red points and lines) obtained from the survey data we described in Design-based aerial surveys, above. We assumed the instantaneous growth rate was constant among sites and years. In principle, Equation 2 permits instantaneous growth rates to vary in space and time. However, data required to identify spatio-temporal differences in growth rates were exceedingly large for our application, and therefore, we estimated a constant (average) growth rate. Finally, the second derivatives in space model the effect that varying spatial densities have on local abundance due to random motion.
The PDE in Equation 2 requires specification of an initial condition, an estimate of the abundance intensity λ(s, t) for t ≤ 1993. We where θ is a scale parameter controlling the initial density, κ is a dispersion parameter controlling the radial distance of the initial density, and d is the epicentre where sea otters were detected in 1993. We used a reflective spatial boundary condition for λ(s, t) at sites adjacent to terrestrial environments, because sea otters are typically restricted to marine environments (see Appendix S1 in Supporting Information for additional details).
To complete the Bayesian hierarchical specification of our model, we assigned prior models for β, γ, p t , θ, κ, and τ. We specified β ~ Normal(0, σ 2 I), where I is the identity matrix. We selected σ 2 optimally, using regularization (see Parameter Estimation, Model Selection and Model Validation, below). We assumed that the instantaneous growth rate of sea otters in Glacier Bay was between 0 and 1, and therefore specified γ ~ Beta(1,1). We specified the prior for detection probability as p t ~ Uniform (0.5,0.95). We expected to learn about annual detection probability via our posterior distributions during years with ISU data (i.e., the years the design-based surveys occurred). However, we did not expect to learn about detection probability during years in which only the distribution surveys were conducted (i.e., the years when ISU data were not collected) and expected our posterior distributions to reflect our prior distributions. We specified θ ~ Normal + (500, 200 2 ), and κ ~ Normal + (10, 100 2 ), for our initial condition parameters, where Normal + indicates the zero-truncated normal distribution. Finally, we specified the overdispersion parameter τ ~ Uniform (0,1000). We chose 1000 as the upper bound of the uniform distribution so that it was sufficiently large to be greater than any realistic value for τ.

| Derived parameters related to colonization and distribution
We used our model to derive parameters related to sea otter colonization dynamics. These included occupancy dynamics and the colonization front. We defined occupancy probability ψ i,t as the probability the latent true abundance n i,t was greater than zero.
Given n i,t ~ NB(λ i,t , τ), the occupancy probability is We estimated the posterior distribution of the colonization front by calculating the probability each site i would be at the front of colonization in year t. We defined the front of colonization, χ t for t = 1993, …, 2012, as the occupied site s i |n i,t > 0 that maximized the distance between the site location s i and epicentre d, χ t = max|s i,t -d|, for all i that n i,t > 0. In cases where land obstructed the path between (3) log ( (s)) = 0 + 1 depth(s) + 2 dist(s) + 3 (slope(s) × depth(s)) + 4 shore(s).
(s,t = 1993) = e −|s−d| 2 2 a site and the epicenter, we measured distance using the shortest path not obstructed by land. Because n i,t is a random quantity, the colonization front is random, and the probability that the ith site was the colonization front was calculated using for k = 1,…, K Markov chain Monte Carlo (MCMC) iterations, and where I { (k) t =i} is an indicator variable that equals one when (k) t = i, and zero otherwise.
There is rich literature in theoretical ecology related to estimating spread rates of biological invasions using reaction-diffusion equations similar to equation 2 (e.g., Andow, Kareiva, Levin, & Okubo, 1990;Fisher, 1937;Shigesada & Kawasaki, 1997;Skellam, 1951). In a homogeneous environment (i.e., δ(s) and γ are constant), spread rates approach 2 √ (Fisher, 1937). Similarly, evidence suggests that spread rates approach 2 √̄̄, asymptotically, when ̄ and ̄ are properly homogenized values of diffusion and growth rates (see Appendix S1 in Supporting Information). We calculated the spatially explicit asymptotic spread rate using means of the posterior distributions for the parameters ̄ and ̄, for comparison to our Bayesian approach described in the previous paragraph and Figure 5. A map of spatially explicit asymptotic spread rates (and a histogram of values) is shown in Figure 6.

| Parameter estimation, model selection and model validation
We computed posterior distributions of model parameters using MCMC methods. For each model fit, we obtained 3 chains of 40,000 MCMC iterations and visually assessed the chains for convergence. We used Bayesian regularization and K-fold cross-validation for model selection to optimally (under the mean absolute difference score function) select the shrinkage parameter, σ 2 , controlling the relative importance of each environmental covariate on motility (Appendix S1 in Supporting Information). We retained the value of σ 2 and the associated parameter estimates from our top model for inference. We used Bayesian p-values for model validation to compare observed data (i.e., data withheld from model fitting) to model-predicted observations sampled from the posterior predictive distribution (Conn et al., 2018;. Nearly all of the observed data (i.e., >99%) fell within the 95% credible intervals of the posterior predictive distributions from our final model, suggesting no lack of model fit. Maps comparing model-based estimates of mean occupancy probability and occupancy data collected during design-based surveys are provided in Appendix S1 in the Supporting Information.

| RE SULTS
The estimated mean intrinsic rate of population growth was 20% per year ( Figure 2, Table 1 Our analyses suggested that sea otters in Glacier Bay preferred shallow areas, close to land, with steep bottom slope, and a relatively small shoreline complexity index (e.g., islands; Figure 4, Table 1), moving quickly throughout Glacier Bay to locate these areas. Our regularization and model selection procedure effectively "shrinks'' estimates of β 1 , …, β 4 with low predictive power to zero . Thus, because all covariates are centred and scaled, larger values of β indicate relatively higher predictive power of the associated covariate, compared to the other covariates (Table 1, Supporting information Figure S1.1).
In 2012, sea otters were observed occupying areas near the Gilbert Peninsula, in mid-Glacier Bay (Figure 1), about 54 km from where they were observed in 1993, suggesting a minimum mean TA B L E 1 Posterior mean and 95% credible intervals for estimated parameters in the sea otter occupancy-abundance model. Covariates associated with each parameter were centered and scaled, except depth (β 1 ) which was an indicator variable equal to one when depth <40 m and zero otherwise  1993 1994 1995 1996 1997 1998 1999

| Modelling the successful return of sea otters to Glacier Bay
Modelling dynamic ecological processes, such as the colonization of an apex predator, using mechanistically motivated PDEs, provides scientific insight that would not be possible using only design-based estimators of abundance (Hefley et al., 2017;Williams et al., 2018). For example, population growth rates of sea otters in Glacier Bay estimated using design-based estimates of abundance reported in Figure 2 suggest annual growth rates of 44% per year between 1999 and 2012 (Bodkin, 2015;Esslinger & Bodkin, 2009). The maximum theoretical growth rate of sea otters is between 20% and 24% (Estes, 1990), suggesting a larger growth rate requires both intrinsic population growth within Glacier Bay and immigration from outside the bay (Bodkin, 2015;Esslinger & Bodkin, 2009). Further, the estimated 44% growth rate does not address uncertainty in abundance estimates, is sensitive to the estimated initial population size and does not leverage other sources of data (e.g., the distribution data). By explicitly modelling both intrinsic population growth and movement of sea otters, as well as initial population size using multiple sources of data, we estimated that sea otter intrinsic growth rate (excluding immigration and emigration) in Glacier Bay was 20% between 1993 and 2012.
By comparison, intrinsic growth rate estimates of northern sea otter populations, below equilibrium density, range from about 8-24% Estes, 1990).
An intrinsic growth rate of 20% resulted in annual abundance estimates that aligned closely with the design-based abundance In a 20-year period, sea otter distribution in Glacier Bay went from non-existent, to colonizing nearly the entire bay at a mean maximum spread rate ranging from six km/year early during colonization, to one km/year near the end of colonization (Figure 4,5).
The asymptotic spread rates calculated using 2 √̄̄, where ̄ and ̄ were the estimated posterior means of the homogenized values of diffusion and growth rate, ranged from 1.5 to 4.5 km/ year ( Figure 5). For comparison, Lubina and Levin (1988) reported spread rates of California sea otters between 1.74 km/year and 3.5 km/year, and in a terrestrial system where apex predator colonization rates have been examined, wolves (Canis lupus) recolonized the Greater Yellowstone Ecosystem at a rate of 9.78 km/ year (Hurford, Hebblewhite, & Lewis, 2006;Smith, Murphy, & Guernsey, 2000). Differing rates of colonization spread through time ( Figure 5,6)

| Why were sea otters so successful in Glacier Bay?
By 2012, sea otters were one of the most abundant marine mammals in Glacier Bay (see Mathews & Pendleton, 2006;Mathews et al., 2011;Womble et al., 2010;Gabriele et al., 2017, for estimates of other marine mammal abundance). Ecological theory posits that three attributes contribute to the relative success of an invading or (re)colonizing species: (a) availability of resources, (b) natural predators and (c) the physical environment (Shea & Chesson, 2002). We discuss these three attributes and their relation to sea otter colonization in Glacier Bay, in the following paragraphs.

| Resource availability
We found that residence time, and therefore occupancy and abundance, was positively correlated with areas <40 m in depth, close to shore, had steep bottom slopes, and a relatively small shoreline complexity index. These physical attributes were likely associated with access to abundant prey that sea otters were able to exploit upon their arrival, suitable resting habitats, in addition to facilitating sea otter social behaviour (Kenyon, 1969;Monson, Estes, Bodkin, & Siniff, 2000). Recent deglaciation and subsequent ecological succession of Glacier Bay occurred in the absence of predation pressure by sea otters, resulting in abundant invertebrate prey that fuelled sea otter population growth.
Sea otter prey availability was examined in Glacier Bay during early stages of colonization, and in nearby control sites outside Glacier Bay (Bodkin et al., 2001). Control sites had been occupied by sea otters for ≥10 years. Sea otters' preferred prey were clams, and clam density was 3-10 times higher in lower Glacier Bay than in control areas, clam biomass was 5-12 times higher in lower Glacier Bay than in control areas, and mean size of preferred clam species was nearly twice as large within Glacier Bay compared to control areas. These results indicated that the mid-trophic-level consumers preferred by sea otters grew in size, abundance and distribution in the absence of apex predation, facilitating colonization when sea otters arrived. It is also evident that, had sea otters been present in southeastern Alaska during deglaciation and subsequent ecological succession, there would not have been the extended period of predator release, and the resulting benthic invertebrate guild may have been much different than what was present in the late 20th century.

| Absence of natural predators and human harvest
Resource availability was not the only factor responsible for sea otter growth rates observed in Glacier Bay. Densities of subtidal clams in other areas across southeast Alaska have been documented as being comparable to those observed in Glacier Bay (Kvitek & Oliver, 1992), but these areas have generally demonstrated more modest population growth rates (Esslinger & Bodkin, 2009). Similarly, Estes et al. (1998) observed increasing prey availability (i.e., sea urchin density) concurrent to sea otter population declines in western Alaska. In both these cases, sea otter populations were regulated by top-down mechanisms, from harvest by humans in the first case (Bodkin et al., 2007) and from predation by killer whales in the second case (Estes et al., 1998). Sea otters are harvested outside of Glacier Bay, but are protected from human harvest within Glacier Bay. Further, although killer whales occurred in Glacier Bay during our study and have been observed preying on sea otters throughout much of their range (Hatfield, Marks, F I G U R E 6 Left: Spatially explicit asymptotic spread rates (km/year) calculated using 2 √̄̄, where ̄ and ̄ were the estimated posterior means of the homogenized values of diffusion and growth rate, respectively (see Appendix S1 in Supporting Information). The dimension of each grid cell in the map is 4000 × 4000 m 2 , and coordinates are in metres ( Tinker, Nolan, & Peirce, 1998), there was no evidence that sea otters were limited by killer whale predation in Glacier Bay (Matkin, Straley, & Gabriele, 2007). Thus, in the absence of human interference and natural predators, and the presence of sufficient resources, apex predators can be highly successful in colonizing areas in which they were formerly absent or previously extirpated.

| Physical environment
Changes in the physical environment (due to changes in land use or global climate change) are often responsible for decreased apex predator abundance (Brashares, Prugh, Stoner, & Epps, 2010 (Comiso, 2002), and several marine species that seasonally occupy Arctic and subarctic habitats are likely to encroach into northerly latitudes, and compete with extant Arctic species (Moore & Huntington, 2008). Our results indicate that sea otters are capable of range expansion in the presence of newly available habitat, and could expand to more northern latitudes following the loss of perennial sea ice, and subsequently, reshape the nearshore marine ecosystem.

| Future change to the nearshore marine ecosystem in Glacier Bay
Colonization of keystone apex predators into communities and associated food webs often pushes ecosystems to alternative ecosystems states, a phenomenon that is well documented for sea otters (Estes et al., 2011). In rocky and soft-sediment marine communities, sea otters generally have a strong, top-down influence on infaunal prey communities and kelp forests (Estes & Duggins, 1995;Kvitek & Oliver, 1992;Weitzman, 2013). The growth and expansion of sea otters within Glacier Bay has resulted in the decline in density and abundance of sea otter prey species in several areas (Weitzman, 2013). A reduced prey base affects other predators including octopus, sea stars, fishes, birds and mammals, which can trigger ecosystem level changes (Bodkin et al., 2001). It is clear that future change in the Glacier Bay ecosystem is forthcoming, and sea otter abundance will reach a threshold, at which point they will face regulating forces. The question remains as to which force will regulate the sea otter populations in Glacier Bay? There is evidence that sea otters have already altered prey abundance (Bodkin et al., 2001;Weitzman, 2013), suggesting resources may become limited, but the pattern of evidence elsewhere in Alaska suggests apex predation by killer whales can severely limit sea otter populations. Perhaps the most likely scenario, supported by historical evidence of sea otters across their range, combined with local knowledge of resource availability surrounding Glacier Bay, is that, after prey becomes limiting, increased emigration and/or decreased survival will result in reduced abundance of sea otters in Glacier Bay.
Critical to the aetiology of conservation biogeography is an understanding of past and present ecosystem states. We provide information chronicling the spatio-temporal population spread of sea otters in an ecosystem that was completely covered by glacial ice analysed the data and led the writing. All authors contributed to manuscript revisions and gave final approval for publication.

S U PP O RTI N G I N FO R M ATI O N
Additional supporting information may be found online in the Supporting Information section at the end of the article.