Viability of a small, geographically-isolated population of Beluga Whales, Delphinapterus leucas: effects of hunting, predation, and mortality events in Cook Inlet, Alaska.
Beluga whales, Delphinapterus leucas, in Cook Inlet (lat. 59[degrees]-61.5[degrees]N, long. 149[degrees]-154[degrees]W), Alaska, make up a small, genetically distinct population that appears to have strong site fidelity to the inlet year-round (O'Corry-Crowe et al., 1997, 2002; Rugh et al., 2004; Hobbs et al., 2005; Shelden et al., 2015). The Cook Inlet beluga whale (CIBW) population declined dramatically in the 1990's, and continued a steady decline through 2012 (Hobbs et al., 2015a), raising concern about its risk of extinction and subsequently leading to listing this distinct population segment as endangered under the U.S. Endangered Species Act (ESA) in October 2008 (NOAA, 2008). Should the CIBW population go extinct, the closest population in Bristol Bay (BB) is 1,500 km away by sea and separated by the Alaska Peninsula that extends three degrees of latitude south of the southern limit of the bay. Similar to CIBWs, satellite-tagged BB whales remained within the bay year-round (Citta et al., In press), thus, it is highly unlikely that beluga whales would repopulate Cook Inlet in the foreseeable future. Extinction of the CIBW population would result in a permanent loss of range for the beluga whale species.
Alutiiq Eskimos and Dena'ina Athabaskan Indians have occupied the coastal areas surrounding Cook Inlet since prehistoric times (de Laguna, 1975). These hunting societies utilized many marine resources including beluga whales. During the 20th century, the Dena'ina in Tyonek (a small village on the west side of Cook Inlet) and Eskimo whalers from communities outside of Cook Inlet hunted beluga whales for subsistence, and there was also periodic sport hunting and large-scale commercial hunts by non-Native hunters (Mahoney and Shelden, 2000). Commercial and sport hunting ended with the introduction of the U.S. Marine Mammal Protection Act (MMPA) in 1972. Subsistence hunting by Alaska Natives was allowed under the MMPA and monitoring by the Alaska Department of Fish and Game (ADFG) indicated that the practice continued from 1972 to the present (Mahoney and Shelden, 2000). During the 1960's, 70's, and early 80's, ADFG conducted a number of aerial surveys that covered parts of Cook Inlet, documenting distribution and numbers of beluga whales (Shelden et al., 2015). The highest estimate of 1,292 beluga whales reported during these surveys was based on counts made in August 1979 (Calkins (1)).
In June 1993, the National Marine Fisheries Service (NMFS) began comprehensive, systematic aerial surveys of the beluga whale population in Cook Inlet (Rugh et al., 2000, 2005, 2010; Shelden et al., 2013, 2015). Survey results showed a decline in abundance of nearly 50% between 1994 and 1998, from an estimate of 653 whales to 347 whales (Hobbs et al., 2000a). Concern over the high level of human-caused mortality on this whale population prompted NMFS to designate it as depleted under the MMPA (NOAA, 2000) and to regulate the Native Alaskan subsistence hunt (Mahoney and Shelden, 2000).
With a limited hunt between 1999 and 2014 (a total of five whales taken), it was anticipated that the population would begin to recover though this has not been the case (Hobbs et al., 2015a). Hunters indicated a preference for large, white-skinned beluga whales (presumably adults) (Mahoney and Shelden, 2000; Huntington, 2000); therefore, it is conceivable that the population had a deficit of reproductive-age females, which could be inhibiting recovery through reduced calf production (Hobbs et al., 2015b). However, the hunt may not be the only risk factor behind the continued decline of this population.
In addition to the decline in numbers, the population has undergone a contraction within its range (Rugh et al., 2010; Shelden et al., 2015). Summer surveys during the 1970's found beluga whales distributed throughout much of the upper inlet and into the lower inlet around Kalgin Island. Since the mid-1990's, from 96% to 100% of beluga whales now congregate in shallow areas near river mouths in the upper inlet during the summer months (Rugh et al., 2010; Shelden et al., 2015).
It is unknown if this contracted distribution is a result of changing habitat (Moore et al., 2000), prey concentration, or predator avoidance (Shelden et al., 2003), or can simply be explained as the contraction of a reduced population into a small number of preferred habitat areas (Goetz et al., 2007, 2012), such as in a "basin" model (MacCall, 1990). While the recent trends in abundance and range are well documented, little is known about other mechanisms influencing the recovery of this beluga whale population.
The CIBW population may be less resilient to natural perturbations or anthropogenic impacts because of its small size and isolation. With such a small population in a relatively restricted area, a substantial portion of the population could rapidly be exposed to events such as infectious disease outbreaks, volcanic eruptions, fish run failures, and toxic spills (Moore et al., 2000; Vos and Shelden, 2005).
Additionally, the grouping behavior of the population could potentially magnify exposure to even very localized anthropogenic and environmental hazards. During the summer months, Cook Inlet beluga whales tend to be found in 2-10 groups of a few individuals to over 200 whales in a single group (Hobbs et al., 2015a; Shelden et al., 2015) in areas such as Knik Arm, the Susitna Delta, and Chickaloon Bay-Turnagain Arm (Fig. 1). If a group of 200 whales were exposed to a toxic spill, this would represent over half of the current population.
The population's small size and grouping behavior also mean that a relatively large percentage of the population could be involved in a single mass stranding (Vos and Shelden, 2005), or entrapped in ice similar to events documented in other beluga whale populations (Siegstad and Heide-Jorgensen, 1994; Heide-Jorgensen et al., 2002), magnifying the effect on the population's recovery. Given the population's restricted range, particularly in summer, declines in local fish stocks, such as Pacific salmon, Oncorhynchus spp., runs in the rivers in Cook Inlet (Eggers and Irvine, 2007; Dischner (2)), could cause nutritional limitation in the population. Predation may also be a factor in the recovery of the population; between 1999 and 2014 a total of 10 beach-cast and floating carcasses reported to the NMFS Alaska Regional Office (NMFS3) were determined to be related to killer whale, Orcinus orca, predation. Three of the carcasses were lactating females, so an additional 3 calves were thought to have died though carcasses were not found (NMFS (3)).
Therefore, potential factors for the delay in recovery include reduced fecundity because the mature female segment of the population is depleted, reduced fecundity or survival due to reduced prey, predation by killer whales, and risks associated with a contracting range and grouping behavior of the whales. To examine these issues, we conducted a population viability analysis (PVA) by developing a detailed population model that included age and sex structure as well as small population effects. Small population effects taken into account included demographic stochasticity, hunt mortality, density-dependent and density-independent effects, constant mortality effects (e.g., predation), unusual mortality events (e.g., catastrophes), and risks associated with the grouping behavior of these whales.
The population model implicitly considers the time lags inherent in long-lived populations where sexual maturity does not occur for many years (Litzky, 2001), and it was projected into the future to examine extinction risk under different risk factor scenarios for CIBWs. The rather detailed population model and Bayesian framework were programmed in FORTRAN and used existing data from the Cook Inlet population and similar beluga populations.
The parameters of the model were estimated by fitting the population model to the abundance time series and hunt data. Variations of the model, termed model scenarios, were developed to investigate a number of hypotheses for the delay in recovery and these hypotheses were statistically compared using Bayesian model selection, which provides the probability of each model scenario (or hypothesis) given the data, through the use of the Bayes factor. We also examine future recovery and extinction risk under each hypothesis.
We considered the following hypotheses (H) explicitly:
H1) The population has not yet begun to recover but it will under the status quo--there has been a delay (time lag) in recovery due to a depletion of mature females, and the population will begin to increase once it rebuilds the mature female segment of the population.
H2) The subsistence hunt was not the sole cause of the decline observed since 1994 and the population has not begun to recover because an as yet unidentified population wide stressor (e.g., nutritional stress, disease, contaminants, noise, loss of key habitat) has caused a decrease in fecundity and/or survival, resulting in a negative growth rate. The population will not begin to recover in the future without a change in the fecundity or survival of the population.
H3) The population has not begun to recover because the reduction in population size from hunting has led to it falling into a "predator pit"--a numerically constant level of predation was sustained by the population when it was >1,000 animals, but at the current population level, predation prevents recovery.
H4) The population has not begun to recover because the level of predation from killer whales coincidentally increased around the time of the population decline (e.g., due to a prey shift by the killer whale population).
H5) The population has not begun to recover because the reduction in population size from hunting and subsequent retraction in range has led to a greater proportion of the population being vulnerable to "catastrophic" mortality events, such as disease outbreaks, oil or toxics spills, volcanic eruptions, or failure of several Pacific salmon runs in one year.
H6) The population has not begun to recover because the reduction in population size from hunting and subsequent retraction in range has led to fewer social groups in the population, such that a greater proportion of the population is vulnerable to mortality events that affect all or a large fraction of a single social group.
H7) The population has not begun to recover because of a combination of killer whale predation, per capita effects, and increased catastrophic or group mortality following population decline.
H8) The population has not begun to recover because the carrying capacity of Cook Inlet for belugas has declined during the population decline.
We use the Bayes factor to compare model scenarios and test hypotheses. The posterior sample of each model scenario was projected forward 100 years to estimate the probabilities of population recovery, increase, decline, and extinction under each hypothesis. The results of the model comparison were then used to identify the best model scenarios for assessing population viability.
The PVA was conducted by fitting a population model for the CIBW population to available abundance data using Bayesian statistical methods. In a Bayesian analysis, prior distributions for the model parameters are combined with a likelihood function for the data to give posterior probability distributions for the parameters, from which inference is based (Gelman et al., 1995; Ellison, 1996; Wade, 2000).
Prior distributions for the model parameters were specified using information from the CIBW population or, if necessary, from other beluga populations. The population model was initiated in 1979 (when the population was thought to be near carrying capacity) which allowed 15 years for the age distribution to accommodate to the presumed hunting removals before the parameters of the model were estimated by fitting the model to a time series of abundance estimates for the years 1994-2014 (Table 1). The population model was additionally projected 100 years into the future from 2014.
The analysis can be viewed as having two stages, a parameter estimation stage and a future projection stage, but it was conducted as a single integrated analysis. This ensures that the PVA projections are fully consistent with the available data (Wade, 2002). Additionally, Bayesian model selection methods were used to compare how well different model scenarios fit the data. The interpretation of the results focuses on the estimated probabilities of extinction or recovery and realized rates of growth or decline.
The methods below are presented in three stages which proceed from the general to the specific. First the population model is described and the various components of the survival and reproductive models are developed mathematically. This is followed by a section which presents the parameters and prior distributions specific to the modeling of the CIBW population. Finally, the methods for developing the posterior distributions, conducting the inference, and testing of hypotheses are described. Comparisons between the model development and the parameter values and prior distributions used in the analysis are presented in the Parameter Estimation section. The relationship between the hypotheses model scenarios and the model components included in each scenario are presented in the Model Scenario Comparison and Selection section.
An age and sex structured population model was developed using life history and population parameters from Cook Inlet and other beluga whale populations (Table 2).
Age and Sex Structure
Age-classes included each year up to maturity to account for the time lag from birth to sexual maturity (Litzky, 2001), and the preference of Native subsistence hunters for adult animals. Females and males were modeled separately to incorporate sex structure into the model and allow for unequal hunt of each. For both males and females, all adults greater than the age of sexual maturity (8+ growth layer groups (GLGs) Table 2 were lumped together for convenience.
Demographic and Environmental Stochasticity
The numbers of individuals in each age- and sex-class were tracked as integers. Births and deaths were modeled as discrete integer events using binomial distributions with the expected birth rate and survival rate, respectively. Births were modeled as a binomial draw using the birth rate and the total number of mature females. Survival from one age and sex class to the next was modeled as a binomial draw using the survival rate and the number of individuals in that age and sex class. The use of a binomial distribution in the population model incorporates demographic stochasticity, the random variations in the number of individuals that happen to die or reproduce in a given year even when the expected rate remains constant (Begon et al., 1996:927). Recorded takes from the Native subsistence hunt were partitioned among adult males, adult females, and the older age classes of immature animals of both sexes. Recorded takes were directly subtracted from each class (Table 1), after modifications to account for uncertainty, such as allocation to sex and struck and lost whales (explained below).
No environmental time series and mechanism has been identified as impacting survival or fecundity of the CIBW population, however demographic stochasticity does not explain the variation observed in the calf index (Hobbs et al., 2015b) or the annual numbers of deaths recorded (NMFS (3)). To account for this observed variation we included a correlated random variation scaled to the variation in the two data sets.
The CIBW population was projected annually for females (f) and males (m) as:
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (1)
Where a is age and B(n, p) is a binomial distribution with n trials and probability p, the variables s and b represent survival and birth probabilities, respectively, which are elaborated below.
Density dependence was included in the birth rate in the form of the commonly-used generalized logistic function (Breiwick et al., 1984).
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (2)
K represents the population size at carrying capacity or, in other words, the population size that the population would reach if there is no reduction in birth rate or survival rate except due to density dependence and there is no environmental variation. The parameter z controls the shape of the density-dependence, where values of z > 1 mean that most of the density-dependent response occurs close to K. [[sigma].sub.b] is a scaling factor to the effect of environmental variation on birthrate relative to the effect on survival, [[epsilon].sub.t] is a stationary, correlated, random environmental deviation with mean = 0, variance = [[sigma].sub.2], and correlation = [rho] (Morris and Doak, 2002:139) and [[omega].sub.t] is a normal random deviate with mean = 0 and variance = [[sigma].sup.2], the same series was used for both birth rate and survival. [B.sub.c] is a constant fractional reduction in per capita birth rate such as would occur as a result of chronic under nutrition or high contaminant loads; [s.sub.0] is the survival rate for the unaffected population when it is small and [S.sub.c] is the survival rate reduction from the per capita effect (described below). [B.sub.c] was set to 1 except in scenarios with the per capita effect where [B.sub.c] was the inverse of the proportional increase in mortality so that this effect was split approximately equally between mortality and reproduction. Note that without environmental variation, when the population is at carrying capacity so that [N.sub.t] = K, then [b.sub.t] = [b.sub.k] [B.sub.c], and when the population is small so that [N.sub.t] is close to zero, then [b.sub.t] = [b.sub.0] [B.sub.c].
The variables [s.sub.0,t], [s.sub.a,t] [s.sub.f,mat,t], and [s.sub.m,mat,t] are annual survival rates by age and sex. They are the products of two components: 1) [S.sub.t], a population wide survival rate that is independent of age and sex, and determined by density dependence; and 2) [S.sub.h,age-sex,t], an age and sex dependent survival rate determined by hunting activity (h) on adult males, adult females, and immature (j) whales.
[s.sub.0,t] = [S.sub.t] [s.sub.f,mat,t], which accounts for a calf's dependence on its mother; (3a)
[s.sub.a,t] = [S.sub.t], for a = 1 to 3; (3b)
[s.sub.a,t] = [S.sub.t] [S.sub.h,j,t], for a = 4 to mat-1; (3 c)
[s.sub.f,mat,t] = [S.sub.t] [S.sub.h,f,t], for adult females; and (3d)
[s.sub.m,mat,t] = [S.sub.t] [S.sub.h,m,t] for adult males (3e)
The hunt mortality was modeled as the sum of the recorded landed whales with different rates for adult males, adult females, and immature whales. Additional mortality was added to account for whales that were struck and lost (injured or killed during a hunt but not retrieved and landed), specified by year. The total hunt mortality ([H.sub.t] = landings + struck and lost) was allocated to age and sex class using a binomial distribution and parameters controlling the expected age and sex bias in the hunt:
[H.sub.j,t] = B([H.sub.t], Pr(Hunt Immature)) (4a)
[H.sub.m,t] = B([H.sub.t] - [H.sub.j,t], Pr(Hunt Male)) (4b)
[H.sub.f,t] = [H.sub.t] - [H.sub.j,t] - [H.sub.m,t] (4c)
Where t is the year, [H.sub.j,t] is the number of immature whales killed, [H.sub.m,t] is the number of mature male whales killed, and [H.sub.f,t] is the number of mature female whales killed. The Pr(Hunt) were chosen from a Beta distribution prior for each case of the model, as described below. The numbers of whales killed were translated into survival rates from the hunt as:
[NV.sub.t] = [[summation].sup.mat-1.sub.a=4] [f.sub.a,t] + [m.sub.a,t] (5a)
[S.sub.h,j,t] = 1 - [H.sub.j,t]/[NV.sub.t] (5b)
[S.sub.h,m,t] = 1 - [H.sub.m,t]/[m.sub.mat,t] (5c)
[S.sub.h,f,t] = 1 - [H.sub.f,t]/[f.sub.mat,t] (5d)
representing survival of immature whales, mature males, and mature females, respectively. [NV.sub.t] is the total number of immature whales vulnerable to hunting. The hunt survival parameters are formulated as survival rates so that if there is no hunt mortality the rates all equal 1.
Population Level Survival Rate
Density dependence was included in the population level component of survival. Variations of the survival model were also specified to allow for (depending upon the model scenario) changes in survival from other factors. The other factors include 1) a constant decline in the survival rate (a "per capita impact") such as might occur as a result of chronic under nutrition or high contaminant loads, 2) a constant numeric decline in survival (e.g., one additional whale per year) such as could occur from density-independent predation by killer whales, 3) occasional unusual mortality events ("catastrophes") that affect a percentage of the entire population, and 4) occasional mortality events ("group mortality") that affect an entire social group, which is a form of an Allee effect (Allee et al., 1949) as this mortality has a greater effect as the population becomes smaller and has fewer social groups. The population level survival rate is:
[S.sub.t] = [[s.sub.0] - ([s.sub.0] - [s.sub.k])[([N.sub.t]/K).sup.Z] + [[epsilon].sub.t]] [S.sub.c][S.sub.p,t][S.sub.e,t][S.sub.g,t] (6)
Where [S.sub.c] is the constant per capita impact, [S.sub.p,t] is the predation effect, [S.sub.e,t] is the catastrophe effect, and [S.sub.g,t] is the group mortality effect (details on calculations of these factors are found below). [[epsilon].sub.t] is the same environmental series used in the birth rate. Note that these additional mortality factors are formulated as multiplicative adjustments to [S.sub.t] that act independently. These specific survival models can then be employed together or separately.
Killer Whale Predation Mortality
Predation from killer whales was modeled as a density-independent survival factor, meaning the expected number of additional deaths will remain the same regardless of the CIBW population size. This assumes that killer whales prey on belugas irrespective of the size of the population (i.e., there is no numerical or functional response by the predator). Killer whale predation mortality was modeled in this way because predation on CIBWs would represent such a small fraction of the annual diet of mammal-eating killer whales in the region that there would be no predation response to the density of CIBWs.
This means that predation mortality was assumed to be determined by factors such as the frequency of killer whale visits to the upper inlet. By modeling predation this way, the same number of predation deaths was a larger percentage of the population size when the population was smaller, and increased the risk to the population. Eventually, a population may decline to the point where predation cannot be sustained and thus prevents its recovery, this has been termed a "predator pit."
Predation mortality from killer whales was specifically modeled as
[S.sub.p,t] = [[N.sub.t] - [C.sub.s][C.sub.p]]/[N.sub.t] (7)
where [C.sub.s] is estimated as the average number of observed deaths from killer whales during the time period 1999 to 2014, and [C.sub.p] represents a scaler to allow for decreasing or increasing the level of killer whale predation in the model scenario (e.g., a value of 2 means that twice as many predation deaths occurred than were observed from 1999 to 2014). Note that [S.sub.p,t] was also constrained to be no less than 0 for the case where the constant number of expected killer whale predations was greater than the population size.
Catastrophes (Unusual Mortality Events)
Occasional unusual mortality events, including but not limited to disease outbreaks, mass strandings, volcanic activity, toxic spills, and failure of Pacific salmon runs, could potentially impact a substantial portion of the CIBW population and were, therefore, important to consider. Catastrophes were modeled as mortality events where an additional specified fraction of the population died in a given year:
[S.sub.e] = 1 - [M.sub.e] B(1, [P.sub.Me]) (8)
where [M.sub.e] is the probability of mortality during an unusual mortality event (e.g., where a value of 0.1 means that 10% of the population was expected to die, in addition to the density-dependent mortality that takes place); and [P.sub.Me] is the binomial probability of an unusual mortality event occurring in a given year. The catastrophe survival factor models random events that affect the entire population. If [P.sub.Me] equals 0, no catastrophic events occur because [S.sub.e] = 1.
Given that a large percentage of the CIBW population may aggregate in a single behavioral group at one time, another way to model an unusual mortality event would be to have it affect a single social group. Events that could affect an entire social group at the same time include entrapment in ice, becoming stranded in a shallow area at low tide, a local toxic spill, or a stranding related to an acoustic event.
A group mortality event can be viewed as a different form of a catastrophic mortality event, but where the size of the social group and the mortality rate within the group determines the number of whales that die, rather than being specified as a fixed proportion of the population. Therefore, catastrophes and group mortality were not included together in any model scenario but were instead viewed as alternative ways of modeling unusual mortality events.
To model this type of mortality affecting all or part of a single social group, a group survival factor is modeled as:
[S.sub.g,t] = 1 - ([M.sub.g][G.sub.t]/[N.sub.t]) B(1, [P.sub.g]) (9)
where [P.sub.g] is the probability of an event occurring in a given year that would lead to a group mortality; [M.sub.g] is the probability of mortality for individuals in the group affected by a group mortality event; and [G.sub.t] is the size of the affected group drawn for each event from the observed distribution of group sizes (truncated at [N.sub.t]). Whether or not an event occurs in a given year is determined by an annual draw from a binomial distribution with probability of occurrence in a given year of [P.sub.g] such that either one or no event occurs in that year.
For this factor, groups were collections of belugas in a locale and may include several social groups in proximity to each other or possibly one large group. If the group size drawn was equal to the population size, then [G.sub.t] was replaced with [N.sub.t] and the Group Mortality model was the same as the Catastrophe model. The risk of an event remained the same but the group size changed. The individual risk in this model at a particular population size is proportional to the ratio of the average group size in the population to the population size, which increases as abundance declines because there are fewer groups, so the chance of being in the group affected increases.
This model implies that only one group is affected by any event so, for example, a localized toxic spill at the Port of Anchorage would affect a group in Knik Arm but not impact groups in other areas such as the Susitna Delta or Chickaloon Bay. Similarly, this models a group stranding event that affects only the group trapped by a falling tide. As the population increases and divides into more groups, it becomes less vulnerable to these sorts of events. If [P.sub.g] = 0.0, no group mortality event occurs in any year.
Starting Population Size and Initial Age Distribution
The model parameters were estimated by fitting the population model to abundance data starting in 1994, when NMFS began an annual series of systematic counts with applied correction factors for the CIBW population estimates (Hobbs et al., 2015a). The 1979 abundance estimate was from a single day survey to which Calkins (1) applied a correction factor based on radio-tag data from the Bristol Bay population (Frost et al., 1985). This estimate was used as the basis for the prior distribution for population size in the start year of 1979, the population was then projected forward to 1994 to allow 15 years for the hunt to affect the sex and age distribution of the population. The population was initialized in 1979, with a stable age distribution calculated from the life-history parameters with the survival rate associated with the population size in 1979, and with a population growth rate of 0.0 ([lambda] set to 1.0, see below). Age and sex classes were filled by randomly sampling as a multinomial distribution from the stable age distribution until the initial population size was reached.
Values for model parameters were taken from data available on beluga populations (Table 2). For parameters for which little data were available or inference was intended, prior distributions were devised as described below. Some parameters were fixed at a single value. A "healthy" population model scenario was specified with a growth rate that was considered to be typical of increasing cetacean populations and no additional mortality from per capita effects, predation, catastrophes, or group mortality. Other model scenarios were specified that included additional mortality as a fixed number of deaths or fixed rate or drawn from a prior distribution of deaths or rates from those factors. The fit of the different model scenarios to the data were compared using Bayesian model selection methods (detailed below).
Annual Growth Rate ([lambda])
The population rate of increase ([lambda]) is not among the parameters of the population model above; however, for given values of the life history parameters (e.g., survival, birth rate, age of sexual maturity) in the population model, [lambda] is determined, because they are functionally related (Euler, 1760; Lotka, 1907). We use a simplified version of the population model as follows. Treating s and b, as constant parameters (rather than probabilities), we have a deterministic projection of the expected values of the abundance and individual age and sex classes. Considering only females, we have a recursion model in expected births by year:
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (10)
This is a discrete form of the Lotka renewal equation (Lotka, 1907; Goodman, 1982), which has a solution for a constant rate of increase by replacing [f.sub.0,t] with [[lambda].sup.t]. Dividing through by [[lambda].sup.t] yields:
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (11)
Note that if the rate of increase and two of the three life history parameters are known, the third life history parameter is determined. Consequently, only three of the four parameters can be assigned values and the fourth will then be determined by the other three. For these analyses, we chose values for [lambda] because the inference included a healthy population model with fixed growth rate range. In the Bayesian analysis, we chose to put an uninformative prior distribution on [lambda] (i.e., a uniform distribution) because we were interested in the probability of the future increase or decline of this population. Furthermore, we assigned prior distributions to age at sexual maturity and survival rate; therefore, birth rate was a derived parameter.
The quantities of interest are the maximum rates of increase ([[lambda].sub.0]) of the population, which under compensatory density-dependence, is the annual growth rate near zero population size. For a population that is doing well and has the capacity to increase, [[lambda].sub.0] should be greater than 1.0. However, if the population is in a state of decline, [[lambda].sub.0] may well be less than 1.0. To allow for both possibilities, the prior distribution for the maximum annual growth rate, [[lambda].sub.0], was specified to be a uniform distribution between 0.94 and 1.031 in order to provide a broad, uninformative prior distribution for this parameter.
It is thought that the maximum annual increase for an odontocete cetacean with a life history such as a beluga whale is unlikely to be much greater than about 1.04 or 4% per year (Reilly and Barlow, 1986; Wade, 1998, 2009). The BB population has increased at an estimated rate of 4.8% per year (95% CL: 2.1%-7.5%) (Lowry et al., 2008). Therefore, the value of 1.06 (a 6% increase per year) could be considered a reasonable upper bound for this parameter. However, carcass counts and calving rate data, described in more detail below, provide maximum possible survival and birth rates which because of the functional relationship of these parameters effectively reduces the upper limit for [[lambda].sub.0] to 1.031.
Log-linear regression of the CIBW abundance estimates for 1999-2014 indicate an annual rate of decline of -1.3% (SE = 0.008) (Shelden et al. (4)); consequently, the lower bound of 0.94 (representing a 6% decline per year or the slope of the recent trend minus more than five times the SE of the slope) was set to a value thought to be low enough to capture any possible outcome from the analysis. The results for each model scenario were examined to ensure that the data did not support any greater rate of increase or decline. Model scenarios were rerun with a broader prior distribution where support was indicated.
If the population is healthy and growing, such as the BB population, we would expect the value of [[lambda].sub.0] to fall between 1.02 and 1.031, (i.e., annual growth of 2% to 3.1%) and by definition, the value of [lambda] at carrying capacity, what we call [[lambda].sub.k], for a healthy population would be 1.0 (i.e., no annual growth or decline), causing the population size to stabilize and reach equilibrium near the value K. Each of the alternative scenarios was modeled as a modification to this healthy population model. To account for several different survival or fecundity models, given that [[lambda].sub.0] was allowed to be less than 1.0 in some scenarios, we devise a healthy population growth rate for a small population, [[lambda].sub.H0], with a prior distribution between 1.02 or [[lambda].sub.0] (whichever was larger) and 1.031, which spanned the range of what would be the expected maximum rate of increase for a beluga whale population.
The model scenario with no additional mortality (healthy), specified a prior distribution for [[lambda].sub.0] of a uniform distribution between 1.02 and 1.031. In this case, [[lambda].sub.0] and [[lambda].sub.H0] are equal, with the assumption that the population will increase if below K. If a scenario includes reduced survival and/or fecundity [lambda] is less than [[lambda].sub.H0] and the difference of [[lambda].sub.H0] - [[lambda].sub.0] determined the magnitude of the reduction of fecundity and/or survival under the alternative scenario. The values of [[lambda].sub.H0] and [[lambda].sub.K] = 1 at N = K determine the density-dependent drop in the population growth rate. For example, if [[lambda].sub.H0] is 1.03 and [[lambda].sub.0] is 0.99, density dependent change in the growth rate is -0.03, whereas, the change resulting from reduced survival and/or fecundity is -0.04 at both N = 0 and N = K.
A uniform prior distribution was set for [s.sub.0], the maximum population level survival rate, when the population is near zero and where (under the standard assumption of density dependence) the population is assumed to be healthy and growing at a maximal rate. A range (0.962-0.975) was specified for this distribution to make it an uninformative prior distribution, while avoiding values that were not possible for the prior range of [[lambda].sub.H0] and the range of birth rates.
The upper value for the prior distribution was calculated by estimating the maximum survival rate that could be achieved by the CIBW population. A minimum annual mortality rate was estimated from annual summaries of beach-cast and floating CIBW carcasses reported to the NMFS Alaska Regional Office (e.g., Vos and Shelden, 2005; NMFS (3)). From 1999 through 2014, a total of 148 carcasses were found, of which 10 were attributed to killer whale predation and, as three of the deaths were lactating females, an additional 3 calves were thought to have died though carcasses were not found (NMFS (3)). The remaining 138 carcasses result in an average of 8.6 documented deaths per year over the 16-year timespan (Allen and Angliss, 2013), or 2.5% per year (SE = 0.3%) from a population size that has averaged 346 animals during those years. The lower value for the prior distribution was set to the minimum value that would allow a growth rate of 1.02 with a birth rate of 0.20 per year to encompass the range of possible values for the healthy population growth rate in Cook Inlet.
The value of [S.sub.K] was set to [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] to partition the density dependent effect approximately evenly between the survival and birth rate. The priors for the parameters of the environmental variation were determined from the variation in the annual counts of beach-cast carcasses. The value for [sigma] was drawn from uniform [0.005, 0.01], the square root of the variance of the annual observed mortality rates (carcass count/average population size) less the variance of the demographic stochasticity. The value for [rho] was chosen from a uniform [0.5, 0.8].
Age of Maturity
The prior distribution for the age of maturity ([a.sub.mat]), or the age at first possible birth was set at 9 GLGs based on an average age of first pregnancy of 8.25 GLGs which would result in a first birth at age 9 GLGs (Table 2).
With the population rate of increase, survival rate, and age of maturity specified we solve Equation 11 for birth rate (b) as a function of the other parameters to get Equation 12. Equation 11 was solved with extreme values of b and s to identify the limits for survival and intrinsic rate of increase that will allow b to fall into the specified ranges (Brandon and Wade, 2006). Equation 12 is then solved for b given s, [a.sub.mat], and [lambda]:
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (12)
This equation was used to determine the annual probability of giving birth when the population was small ([b.sub.0]) or large ([b.sub.K]), using the values for [[lambda].sub.H0] and [s.sub.0] or [[lambda].sub.K] = 1 and [S.sub.K], respectively. The values for [s.sub.0] and [S.sub.K] were constrained so that [b.sub.0] and [b.sub.K] fell into the biologically reasonable range for these variables of [0.05, 0.20] and [0.0, [b.sub.0]], respectively. Ranges were derived from the data on pregnancy rates and birth interval from wild beluga populations (Table 2) and a calving rate study conducted in Cook Inlet (Hobbs et al., 2015b), with 0.05 set lower than the lowest value of 0.13 found in the literature to allow for the possibility of poorer fecundity in Cook Inlet and the value 0.20 set at nearly twice the value estimated in Cook Inlet to allow for the possibility that the time period in which the calving rate was estimated was unusually low and to also allow for potential bias in the index itself. The value for [[sigma].sub.b] is set to 33 so that variation in annual birth rate had an equivalent effect on the variation in annual growth as on the variation in survival. The same correlated random series was used for both survival and birth rates.
The prior distribution for carrying capacity (K), the population size that the population would reach if there is no reduction in birth rate or survival rate except due to density dependence and there is no environmental variation, was a uniform distribution U[811, 2056] based on a log-normal 95% confidence interval for an estimated 1,293 beluga whales from an aerial survey conducted in August 1979 by the ADFG (Calkins (1); Hobbs and Shelden (5)). No error was calculated for this estimate; thus, we have assumed that it was no more accurate than the 1994 estimate, and we have applied the CV from the 1994 aerial survey estimate.
While this is the best available estimate of population size before the 1990's, it potentially represents a population that was depleted to an unknown degree by poorly documented and undocumented removals in earlier years of the 20th century (Mahoney and Shelden, 2000). The prior distribution on K represents the uncertainty in the estimate from 1979. For hypothesis H8 we also consider a change in K over time. For these scenarios, we reduce K by a constant amount each year between 1979 and 1999, the values for [K.sub.1999] in 1999 are drawn from a uniform [100, 800]. Carrying capacity continues at a constant level after 1999.
The parameter z determines the maximum net production level (MNPL; the population size where growth in numbers is greatest) as a fraction of K and determines the shape of the growth curve of the population (Taylor and DeMaster, 1993). MNPL has not been determined for Cook Inlet or any other beluga whale population, consequently we adopted a range considered reasonable for cetaceans of 50-80% (MNPL 650 to 1,040 for K = 1,300). MNPL/K is drawn from U[0.5, 0.8], then the value of z is determined by solving 1/[z+1] = [(MNPL/K).sup.Z] iteratively for z. This modeling effort was focused on the behavior of the population at sizes at or below the current population size, which is well below these values of K and MNPL. Therefore, the specific values of K and MNPL will have little or minimal influence on the results.
The prior distribution for the number of whales killed in the hunt was handled differently for two time periods. The total number of whales killed each year were poorly documented before 1994 (Mahoney and Shelden, 2000). Reported landings averaged 10 per year for the years 1987-89 and 1991-93 when a partial survey of Native hunters was conducted (Stanek, 1994). But only in 1993 was there an effort to estimate the complete removals (26 landed beluga whales). No independent study of struck and lost rates was conducted during the Cook Inlet hunt.
Struck and lost numbers reported by hunters during 1979-93 averaged 25% of the reported landings, resulting in an average of 12.5 whales reported killed in those years. In 1993 the estimated kill was 32.5 whales. Therefore, to allow for the uncertainty during this time period, the prior distribution for the number of animals killed each year for the years 1979-93 ([H.sub.79-93]) was specified as a uniform distribution from 10 to 40 whales to span the range of possible take levels.
NMFS first documented the entire take (landings, and struck and lost) by the subsistence hunt starting in 1994 (Mahoney and Shelden, 2000), so the data for years 1994-98 are believed to be accurate. Hunting dramatically declined starting in 1999 with the moratorium on the hunt and with subsequent harvest management plans, and recorded landings and struck and lost whales from the hunt are thought to be accurate for the years 1999-2014 (NMFS (3)) (Table 1). Therefore, the prior distributions for total kill (landings plus struck and lost) in the years 1994-2014 were fixed at the reported values. In years when more than one value was reported, or a range was reported, the total killed in the hunt was drawn from a uniform distribution with limits at the smallest and largest value (see Table 1, i.e., 1995 U[68,74]; 1996 U[98,147]; 1997 U[65,75]).
The youngest animal documented in the hunt was 4 GLGs (Mahoney and Shelden, 2000), so the first immature age vulnerable to the hunt was set to age 4. Although the hunters are thought to have a preference for taking large white adults, some gray immature beluga were also taken in the hunt. A review of hunted whales listed in Table 2 in Mahoney and Shelden (2000:131) shows 33 individuals for which age or length was known and sex determined, 7 of which were immature (i.e., smaller than 340 cm, not pregnant or lactating, and/or younger than 9, with the youngest at 4), of these immature whales, four were females.
While 33 is a relatively small sample of the more than 250 whales killed in the hunt between 1994 and 1998, it does provide a basis for an informed prior for the fraction of the hunt that was immature animals. The standard Beta distribution, Beta([alpha], [beta]), is considered the likelihood distribution for the probability of success estimated from binomial data with [alpha]-1, successes and [beta]-1, failures (Johnson and Kotz, 1970). The distribution for the probability that an individual taken in the hunt is immature is then Beta(8, 27), with the parameters derived from the values 7 immature captures (successes) and 26 mature animals (failures) recorded in the hunt data. Therefore, the informed prior distribution for the probability that an animal taken in the hunt was an immature, Pr(Hunt Immature), was drawn from Beta(8, 27), truncated at the 2.5th and 97.5th percentiles.
For mature whales, the hunt was biased towards males. Of the 26 mature whales for which sex was determined, 16 were males and 10 were females. The male bias may be due to the fact that when a calf was seen closely associated with a whale (presumably the mother), the hunters would break off pursuit, creating a bias toward taking males (Huntington, 2000). The Beta distribution was also applied; therefore, the informed prior distribution for sex-bias of mature belugas in the hunt, Pr(Hunt Male), the probability that a mature animal taken in the hunt was a male, was drawn from a standard Beta distribution, Beta(17, 11), truncated at the 2.5th and 97.5th percentiles.
Modifications to Birth and Survival Rates
To account for the observed growth rate of the population, we compared the effects of four different survival factors, each of which is formulated as a modification of the survival rate, [S.sub.c], [S.sub.p,t], [S.sub.e,t], and [S.sub.g,t], and one modification of the birth rate, [B.sub.c]. Of these, only [S.sub.p,t], can be parameterized from data for the Cook Inlet beluga population and is either equal to 1.0 or the value determined by the parameters. For the others, we use the [[lambda].sub.0], ([S.sub.x] [S.sub.p] [s.sub.0]) and [b.sub.0] in equation (11), solving iteratively for [S.sub.x], where [S.sub.x] is the product of [S.sub.c], E([S.sub.e,t]), and E([S.sub.g,t]). Note that these modifications do not include density dependence so it was unnecessary to solve using N = K, [s.sub.k] and [b.sub.k]. E(X) indicates the expected value of a stochastic process X. [S.sub.p] is the value of [S.sub.p,t] when the population size is 350 belugas.
When more than one modification was used in a model scenario, the value for [S.sub.x] is partitioned into survival effects which multiplied together to equal [S.sub.x]. This was done using random partitions between zero and one that sum to one as exponents of [S.sub.x]. For example, if 3 partitions are needed two random numbers between 0 and 1 are drawn and then sorted into order with 0 and 1. The first partition is equal to the value of the lower number, the second is the difference between the two random numbers, and the third is the difference between the larger random number and 1, so that product of the three partial powers of [S.sub.x] is [S.sub.x].
If the model scenario includes killer whale predation, which has fixed parameters, this has been accounted for and [S.sub.x] is partitioned to account for the other survival effects if it is less than 1. In a model scenario with the per capita effect, both survival rate and birth were modified, so first the partitioning described above is completed then [S.sub.c] is raised to the power 0.5 and then [B.sub.c] = [1 - [s.sub.0]]/[1 - [s.sub.0][S.sub.c]] so the decrease in birth rate is proportional to the increase in mortality.
Per Capita Survival and Fecundity Effects
When these survival and fecundity effects are included in a model scenario they are used together. The values for [B.sub.c] and [S.sub.c] are based on the necessary partitioning of the value of [S.sub.x]. Both [B.sub.c] and [S.sub.c] are constrained to be greater than zero and less than or equal to 1. These are intended to model effects that impact the population on a per capita basis such as reduced frequency of foraging events resulting in chronic under nutrition, or significant contaminant loads resulting in poorer survival and fecundity.
Killer Whale Predation Mortality
Killer whale predation was modeled as a constant expected number of individuals. [C.sub.s] is estimated as the average number of observed deaths from killer whales during the time period 1999-2014. The data available on killer whale predation in Cook Inlet are based on recovered carcasses, and determining the actual number of deaths depends on knowing the discovery rate of carcasses. To get around this problem, [C.sub.s] is calculated from the expected number of total deaths (the total mortality rate in the population model times the population size), and from the proportion of carcasses determined to be the result of killer whale predation.
During the period 1985-2002, known beluga whale deaths resulting from killer whales averaged about one per year (Shelden et al., 2003). More complete records are available for the period between 1999 and 2014, when 148 beluga carcasses were discovered and of these 10 were determined to be the result of killer whale predation (NMFS (3)), or 6.8% of all observed carcasses. Therefore, if we assume that 6.8% of all deaths between 1999 and 2014 occurred from killer whale predation, we can estimate the number of predation deaths as [C.sub.s] = (0.068) (1-[s.sub.346]) (346) where, [s.sub.346] is the density dependent survival rate at N = 346, which was the average population size during the period 1999-2014.
These numbers may be low relative to actual deaths because belugas that are nearly or entirely consumed will not be recovered; indeed, three of the carcasses were lactating females but no calf carcasses were recovered suggesting that the calves were consumed. The parameter [C.sub.p] represents a scaler to allow for decreasing or increasing the level of killer whale predation in the population model (e.g., a value of 2 means that twice as many predation deaths occurred than were observed from 1999 to 2014). Note that when [C.sub.p] = 0, there is no mortality from predation as [S.sub.p,t] = 1. The value of [C.sub.p] was set to either 0.0 (no predation), 1.0 (predation equal to the number estimated from carcass data collected from 1999 to 2014), or 2.0 (predation equal to twice the number estimated from carcass data such as associated calves or carcasses that were entirely consumed).
Catastrophes (Unusual Mortality Events)
No population-wide mortality event has been documented for Cook Inlet belugas, and the events that we are considering here, such a disease outbreaks, oil spills, volcanic eruptions, or failure of Pacific salmon runs, are infrequent. Examples are available for other marine mammal populations, such as recent disease outbreaks (6) in bottlenose dolphins, Tursiops truncatus, on the east coast of the United States and in ice associated seals, Phoca spp., in the Arctic waters of Alaska.
Oil spills, such as the Exxon Valdez in Prince William Sound, Alaska, and recently the Deep Water Horizon spill in the Gulf of Mexico, have resulted in deaths of a number of species of marine mammals (Matkin et al., 2008; Schwacke et al., 2014). Between 1976 and 2004, nine documented major oil spills occurred in Cook Inlet, ranging in volume from 5,700 gallons to 395,640 gallons. (7) In addition, oil from the Exxon Valdez spill entered Cook Inlet (8), resulting in an average of one spill per three years.
There are seven known volcanoes adjacent to Cook Inlet, four of which are sufficiently active to warrant monitoring by the Alaska Volcano Observatory (Miller et al., 1998). In addition, there are a number of volcanos south and west of Cook Inlet that could produce significant ash falls in Cook Inlet. While eruptions are infrequent, they can pose a substantial hazard when they occur.
Finally, anadromous species, such as Pacific salmon and eulachon, Thaleichthys pacificus, represent a concentrated feeding opportunity that the beluga population may depend on as an important component of their annual nutrition. While individual fish runs are known to vary from year to year, a coincidental failure of several runs in one year could have a substantial impact on beluga survival and reproduction. For example, other marine mammal populations, such as the southern resident killer whale population, have shown reduced survival and fecundity rates in years in which salmon returns to the Frazer River, B.C., Can., are low (Ford et al., 2010; Ward et al., 2009).
For catastrophic events, the expected or average annual survival rate E([S.sub.e]) has been set above for each population run either in the partitioning of additional survival effects or directly depending on the model scenario. The value of E([S.sub.e]) = (1-[M.sub.e] [P.sub.Me]). To determine the mortality rate per event, [M.sub.e], we use a broad prior with [M.sub.e] drawn from U[0.10, 0.50] for each population run. The frequency of events, [P.sub.Me], is then calculated as (1-[S.sub.e])/[M.sub.e] for that run. The values for [M.sub.e] and [P.sub.Me] remained fixed through each population run which resulted in an expected or average annual unusual mortality rate of (1-[S.sub.e]) for the population run.
The group size [G.sub.t] is drawn at random from the distribution of groups observed in Cook Inlet during aerial surveys from 1994 to 2014, truncated at the current population size (Fig. 2). Groups within 5 km of each other are treated as one. Note that for the current population size, the average observed group size represented about 20% of the total population (Fig. 3).
An example of one type of group mortality event that was observed for the time period was a stranding of 46 animals in 2003 in which 5 animals died (Vos and Shelden, 2005), representing a group mortality rate for the event of 5/46 = 0.11. Under the assumption that all group strandings were observed for the years 1994-2014, the probability of a group mortality event in this example would then be 1/21 = 0.05. This is but one example and not the only possibility for a group mortality event. The expected or average annual survival rate E([S.sub.g]) has been set above for each population run either in the partitioning of additional survival effects or directly depending on the model scenario. Thus, we specified a broad prior distribution for the group mortality rate, [M.sub.g], of U[(1-[S.sub.g])/0.20, 1.0], with the lower limit determined by the expected mortality rate if an event occurred every year and affected 20% of the population. A value for Mg was drawn for each population run then the value for [P.sub.g] is set to (1-[S.sub.g])/ ([M.sub.g] * 0.20) for that run so that the expected mortality rate for the population is (1-[S.sub.g]).
Population Size in 1979
To set up the initial population size and age structure in 1979 ([N.sub.1979]), [N.sub.1979] was drawn from a uniform distribution ranging from 800 to K belugas. The age structure for each sex was then drawn as a multinomial from [N.sub.1979] using the probabilities from a stable age and sex distribution with the density-dependent survival rate for [N.sub.1979].
The parameters of the population model were estimated by fitting the model to population abundance estimates available from 1994 to 2014 and the ratios of adult males and females and immature animals observed in the subsistence hunt. Aerial surveys have been conducted each June-July from 1994 to 2014 using essentially the same data collection methods through the entire time series (Hobbs et al., 2000a, b; Hobbs et al., 2015a). The subsistence hunt was subsampled as discussed above based on 7 immature belugas, 10 mature females, and 16 mature males. The likelihood function used was
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (13)
where [L.sub.i] is the relative likelihood of the ith population projection; T(X, DF Y) is the density of Student's t-distribution at X with Y degrees of freedom; [N.sub.t,i] is the population size of the ith projection in year t; [[bar.N].sub.t] and CV([[bar.N].sub.t]) are the estimated abundance (point estimate) and associated coefficient of variation in year t; [H.sub.j] is the sum of [H.sub.j,t], t = 1993 to 1998; [H.sub.f] is the sum of [H.sub.f,t], t = 1993 to 1998; and [H.sub.m] is the sum of [H.sub.m,t],t = 1993 to 1998. [bar.[epsilon]] is the average deviation of the environmental variation from 1994 to 2014 and SE([bar.[epsilon]]) is the standard error for [bar.[epsilon]] when [sigma] = 0.01.
The Student's t-distribution was chosen as the likelihood function for the abundance estimates by comparing the fit of several different likelihood functions to the distribution of bootstrap abundance estimates that results from analysis of these survey data. The Student's t-distribution with 10 degrees of freedom was chosen because it provided a better fit than a Student's t with 5 degrees of freedom, a gamma distribution, a log-normal distribution, or a normal distribution.
Each individual population projection was fully defined by the 18-parameter vector: [[lambda].sub.0], [[lambda].sub.H0], [[lambda].sub.k], [S.sub.0], [S.sub.K], [a.sub.mat], [N.sub.1979], K, [K.sub.1999], z, [H.sub.t], Pr(Hunt Immature), Pr(Hunt male), [C.sub.p], [P.sub.Me], [M.sub.e], [P.sub.g], and [M.sub.g] (Table 3). Note that two of the parameters ([a.sub.mat], [C.sub.p]) were set at fixed values for any particular model scenario. While this is a substantial number of parameters, the purpose was not to estimate posterior distributions for all of these parameters to draw inferences directly, but instead include sufficient detail in each model scenario to test the ability of particular model scenarios to reproduce the existing abundance data and determine the consequences for population projections. The time lags inherent in age-structured populations are a key element of this inference. For clarity, other derived parameters (meaning they are functions of the estimated parameters) that are referred to in the methods are listed in Table 4.
A Metropolis-Hastings Markov chain Monte Carlo (MCMC) algorithm (King et al., 2010) was used to generate samples of parameter vectors (and output quantities of interest) from the posterior distribution; these samples are then used to approximate the posterior distribution. This algorithm generates a random walk through the parameter space by selecting a test set of parameters from the neighborhood of the current parameter set, estimating the likelihood for the test parameter set, then applying a stochastic acceptance test to determine if the chain updated to the new parameters or remained at the current parameters.
All parameters were updated at once, drawn from distributions centered on each parameter value and spanning one tenth of the prior of each parameter. Where a distribution other than the uniform distribution is used as the prior, the cumulative probability was used as the uniform prior then transformed. Ranges of 0.01, 0.10, and 0.20 were also tested, but 0.05 was found to give the best performance, with higher acceptance rates similar to 0.01 but larger jumps.
For each model scenario, three chains of 10,000,000 trials were generated after an initial burn in of 6,000. Every 1,000th trial was retained to create three sets of 10,000 samples from the posterior distribution. When the algorithm repeated a parameter set, i.e. more than 1,000 trials occurred without a jump, the chain was forced to jump to the next highest likelihood point in the previous 1,500 trials and a burn in of 2,000 trials was completed before accepting the next point. This process insured that the algorithm made a thorough sampling of the posterior distribution. Inference was based on the combined sample of 30,000. We tested convergence of the three chains using the method of Brooks and Gelman (1998), which uses the ratio of the covariance matrix of the parameter sets from the individual chains to the covariance among their means to estimate the potential scale reduction factor (PSRF). Values of PSRF close to 1.0 indicate that each chain is a representative sample of the posterior distribution.
The population sizes for the projections from 1994 to 2164 were retained as output quantities of interest, covering the period for which we have data, 1994 to 2004, and projecting 150 years into the future. A population with one or zero individuals or only one sex was considered extinct, and its population size was set to 0 for the rest of the trajectory.
Model Scenario Comparison and Selection
Model scenarios were specified that included different prior distributions for [[lambda].sub.0], and different fixed values for the additional mortality factors (Table 5). A model scenario (Model A) was specified that assumed the population was healthy (or recovering) by setting the prior distribution for U[1.02, 1.031]. All other model scenarios are modifications to this healthy population scenario with [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] drawn from U[1.02, 1.031] and [[lambda].sub.0] drawn from U[0.94, 1.031] or specified by other parameters with additional mortality or reduced fecundity modeled, as described above. Model B is a scenario that included per capita mortality and fecundity factors, and had a prior distribution for [lambda], of U[0.94, 1.031]. Models C-E considered the effect of killer whale predation on an otherwise healthy population, in Model C, [C.sub.p] is set to 1.0, in Model D, [C.sub.p] is set to 2, and in Model E, [C.sub.p] changes from 1 to 2 in 1999 so that the increase in killer whale mortality coincides with the change in Native subsistence hunting practices.
The remaining models, F-O, had a prior distribution for [[lambda].sub.0] of U[0.94, 1.031], and where killer whale mortality is included, [C.sub.p] is set to 1. Model F accounts for the reduced growth rate with both killer whale predation and per capita effects. Models G-L account for reduced growth rates with stochastic mortality effects. Model G considers the catastrophic mortality by itself. Model H uses the group mortality alone. Model I is catastrophic mortality with killer whale mortality. Model J is group mortality with killer whale mortality. Model K partitions reduced growth among catastrophic mortality, killer whale predation, and per capita effects. Model L is similar to Model K but with group mortality instead of catastrophic mortality. Models M, N, and O are the same as Models A, C, and F, respectively, but with a value for K that declined from 1979 to 1999 and remained low after 1999.
Model scenarios were compared using the Bayes factor (Kass and Raftery, 1995; Wade, 2002), where [H.sub.xy] represents the Bayes factor comparing model scenario x to model scenario y. If model scenarios are considered to have equal prior probability (prior to examination of the data), then the Bayes factor gives the posterior odds ratio for one model scenario compared to another (Kass and Raftery, 1995). In other words, a Bayes factor of [H.sub.xy] = 3 means model scenario x is 3 times more probable, in light of the data (i.e., the abundance estimates), than model scenario y.
Kass and Raftery (1995) provide a useful scale for interpreting the Bayes factor by converting it to twice the natural logarithm of the Bayes factor. With the larger of the two posterior probabilities in the numerator (model scenario x), a value less than 2 indicates no substantial difference between the model scenarios, and values between 2 and 6 are interpreted as positive evidence for one model scenario over the other. Values between 6 and 10 are strong evidence, and values greater than 10 are very strong evidence for model scenario x over model scenario y (Kass and Raftery, 1995). We employ the harmonic mean identity to estimate the relative Bayes factor as the harmonic mean of the distribution of the log likelihoods from the MCMC sample using only data from the posterior sample. The harmonic mean is a poor estimator on its own as it is highly sensitive to the few very low likelihood samples, which the MCMC algorithm will tend to include (Kass and Raftery, 1995). To resolve this we calculated the harmonic means of each posterior data set truncated at various low likelihood values and selected the quantile at which these low likelihood samples begin to dominate the inference for some of the scenarios. Inference was then done using the truncated posterior data set. The value for comparison is then twice the difference between the log Bayes factors of the two models or, in other words, twice the difference of the log of the harmonic means. From the three separate chains, we estimated a standard error for each 2 log Bayes factor estimate. Hypothesis testing was done by comparing the Bayes factor, and examining the posterior distributions.
Results and Discussion
To illustrate the consequences of the different mortality effects in relation to population size, examples of expected annual survival rates by population size are given in Figure 4. Where survival rate is determined only by density dependence, survival increased as the population declined (Fig. 4, solid line), so that the highest survival rate (97% in this example) occurs as the population approaches zero.
When the catastrophic mortality effect is added, because it is independent of population size, in this example there is a 1% reduction in survival throughout the range of population sizes. When killer whale predation or group mortality effects are added to the population model, survival still increases as the population declines at moderate population sizes, but below some threshold population size, survival declines as population numbers decline (Fig. 4). This represents depensation, or an Allee effect, which is known to increase extinction risk for small populations (Courchamp et al., 1999).
Consequently, if demographic stochasticity environmental variability or mortality events pushed the population below the threshold, the population continues to decline to extinction. For example, the group mortality effect intensifies as the population declines because there are fewer groups, so each individual is more likely to be in a group that is affected. When the group mortality effect is included with the constant mortality effect of predation, the threshold effects are combined (Fig. 4, gray dashed line).
Population Model Results
The MCMC numerical integration method makes it possible to calculate a posterior distribution for any output quantity of interest. For each of the model scenarios, we calculated posterior distributions for population size in 50 and 100 years, as well as for the population growth rate over the next 20 years (2014-34). All model scenarios showed a range of outcomes (Fig. 5; Table 6), which could include a decline to extinction, an increase to K, and intermediate values between extinction and K. For the purposes of these results, we define recovery as a population size exceeding 780 whales or 60% of 1,300, the largest historical CIBW abundance estimate that is based on survey data (Calkins (1)). However, the probability of each outcome depended on the parameters of the individual model scenarios.
In the healthy population scenario (Model A, Fig. 5A) there was a 99% probability the population would increase during the next 100 years, with a median value above 1,000 by 2080 and 0% probability of extinction in 100 years. In the per capita mortality scenario (Model B, Fig. 5B), there was an over 50% probability that the population would decline (Table 6), with limited probability of recovery to 780 belugas. This contrast in outcomes was largely predicted by the distribution of population trends over the years 2014-2034, where the Model A growth rate over the next 20 years is nearly all increasing (i.e. >0.0) while Model B on average is near zero and has much of its distribution below 0.0 resulting in population decline (Fig. 5A,B; Table 6).
Models C-E include constant mortality from predation as an additional source of mortality in the healthy model scenario (Table 5; Fig. 5C-E). In all three models this resulted in a small quantitative change such that the population was estimated to increase more slowly and take longer to recover, but, as in Model A, there was a high probability the population would increase and not go extinct, suggesting that the threshold effect was not acting on these population trajectories in most cases. The 20-year average trend (Table 6) for Models C-E, does decrease from Model A, reflecting decreased survival resulting from predation.
Model F combines the per capita mortality with the predation mortality at the observed level showing a slightly broader range of outcomes, but a 6% probability of extinction in 100 years (Fig. 5F; Table 6). Models G and H both include additional mortality as stochastic events resulting in increased uncertainty in the outcomes at 100 years, with the main difference that the threshold effect of the group mortality resulted in a greater risk of extinction (Fig. 5G, H; Table 6). The four models (I-L) are mixtures of the different mortality effects (Table 5). Including predation mortality with the catastrophic mortality, Model I results in similar uncertainty in the outcomes but shifts the range of outcomes toward increased risk of decline and extinction (Fig. 5I; Table 6), while including predation mortality with the group mortality, Model J has the highest risk of extinction of the scenarios (Fig. 5J; Table 6). Models K and L are mixtures of Models I and J with Model B (Table 5) and showed behavior more similar to Models K and L, respectively (Fig. 5K, L; Table 6). The three models with reduced carrying capacity all behaved similarly except that the threshold effect of predation and the per capita decline acting together had an extinction risk 3% in 100 years and predation alone had an extinction risk of 1% indicating that the lower carrying capacity reduced the risk of extinction compared to Models C and F.
Model Scenario Comparison with the Bayes Factor
The inference based on the natural logarithm of harmonic means of the posterior samples of likelihoods truncated at quantiles above 0.01 remained consistent while the inference became quite variable at quantiles below 0.01 thus, the 0.01 quantile was selected as the truncation point (Fig. 6). Samples with likelihoods less than that value of the 0.01 quantile, the lowest 300, in each data set were discarded and inference was conducted using the remaining 29,700 samples. The standard errors for the estimated 2ln(Bayes factors); were all less than 1.0 and the PSRF values were all less than 1.5 indicating that the three chains were consistent with each other and covered the posterior distribution. The alternative model scenarios with the per capita effects or the reduced carrying capacity, Models B, F, M, N, and O had substantially higher posterior probability than Model A, the healthy population scenario (Table 7). All of these alternative model scenarios had 2ln(Bayes factor) values greater than 20. Therefore, there was substantial evidence that the data supported these alternative model scenarios more than the healthy population scenario. Models D, E, H, K, and L also had higher posterior probability than Model A. The scenarios with the highest posterior probabilities were the declining carrying capacity Models M, N, and O.
Using Models B, F, M, N, and O as the basis of comparison, we see that there is positive evidence to prefer these models over all others and, in most cases, there is strong or very strong evidence supporting these models. While the probabilities are similar, the risk of extinction is different, ranging from 0% in 100 years to 6% in 100 years, with differences resulting from the threshold effect of the predation mortality and the declining growth rates in the per capita model. Also, these models predict low probabilities of recovery to 780 belugas. Thus, a better understanding of the change in predation rates with CIBW population size, factors affecting carrying capacity and stressors impacting survival and fecundity are essential to improving the estimates of extinction risk and identifying a path to recovery.
Hypothesis comparison and ranking using the Bayes factors provide alternative views to the same analysis. Under Hypothesis H1, the population has not yet begun to recover but it will under the status quo--there has been a delay (time lag) in recovery due to a depletion of mature females and the population will begin to increase once it rebuilds the mature female segment of the population. The data do not support this hypothesis. The Bayes factor comparison estimates Model B is 8.4x[10.sup.5] times more likely than Model A (Table 7). In other words, the proposed mechanism in Model B is supported, while a deficit of reproductive age females from the biased hunt cannot explain the lack of recovery in the population. Thus, we must conclude that the current lack of recovery is not due to time lags that will resolve in the future.
Under Hypothesis H2, the subsistence hunt was not the sole cause of the decline observed since 1994 and the population has not begun to recover because an as yet unidentified population-wide stressor (e.g., nutritional stress, disease, contaminants, noise, loss of key habitat) has caused a decrease in fecundity and/or survival, resulting in a negative growth rate. The population will not begin to recover in the future without a change in the fecundity and or survival of the population. The data support this hypothesis because Model B was 8.4x[10.sup.5] times more likely than Model A (Table 7), and there was over a 50% probability that the population will continue to decline (Table 6).
Under Hypothesis H3, the population has not begun to recover because the reduction in population size from hunting has led to it falling into a "predator pit" (a numerically constant level of predation was sustained by the population when it was >1,000 animals, but at the current population level, predation prevents recovery). The data do not support this hypothesis. Models D and E were more likely than Model A (Table 7), suggesting that the mortality rates from predation are likely underestimated by carcass counts. However even when the predation mortality risk was doubled (Model D), the anticipated average growth rate over the next 20 years was reduced only by 0.3% from the average for Model A, and the probability of decline remained only 2% (Table 6).
Under Hypothesis H4, the population has not begun to recover because the level of predation from killer whales coincidentally increased around the time of the population decline (e.g., due to a prey shift by the killer whale population). The data do not provide support for this hypothesis for the same reasons stated under H3, while Model E is more likely than Model A, it does not reduce the growth rate sufficiently to prevent recovery.
Under Hypothesis H5, the population has not begun to recover because the reduction in population size from hunting and subsequent retraction in range has led to a greater proportion of the population being vulnerable to "catastrophic" mortality events such as disease outbreaks, oil or toxics spills, volcanic eruptions, or fish run failures. The results do not support this hypothesis in that Model G was 0.0027 less likely than Model A (Table 7), largely the result of the greater variability of growth from year to year.
Under Hypothesis H6, the population has not begun to recover because the reduction in population size from hunting and subsequent retraction in range has led to fewer social groups in the population, such that a greater proportion of the population is vulnerable to mortality events that affect all or a large fraction of a single social group. The results do support this hypothesis in that Model H was 110 times more likely than Model A (Table 7).
Under Hypothesis H7, the population has not begun to recover because of a combination of killer whale predation, per capita effects, and increased catastrophic or group mortality following population decline. The results support this hypothesis in that some of the model scenarios with mixed mortality effects, Models F, K, and L, were more likely than Model A. Models I and J were less likely than Model A suggesting that the mixture with the per capita effect was key, and all of the mixed models were much less likely than Model B (Table 7).
Under Hypothesis H8, the population has not begun to recover because the carrying capacity has declined to a low level. The results support this hypotheses in that Models M, N, and O were substantially more likely than Model A (Table 7). Through hypothesis comparison, we demonstrated that several of the mortality models were more likely than the assumption of a healthy population, but this did not distinguish among the models except that the constant or predation mortality model was not sufficient by itself. The Bayes factors provided a means to rank the model scenarios and select the most likely, the per capita mortality or reduced carrying capacity with or without the threshold effect of constant or predation mortality.
There was some variability in the posterior distributions for the 1994 abundance, with modes ranging from 590 for Model A to 700 for Model O. Model scenarios with higher average growth rates had lower modes while models with negative average growth rates had higher modes (Fig. 7; Table 6). The three highest modes were Models M, N, and O which had declining carrying capacity. All of the results fell well within the sampling distribution for the 1994 abundance estimate, indicating that it was not very selective in the likelihood. Posterior distributions for the 2014 abundance distributions closely matched the sampling distribution of the 2014 abundance estimate in shape, but were somewhat lower, while the healthy population scenario and killer whale predation model scenarios favored higher values, with modes of 50 belugas greater than the 2014 abundance estimate (Fig. 8). The stochastic mortality effects broadened the distributions for Models G-J but when the per capita effect was included, Models K and L, the posterior distributions were more similar to Models B and F.
As discussed above, the current lack of recovery in the population cannot be explained by a deficit of adult females caused by the subsistence hunt, as the healthy population scenario was shown to have a much poorer fit to the data. This occurs because the healthy scenario is forced (via the prior distribution) to have a positive intrinsic rate of increase, and this results in an increasing population trajectory (Fig. 9), but the trend of the abundance estimates indicates the population is not growing. Note that the posterior distribution of the intrinsic growth rate for Model A is skewed toward its lower bound indicating that the lower limit of the prior distribution was informative. The constraints on the upper limit of survival rate and birth rate limited the intrinsic growth rate to be less than 1.031, however the posterior distribution of Model A had declined to low probability so that these constraints were not informative to Model A.
A comparison of the distributions of the number of females and the fraction of females in the population trajectories for each model, shows a moderately depleted adult female fraction in the Model A trajectories during the years 1997-2003, but this is fully recovered by 2014 (Fig. 10). The results for Model B indicate that the male fraction was more strongly depleted than the female fraction but this is also resolved by 2014. In other words, while the deficit of reproductive age females is a plausible mechanism for slowing recovery, the trajectory for the numbers of adult females in Model A (Fig. 10) shows an upturn within the first few years after the end of uncontrolled hunting in 1999. This result for adult females is corroborated by the abundance trajectories for these years where the Model A trajectories begin to increase after 1999 (Fig. 11), indicating that any delay in recovery that would have occurred from a deficit of adult females would have lasted for only a few years.
One interesting outcome is that the best fit for Models A, C, D, and E occurred by estimating that a larger fraction of adult females were taken in the hunt (Fig. 12). This is a case of two sources of data conflicting. For example, Model A imposes an increasing population growth rate, which conflicts with the observed abundance trend. The growth rate can be partially offset for a few years (Fig. 11) by having a higher proportion of females taken in the hunt than actually observed, as this will slow the model population increase.
Model scenarios with a catastrophe or a group mortality effect had greater variability, such that there was a broader range of population trajectories possible. This is evident in the probabilities of recovery and extinction, and the range of growth rates from 2014 to 2034 (Fig. 5G, H; Table 6). While this variability might allow the trajectories to match the point estimates of abundance closely, it also could push the trajectory away from the abundance time series with one or two large events.
In both Models G (Fig. 13) and H (Fig. 14), the posterior distributions of the mortality rate for an event showed a preference for lower expected mortality rates (Figs. 13c, 14c). For catastrophic mortality, the histogram indicates a preference for the lower range of 10-20% mortality per event (Fig. 13a). For the group mortality, the mode is near 25% mortality per event but the distribution is quite broad (Fig. 14a); however, this is applied to a group that on average is 20% of the population and averages approximately 5% mortality at the population level. Even at 100% mortality, this only represents a 20% mortality to the population. Thus, the two survival models support similar ranges of mortality at the population level.
The low mortality rate for events limits the variability of the trajectories to the lowest values possible for these model scenarios. The probability of events was then determined by the change in survival rate required to achieve the change in growth rate from the healthy population scenario. For the catastrophic mortality scenario, this resulted in a fairly strong preference for lower probability of events, indicating that the selection in Model G was stronger on the event probability rather than the mortality rate (Fig. 13 a, b). In contrast, for the group mortality scenario in Model H, the probability of an event was less strongly selected and resulted in a posterior distribution mode near 20% but ranging from 0% to 100% (Fig. 14b). When these are combined into the expected annual mortality, the posterior distributions have lower medians, with the median and distribution for the expected annual group mortality about twice as broadly distributed (Figs. 13c, 14c).
In general, when the predation mortality effect with [C.sub.c] = 1 was added to the healthy population scenario (Model A vs. Model C), there was little change in the population outcomes, but it improved the fit to the data. However, when [C.sub.c] = 1 was included in the model scenario with per capita effects on fecundity or survival (Model B vs. Model F), there was no improvement in the fit to the data, and, in fact, the fit was reduced, but it added 6% to the probabilities of extinction in 100 years. The group and catastrophic mortality scenarios, with or without predation mortality added, did not provide a good fit to the data by themselves, and received little posterior probability. However, when per capita effects were added, these scenarios improved in fit over Model A, but remained poorer fits than Model B, suggesting that of the four types of effects the per capita effects were the most likely (Table 7).
Little information is available for the model on the carrying capacity of Cook Inlet for belugas and consequently a range of posterior distributions resulted for the different model scenarios. In general, the models without either a decline in carrying capacity or a per capita effect (Models A, C-E, G-J) selected the lower values for K (Fig. 15). Models with a change in carrying capacity (Models M-O) and a per capita effect (Models B, F, K, L, O) selected higher values for K, with Model O, which included both scenarios, showing the strongest selection for high values (Fig. 15). For the three models that allowed a change in K over time, the two models that did not include a per capita effect had posterior distributions for [K.sub.1999] that closely matched the sampling distribution for the abundance estimate in 2014 (Fig. 16). Model O (with the per capita effect and decreasing K) did select the range of the 2014 abundance estimate, but also chose the range between 400 and 800 with nearly equal preference. This difference between the models with and without the per capita effect is also evident in the projections for the next 100 years, with the majority of cases for Models M and N showing the CIBW population remaining stable at a low level while Model O shows a significant number of declining population trajectories. None of the models selected values below 250 indicating that if the carrying capacity has declined it is not less than the current population size.
Extinction Risk and Potential for Recovery
The key purpose for this analysis was to estimate the risk of extinction and to compare potential mortality models for the endangered CIBW population. The five most likely scenarios, Models B, F, M, N, and O, are all much more likely than any of the other model scenarios, therefore, they are the preferred scenarios to estimate extinction risk. These models fit the data well because they all indicate there has been a substantial decline in the population growth rate from an expected healthy population model. In the case of Models B and F, it is through a per capita decline in fecundity or survival, whereas, in Models M, N, and O, it is through a reduction in carrying capacity. The reduction in carrying capacity in turn reduces the per capita survival and fecundity through the density dependence of the population model, as the model essentially estimates the population is close to carrying capacity.
To the five most likely models we add Model K to account for the risk of catastrophic events. While these events are known to occur, there are few data on their impact on beluga populations that could be included in the current analysis. Consequently, the likelihood of this model was not as well supported by the data used. It has, however, the highest likelihood of the models that include unusual mortality events. From these models, we estimate that the risk of extinction in 100 years is between 0% and 14% with the main uncertainty resulting from the strength of the threshold resulting from the constant mortality effect.
In 2008, Hobbs and Shelden (5) estimated the range of extinction risk to be between 1% and 27% in 100 years from a similar set of models. In our current analysis, Models B, F, and K correspond to Models a, d, and h in Hobbs and Shelden (5), respectively. There are no equivalents for Models M, N, and O in Hobbs and Shelden5 and the current analysis includes no equivalents for their Models c, e, and g. Models a, d, and h in Hobbs and Shelden (5) had extinction risks of 1%, 12%, and 26% in 100 years, respectively; consequently, they spanned the range of extinction risks (i.e., 1-27%) in that paper, but are roughly twice the values for the corresponding models in our current analysis.
There are several reasons for these changes. First, we now have 20 abundance estimates, up from the 15 estimates used in the earlier analysis. These additional data supported the selection of annual growth rates close to zero; consequently, there is a more precise estimation of the growth rate parameter with fewer outliers that would be low and trend towards extinction. Second, the mean growth rate for Model B is 0.1%/yr with 53% of cases declining, while the corresponding Model a in Hobbs and Shelden (5) had a mean growth rate of -0.4%/yr and 62% of cases declining, which supported this slightly improved outlook. In other words, with more data the recent trend of the population has been estimated more precisely and it has been found to be roughly stable or slowly declining. Thus, although the population is not increasing as expected, the data indicate the population is not declining precipitously, and, therefore, the probability of extinction is lower than in the previous analyses (i.e., Hobbs and Shelden (5)).
Models K and h were not identical in that in Model h of Hobbs and Shelden (5) the mortality rate was fixed at 20% and the probability of an event was fixed at 5%; consequently, the expected mortality was 1%/yr. In Model K these percentages were allowed to vary with the mortality rate chosen from between 10% and 50%, and the probability selected so that the expected mortality per year met a value selected by the model for that case. Model K also included per capita mortality and predation mortality, so the three were competing. In Model K, the average values of the probability of an event was 3.1% and the expected annual mortality from catastrophic events was 0.7%; therefore, Model K selected fewer events with a lower expected mortality rate than in Model h, which resulted in a decreased risk of extinction. Thus the effect of catastrophic mortality events was reduced in Model K compared to model h in Hobbs and Shelden (5).
The strong selection of the five model scenarios (Models B, F, M-O) indicates that the observed decline of the population is likely the result of a per capita or chronic decrease in survival or reproduction or that carrying capacity has declined. While we have provided examples of long-term habitat changes, reduction of available forage, and anthropogenic effects such as introduced toxins, there is also the possibility that the reduction in the fecundity or survival of the population is the result of reduction in population size itself.
There are many examples of population level consequences beyond the numerical effect of the removal of individuals caused by social disruption from hunting. For example, hunting elephants, Loxodonta africana, led to lower fecundity in several elephant populations (Poole, 1987; Gobush et al., 2008). Part of the explanation for this was the loss of "social knowledge" from the removal of matriarchs for their large tusks (McComb et al., 2001). Similarly, lower fecundity has been seen in a population of bighorn sheep, Ovis canadensis, where rams with larger horns were removed (Coltman et al., 2003). Heavy hunting pressure on wolves decreased the average size of social groups and led to lower natural survival for a variety of reasons (Haber, 1996). Specifically, a review of multiple studies showed that removal of breeding wolves, Canis lupus, led to decreased wolf packs or to the dissolution of packs, and that pup survival was higher in larger packs and was correlated with the presence of auxiliary nonbreeders (Brainerd et al., 2008). Heavy hunting pressure on cougar populations, Puma concolor, was correlated with increased immigration, reduced kitten survival, and reduced female population growth (Cooley et al., 2009).
At least two examples exist showing lower fecundity in cetaceans after hunting and social disruption. Fecundity was very low in Galapagos sperm whales, Physeter macrocephalus, after the population was severely depleted during commercial whaling, apparently due to a lack of mature males (whitehead et al., 1997). Spinner dolphins, Stenella longirostris, subjected to chase and encirclement by yellowfin tuna, Thunnus albacares, purse seiners in the eastern tropical Pacific have experienced a decline in fecundity (Cramer et al., 2008), and this could be due to the loss of breeding males, which represent only a tiny fraction of all mature males (Perrin and Mesnick, 2003). There was also apparently lower calf survival caused by the separation of mothers and dependent calves during the chase (Archer et al., 2001). The CIBW population underwent heavy hunting pressure for at least 20 years (1979-99), which caused the population to decline from approximately 1,300 whales to about 300 whales, and the number of adult females in the population may have declined to less than 100 (Fig. 10). Social disruption from the hunt causing a decline in fecundity or survival in the population is plausible, however, we did not attempt to develop a model specific to this mechanism.
Some examples from other mammalian populations suggest that once a hunt is limited or ended, the population would rebound. Other examples, such as the loss of matriarchs in elephants, dominant breeding males in sperm whales, or optimal pack structure in wolves, suggest there is a lag of 1-2 generations, during which the population continues to decline, before the population recovers. These small population effects and the threshold effects modeled above, may act at population sizes near to the current population size, suggesting that a significantly larger population is necessary to keep the population from becoming extinct in the next 100 years.
The ultimate goal for the CIBW population is not just avoiding extinction but recovering to a sustainable size and distribution. To this end, we can draw two conclusions from the results of this modeling effort. First, if the intrinsic growth rate can be increased above 2%, then recovery is very likely. Second, one or more per capita effects are reducing survival and/or fecundity or substantial carrying capacity has been lost. While there was considerable variation in extinction risk at 100 years among the models, the assumption of an intrinsic rate of growth greater than 2% removed the risk of decline and extinction even when the threshold effect of predation mortality was included. Thus, research intended to support recovery of this population should be directed toward identifying and reversing these per capita effects. Also, a better understanding of the factors determining carrying capacity is required. The possible presence of a threshold in the growth rate resulting in an Allee effect, increases the urgency of these efforts, as reversal of the trend may become far more difficult if the population declines further. Taken as a whole, these modeling results indicate clearly that it is likely that the CIBW population will continue to decline or go extinct unless factors determining its growth and survival are altered in its favor.
Funding for this project was provided to the Alaska Fisheries Science Center's National Marine Mammal Laboratory by the NMFS Marine Mammal Assessment Program as well as through funds and logistical support from the NMFS Alaska Regional Office. Document reviews were provided by Brian Fadely (NMML), G. Pendleton (ADFG), G. Duker, and J. Lee (AFSC Publications Unit), and two anonymous reviewers.
Allee, W. C., O. Park, A. E. Emerson, T. Park, and K. P. Schmidt. 1949. Principles of animal ecology. W. B. Saunders Co, Philadelphia, PA, 837 p.
Allen, B. M., and R. P. Angliss. 2013. Alaska marine mammal stock assessments, 2012. U.S. Dep. Commer., NOAA Tech. Memo. NMFS AFSC-245, 282 p.
Archer, F, T. Gerrodette, A. Dizon, K. Abella, and S. Southern. 2001. Unobserved kill of nursing dolphin calves in a tuna purse seine fishery. Mar. Mammal Sci. 17:540-554. (doi: 10.1111/j.1748-7692.2001.tb01003.x).
Begon, M., J. L. Harper, and C. R. Townsend. 1996. Ecology: individuals, populations and communities, 3rd ed. Blackwell Sci., Oxford, U.K., 247 p. (doi: 10.1002/9781444313765).
Brainerd, S. M., H. Andren, E. E. Bangs, E. H. Bradley, J. A. Fontaine, W. Hall, Y. Iliopoulos, M. D. Jimenez, E. A. Jozwiak, O. Liberg, C. M. Mack, T. J. Meier, C. C. Niemeyer, H. C. Pedersen, H. Sand, R. N. Schultz, D. W. Smith, P. Wabakken, and A. P. Wydeven. 2008. The effects of breeder loss on wolves. J. Wildl. Manage. 72:89-98 (doi: 10.2193/2006-305).
Brandon, J., and P. R. Wade. 2006. Assessment of the Bering-Chukchi-Beaufort Seas stock of bowhead whales using Bayesian model averaging. J. Cetacean Res. Manage. 8(3):225-239
Breiwick, J. M., L. L. Eberhardt, and H. W. Braham. 1984. Population dynamics of western Arctic bowhead whales (Balaena mysticetus). Can. J. Fish. Aquat. Sci. 41(3):484-496. (doi: 10.1139/f84-058).
Brooks, S. P., and A. Gelman. 1998. General methods for monitoring convergence of iterative simulations. J. Comput. Graph. Stat. 7:434-455.
Citta, J. J., L. T. Quakenbush, K. J. Frost, L. F. Lowry, R. C. Hobbs, and H. Aderman. In press. Movements of beluga whales (Delphinapterus leucas) in Bristol Bay, Alaska. Mar. Mammal Sci.
Coltman, D. W., P. O'Donoghue, J. T. Jorgenson, J. T. Hogg, C. Strobeck, and M. Festa-Bianchet. 2003. Undesirable evolutionary consequences of trophy hunting. Nature 426:655-658. (doi: 10.1038/nature02177). Cooley, H. S., R. B. Wielgus, G. M. Koehler, H. S. Robinson, and B. T. Maletzke. 2009. Does hunting regulate cougar populations? A test of the compensatory mortality hypothesis. Ecol. 90:2913-2921 (doi: 10.1890/08-1805.1).
Courchamp, F., T. Clutton-Brock, and B. Grenfell. 1999. Inverse density dependence and the Allee effect. Trends Ecol. Evol. 14(10):405-410. (doi: 10.1016/S0169-5347(99)01683-3).
Cramer, K. L., W. L. Perryman, T. Gerrodette. 2008. Declines in reproductive output in two dolphin populations depleted by the yellowfin tuna purse-seine fishery. Mar. Ecol. Prog. Ser. 369:273-285. (doi: 10.3354/meps07606).
de Laguna, F. 1975. The archaeology of Cook Inlet, Alaska, 2nd ed. Alaska Hist. Soc., Anchorage, [Orig. publ. 1934. Univ. Pennsylvania], 264 p.
Eggers, D. M., and J. R. Irvine. 2007. Trends in abundance and biological characteristics for north Pacific sockeye salmon. N. Pac. Anadr. Fish Comm. Bull. 4: 53-75.
Ellison, A. M. 1996. An introduction to Bayesian inference for ecological research and environmental decision-making. Ecol. Appl. 6: 1036-1046. (doi: 10.2307/2269588).
Euler, L. 1760. A general investigation into the mortality and multiplication of the human species. Transl. N. Keyfitz and B. Keyfitz for Theoretical Pop. Biol. (1970)1:307-314 (doi:10.1016/0040-5809(70)90048-1).
Ford, J. K. B., G. M. Ellis, P. F. Olesiuk and K. C. Balcomb. 2010. Linking killer whale survival and prey abundance: food limitation in the oceans' apex predator? Biol. Lett. 6(1):139-142. (doi: 10.1098/rsbl.2009.0468).
Frost, K. J., L. F. Lowry, and R. R. Nelson. 1985. Radiotagging studies of belukah whales (Delphinapterus leucas) in Bristol Bay, Alaska. Mar. Mammal Sci. 1:191-202. (doi: 10.1111/ j.1748-7692.1985.tb00008.x).
Gelman, A., J. B. Carlin, H. S. Stern, and D. B. Rubin. 1995. Bayesian data analysis. Chapman and Hall, Lond., 526 p.
Gobush, K. S., B. M. Mutayoba, and S. K. Wasser. 2008. Long-term impacts of poaching on relatedness, stress physiology, and reproductive output of adult female African elephants. Conserv. Biol. 22(6):1590-1599. (doi: 10.1111/j.1523-1739.2008.01035.x).
Goetz, K. T., D. J. Rugh, A. J. Read, and R. C. Hobbs. 2007. Habitat use in a marine ecosystem: Beluga whales Delphinapterus leucas in Cook Inlet, Alaska. Mar. Ecol. Prog. Ser. 330:247-256. (doi: 10.3354/meps330247).
--, R. A. Montgomery, J. M. Ver Hoef, R. C. Hobbs, and D. S. Johnson. 2012. Identifying essential summer habitat of the endangered beluga whale Delphinapterus leucas in Cook Inlet, Alaska. Endang. Species Res. 16:135-147. (doi: 10.3354/esr00394).
Goodman, D. 1982. Optimal life histories, optimal notation, and the value of reproductive value. Am. Nat. 119:803-823.
Haber, G. C. 1996. Biological, conservation, and ethical implications of exploiting and controlling wolves. Conserv. Biol. 10:1068-1081. (doi: 10.1046/j.1523-1739.1996.10041068.x).
Heide-Jorgensen, M. P., P. Richard, M. Ramsay, and S. Akeeagok. 2002. Three recent ice entrapments of Arctic cetaceans in West Greenland and the eastern Canadian High Arctic. NAMMCO Sci. Publ. 4:143-148. (doi: 10.7557/3.2841).
Hobbs, R. C., D. J. Rugh, and D. P. DeMaster. 2000a. Abundance of beluga whales, Delphinapterus leucas, in Cook Inlet, Alaska, 1994-2000. Mar. Fish. Rev. 62(3):37-45.
--, J. M. Waite, and D. J. Rugh. 2000b. Beluga, Delphinapterus leucas, group sizes in Cook Inlet, Alaska, based on observer counts and aerial video. Mar. Fish. Rev. 62(3):46-59.
--, K. L. Laidre, D. J. Vos, B. A. Mahoney, and M. Eagleton. 2005. Movements and area use of belugas, Delphinapterus leucas, in a subarctic Alaskan estuary. Arctic 58(4):331-340. (doi: 10.14430/arctic447).
--, K. E. W. Shelden, D. J. Rugh, C. L. Sims, and J. M. Waite. 2015a. Estimated abundance and trend in aerial counts of beluga whales, Delphinapterus leucas, in Cook Inlet, Alaska, 1994-2012. Mar. Fish. Rev. 77(1):11-31. (doi: 10.7755/MFR.77.1.2).
--, C. L. Sims, K. E. W. Shelden, L. Vate Brattstrom, and D. J. Rugh. 2015b. Annual calf indices for beluga whales, Delphinapterus leucas, in Cook Inlet, Alaska, 2006-12. Mar. Fish. Rev. 77(2):40-58. (doi: 10.7755/MFR.77.2.3).
Huntington, H. P. 2000. Traditional knowledge of the ecology of belugas, Delphinapterus leucas, in Cook Inlet, Alaska. Mar. Fish. Rev. 62(3):134-140.
Johnson, N. L., and S. Kotz. 1970. Distributions in statistics: continuous univariate distributions. Vol. 2. John Wiley and Sons, N.Y., 306 p.
Kass, R. E., and A. E. Raftery. 1995. Bayes Factors. J. Am. Stat. Assoc. 90(430):773-795.
King, R., B. Morgan, O, Gimenez, and S. Brooks. 2010. Bayesian analysis for population ecology. CRC Press, Boca Raton, Fla., 442 p.
Litzky, L. K. 2001. Monitoring recovery status and age structure of Cook Inlet, Alaska beluga whale by skin color determination. M.S. thesis, Univ. Wash., Seattle, 76 p.
Lockyer, C., A. A. Hohn, D. W.Doidge, M. P. Heide-Jorgensen, and R. Suydam. 2007. Age determination in belugas (Delphinapterus leucas): a quest for validation of dentinal layering. Aquatic Mammals 33(3):293-304. (doi: 10.1578/AM.33.3.2007.293).
Lotka, A. J. 1907. Relation between birth rates and death rates. Science 26:21-22.
Lowry, L. F, K. J. Frost, A. Zerbini, D. Demaster, and R. R. Reeves. 2008. Trend in aerial counts of beluga or white whales (Delphinapterus leucas) in Bristol Bay, Alaska, 1993-2005. J. Cetacean Res. Manage. 10(3):201-207.
MacCall, A. D. 1990. Dynamic geography of marine fish populations. Univ. Wash. Press, Seattle, 153 p.
Mahoney, B. A., and K. E. W. Shelden. 2000. Harvest history of beluga whale, Delphinapterus leucas, in Cook Inlet, Alaska. Mar. Fish. Rev. 62(3):124-133.
Matkin, C. O., E. L. Saulitis, G. M. Ellis, P, Olesiuk, and S. D. Rice. 2008. Ongoing population-level impacts on killer whales Orcinus orca following the 'Exxon Valdez' oil spill in Prince William Sound, Alaska. Mar. Ecol. Prog. Ser. 356:269-281. (doi: 10.3354/ meps07273)
McComb K., C. Moss, S. M. Durant, L. Baker, S. Sayialel. 2001. Matriarchs as repositories of social knowledge in African elephants. Science 292:491-494
Miller, T. P., R. G. McGimsey, D. H. Richter, J. R. Riehle, C. J. Nye, M. E. Yount, and J. A. Dumoulin. 1998. Catalog of the historically active volcanoes of Alaska. U.S. Geol. Surv. Open-File Rep. OF 98-0582, 104 p.
Moore, S. E., K. E. W. Shelden, L. K. Litzky, B. A. Mahoney, and D. J. Rugh. 2000. Beluga, Delphinapterus leucas, habitat associations in Cook Inlet, Alaska. Mar. Fish. Rev. 62(3):60-80.
Morris, W. F., and D. F. Doak. 2002. Quantitative conservation biology: Theory and practice of population viability analysis. Sinauer Associates, Inc., Sunderland, MA.
NOAA. 2000. Designating the Cook Inlet, Alaska, stock of beluga whale as depleted under the Marine Mammal Protection Act (MMPA). 65 Fed. Regist. 34590 (31 May 2000), p. 34,590-34,597 (avail. online at: https://federalregister.gov/a/00-13371).
--. 2008. Endangered and threatened species; endangered status for the Cook Inlet beluga whale. 73 Fed. Regist. 62919 (22 Oct. 2008), p. 62,919-62,930 (avail. online at: https://federalregister.gov/a/E8-25100).
O'Corry-Crowe, G. M., R. S. Suydam, A. Rosenberg, K. J. Frost, and A. E. Dizon. 1997. Phylogeography, population structure and dispersal patterns of the beluga whale, Delphinapterus leucas, in the western Nearctic revealed by mitochondrial DNA. Mol. Ecol. 6:955-970. (doi: 10.1046/j.1365-294X.1997.00267.x).
--, A. E. Dizon, R. S. Suydam, and L. F. Lowry. 2002. Molecular genetic studies of population structure and movement patterns in a migratory species: the beluga whale, Delphinapterus leucas, in the western Nearctic. In C. J. Pfeiffer (Editor), Molecular and cell biology of marine mammals, p. 53-64. Krieger Publ. Co., Malabar, Fla.
Perrin, W. F, and S. L. Mesnick. 2003. Sexual ecology of the spinner dolphin, Stenella longirostris: geographic variation in mating system. Mar. Mammal Sci. 19:462-483. (doi: 10.1111/j.1748-7692.2003.tb01315.x).
Poole, J. H. 1987. Rutting behavior of African elephants: the phenomenon of musth. Behaviour 102:283-315.
Reilly, S. B., and J. Barlow. 1986. Rates of increase in dolphin population size. Fish. Bull. 84:527-533.
Rugh, D. J., K. E. W. Shelden, and B. A. Mahoney. 2000. Distribution of belugas, Del phinapterus leucas, in Cook Inlet, Alaska, during June/July, 1993-2000. Mar. Fish. Rev. 62(3):6-21.
--, B. A. Mahoney, and B. K. Smith. 2004. Aerial surveys of beluga whales in Cook Inlet, Alaska, between June 2001 and June 2002. U.S. Dep. Commer., NOAA Tech. Memo. NMFS-AFSC-145, 26 p.
--, K. E. W, Shelden, C. L. Sims, B. A. Mahoney, B. K. Smith, L. K. Litzky, and R. C. Hobbs. 2005. Aerial surveys of belugas in Cook Inlet, Alaska, June 2001, 2002, 2003, and 2004. U.S. Dep. Commer., NOAA Tech Memo. NMFS-AFSC-149, 71 p.
--, --, and R. C. Hobbs. 2010. Range contraction in a beluga whale population. Endang. Species Res. 12:69-75 (doi: 10.3354/esr00293).
Schwacke, L. H., C. R. Smith, F. I. Townsend, R. S. Wells, L. B. Hart, B. C. Balmer, T. K. Collier, S. De Guise, M. M. Fry, L. J. Guillette, S. V Lamb, S. M. Lane, W. E. McFee, N. J. Place, M. C. Tumlin, G. M. Ylitalo, E. S. Zolman, and T. K. Rowles. 2014. Health of common bottlenose dolphins (Tursiops truncatus) in Barataria Bay, Louisiana, following the Deepwater Horizon oil spill. Environ. Sci. Technol. 48 (1):93-103. (doi: 10.1021/ es403610f).
Shelden, K. E. W., D. J. Rugh, B. A. Mahoney, and M. E. Dahlheim. 2003. Killer whale predation on belugas in Cook Inlet, Alaska: Implications for a depleted population. Mar. Mammal Sci. 19(3):529-544.
--, --, K. T. Goetz, C. L. Sims, L. Vate Brattstrom, J. A. Mocklin, B. A. Mahoney, B. K. Smith, and R. C. Hobbs. 2013. Aerial surveys of beluga whales, Delphinapterus leucas, in Cook Inlet, Alaska, June 2005 to 2012. U.S. Dep. Commer., NOAA Tech. Memo. NMFS-AFSC-263, 122 p.
--, K. T. Goetz, D. J. Rugh, D. G. Calkins, B. A. Mahoney, and R. C. Hobbs. 2015. Spatio-temporal changes in beluga whale, Delphinapterus leucas, distribution: results from aerial surveys (1977-2014), opportunistic sightings (1975-2014), and satellite tagging (1999-2003) in Cook Inlet, Alaska. Mar. Fish. Rev. 77(2): 1-31. (doi: 10.7755/MFR.77.2.1).
Siegstad, H., and M. P Heide-Jorgensen. 1994. Ice entrapments of narhwals (Monodon monoceros) and white whales (Delphinapterus leucas) in Greenland. Meddr Gronland, Bioscience 39:151-160.
Stanek, R. T. 1994. The subsistence use of beluga whale in Cook Inlet by Alaska Natives, 1993. Alaska Dep. Fish Game, Div. Subsistence, Tech. Pap. 232, 24 p.
Stewart, R. E. A., S. E. Campana, C. M. Jones, and B. E. Stewart. 2006. Bomb radiocarbon dating calibrates beluga (Delphinapterus leucas) age estimates. Can. J. Zool. 84(12):1840-1852. (doi: 10.1139/z06-182).
Taylor, B. L., and D. P DeMaster, 1993. Implications of non-linear density dependence. Mar. Mammal Sci. 9:360-371. (doi: 10.1111/ j.1748-7692.1993.tb00469.x).
Vos, D. J., and K. E. W. Shelden. 2005. Unusual mortality in the depleted Cook Inlet beluga (Delphinapterus leucas) population. Northwest. Nat. 86(2):59-65. (doi: 10.1898/1051-1733(2005)086[0059:UMITD C]2.0.CO;2).
Wade, P. R. 1998. Calculating limits to the allowable human caused mortality of cetaceans and pinnipeds. Mar. Mammal Sci. 14(1):1-37. (doi: 10.1111/j.1748-7692.1998. tb00688.x).
--. 2000. Bayesian methods in conservation biology. Conserv. Biol. 14: 1308-1316.
--. 2002. A Bayesian stock assessment of the eastern Pacific gray whale using abundance and harvest data from 1967 to 1996. J. Cetacean Res. Manage. 4:85-98.
--. 2009. Population dynamics. In W. F, Perrin, B. Wursig, and J. G. M. Thewissen (Editors), Encyclopedia of marine mammals, 2nd ed., p. 913-918. Academic Press, San Diego.
Ward, E. J. E. E. Holmes, and K. C. Balcomb. 2009. Quantifying the effects of prey abundance on killer whale reproduction. J. Applied Ecol. 46: 632-640. (doi: 10.1111/j.1365-2664.2009.01647.x).
Whitehead, H., J. Christal, and S. Dufault 1997. Past and distant whaling and the rapid decline of sperm whales off the Galapagos Islands. Conserv. Biol. 11:1387-1396. (doi: 10.1046/j.1523-1739.1997.96246.x).
Rod Hobbs (email@example.com), Paul Wade, and Kim Shelden are with the National Marine Mammal Laboratory, Alaska Fisheries Science Center, National Marine Fisheries Service, NOAA, 7600 Sand Point Way NE, Seattle, WA 98115-6349.
(1) Calkins, D. G. 1989. Status of belukha whales in Cook Inlet. In L. E. Jarvela and L. K. Thorsteinson (Editors), Proceedings of the Gulf of Alaska, Cook Inlet, and North Aleutian Basin Information update meeting, 7-8 Feb. 1989, Anchorage, Alaska, p. 109-112. U.S. Dep. Inter., Minerals Manage. Serv., OCS Study, MMS 89-0041.
(2) Dischner, M. 2013. Inlet drift fishermen file suit over salmon management. Alaska J. Commer., Feb. I ss. 3 [http ://www.alaskaj ournal. com/Alaska-Journal-of-Commerce/FebruaryIssue-3-2013/Inlet-drift-fishermen-file-suit-oversalmon-management/].
(3) NMFS. 2015. Draft recovery plan for the Cook Inlet beluga whale (Delphinapterus leucas). U.S. Dep. Commer., NOAa, Natl. Mar. Fish. Serv., Alaska Region. Off., Protected Resour. Div., Juneau, AK, 274 p.
(4) Shelden, K. E. W., C. L. Sims, L. Vate Brattstrom, K. T. Goetz, and R. C. Hobbs. 2015. Aerial surveys of beluga whales (Delphinapterus leucas) in Cook Inlet, Alaska, June 2014. AFSC Proc. Rep. 2015-03, 55 p. Available at http://www.afsc.noaa.gov/Publications/ProcRpt/ PR2015-03.pdf accessed 8 Sept. 2015.
(5) Hobbs, R. C., and K. E. W. Shelden. 2008. Supplemental status review and extinction assessment of Cook Inlet belugas (Delphinapterus leucas). AFSC Processed Rep. 2008-08, 76 p. Available at http://www.afsc.noaa.gov/Publications/ProcRpt/pR2008-08.pdf, accessed 17 Sept. 2015.
(6) An unusual mortality event is defined under the Marine Mammal Protection Act as "a stranding that is unexpected; involves a significant die-off of any marine mammal population; and demands immediate response." See http://www. nmfs.noaa.gov/pr/health/mmume/ accessed 14 May 2014.
(7) Alaska Department of Environmental Conservation. 2011. Major oil spills to coastal waters. Spill Prevention and Response. Division of Water. Available at http://www.dec.state.ak.us/water/wqsar/index.htm accessed 19 Mar. 2014.
(8) Map of the Exxon Valdez Oil Spill, Exxon Valdez Oil Spill Trusties Council website http:// www.evostc.state.ak.us/index.cfm?FA=facts. map, accessed 19 Mar. 2014.
Table 1.--Time series of Cook Inlet beluga whale data used in the Bayesian analysis. Estimated abundance was calculated from observer and video data. Hunt landings and struck and lost data were from Mahoney and Shel-den (2000) and NMFS Alaska Regional Office (NMFS, text footnote 3). Where conflicting sources occurred, all are listed and the range of values is used. Note that killed but lost are included with the struck and lost. Year Estimated Abundance Hunt landings abundance CV (struck and lost) 1994 653 0.24 19(2) 1995 491 0.21 60(14), 52(22), 42(26) 1996 594 0.20 49(49-98) 1997 440 0.13 35(30-40), 35(35) 1998 347 0.17 21(21) 1999 367 0.09 0(0) 2000 435 0.14 0(0) 2001 386 0.10 1(0) 2002 313 0.10 1(0) 2003 357 0.08 1(0) 2004 366 0.13 0(0) 2005 278 0.10 2(0) 2006 302 0.10 0(0) 2007 375 0.08 0(0) 2008 375 0.11 0(0) 2009 321 0.11 0(0) 2010 340 0.08 0(0) 2011 284 0.09 0(0) 2012 312 0.13 0(0) 2014 340 0.08 0(0) Table 2.--Review of female beluga life history parameters from the published literature. Growth layer groups (GLGs) obtained from beluga teeth were reported instead of ages. Past studies (e.g., source 1) used two GLG layers to represent one year of growth; however, recent studies support one GLG per year of growth (Stewart et al., 2006; Lockyer et al., 2007; NAMMCO: http://septentrio.uit.no/index.php/ NAMMCOSP/issue/view/236). Parameters Data Age at sexual 7-13 GLGs (mean = 10, no sample size), 5-6 to maturity 11-12 GLGs (mean = 9, n = 33) 9-11 GLGs (mean = 10, no sample size) 8-9 GLGs (0%), 10-11 GLGs (33%), 12-13 GLGs (94%), 16-17 GLGs (100%)(n = 207) 9.1 [+ or -] 2.8 GLGs (n = 23) 50% at 8.25 GLGs (n = 87) Age at color 12 GLGs (minimum age) change (gray to 14 GLGs (minimum from Mackenzie Delta), white) 17 GLGs (minimums from western Hudson Bay) 9-10 GLGs for males, 10-12 GLGs for females Age at 1st 54% at 8-9 GLGs (n = 12 of 22) conception 41 % at 10-11 GLGs (n = 9 of 22) 4% at 12-13 GLGs (n = 1 of 22) 8.27 GLGs (SE = 2.88, n = 87) Age at senescence 42-43 GLGs (arbitrarily assumed by Kleinenberg) 40 GLGs (corpora level off and decline) Pregnancy and With small fetuses: 0-11 GLGs (0.055), 12-21 GLGs birth rates (0.414), 22-45 GLGs (0.363), 46-57 GLGs (0.267), 58-77 GLGs (0.190). With full-term/neonate: 0-11 GLGs (0), 12-21 GLGs (0.326), 22-45 GLGs (0.333), 46-51 GLGs (0.278), 52-57 GLGs (0.182), 58-77 GLGs (0.125) 0.41 (w/small fetuses); 0.56 (w/full term fetuses or neonates) Lifespan 60-61 GLGs 50-53 GLGs >60 GLGs (oldest female estimated at 70+ GLGs) 46 GLGs (male, tooth worn w/no visible neonatal line) 57 GLGs (female) Adult annual 0.9064 (average based on mean annual mortality survival rate = 0.0936) 0.91-0.92 0.842 and 0.905 (assuming 2GLGs/yr vs. 1 GLG/yr) 0.96-0.97 0.935 Immature annual 0.905 (for neonates in first half year of life, survival mortality rate = 0.095) 0.955 (based on pilot whale net recruitment) Reproductive rate 0.13 (ratio of calves to adult females, modeled) 0.143 (ratio of calves to adult females) 0.114-0.117 (ratio of calves to whales) 0.104 (a model population of 1,000 that included 94 calves) 0.097 (ratio of calves to whales) 0.08-0.10 (ratio of calves to whales) 0.12 (ratio of calves to whales) 0.056-0.10 (ratio of calves to whales) 0.08-0.14 (ratio of calves to whales) 0.08 (unknown) Lactation period At least 2 years 21 months on average 23 months (range:18-32 months) Calving interval 3 years 2-3 years >2 years Parameters Source Age at sexual 2A maturity 1A 3A 4 5 Age at color 1 change (gray to 2 white) 2 5 Age at 1st 3 conception 3 3 5 Age at senescence 1 5 Pregnancy and 3 birth rates 5 Lifespan 1 2B 3 5 5 Adult annual 3 survival 6.7 8 9 10 Immature annual 2 survival 11 Reproductive rate 2 2 2 3 7 11 12 13 14 15 Lactation period 1 2C 7A Calving interval 1, 2D, 3B 5 7B (1) Brodie, P. F. 1971. A reconsideration of aspects of growth, reproduction, and behavior of the white whale (Delphinapterus leucas) with reference to the Cumberland Sound, Baffi n Island, population. J. Fish. Res. Bd. Can. 28:1309-1318. A. [Canada] Cumberland Sound, Baffi n Island, population, n=124 animals (86% captured in nets which biased the sample toward females with newborns), ASM excluded one immature animal age 15 GLGs, sample sizes not provided though Fig.3 appears to show 51 females in the sample. (2) Sergeant, D. E. 1973. Biology of white whales (Delphinapterus leucas) in western Hudson Bay. J. Fish. Res. Bd Can. 30:1065-1090. A. [Canada] Churchill and Whale Cove in western Hudson Bay, additional information from the Mackenzie Delta, Beaufort Sea and data collected by Khuzin (sample of 33) in the Kara/Barents seas, Russia. B. Found differences in maximum age based on sampling technique. Life span of netted whales tended to be lower (40 GLGs at Whale Cove) than those selected and harpooned (50 GLGs at Churchill, 53 GLGs at Mackenzie Delta). Similar results were reported by Brodie (1971) for whales netted in Cumberland Sound (40 GLGs). C. based on length of gestation (14 months) x 33 lactating/ 22 pregnant whales. D. In 7 of the 29 pregnant females examined from Whale Cove, lactation was still occurring and for some analyses a 2 year calving cycle was assumed for 25% of the adult female population (p. 1084). Concluded "overlap of pregnancy and previous lactation is infrequent so that calving occurs about once in three years. (3) Burns, J. J., and G. A. Seaman. 1986. Investigations of belukha whales in coastal waters of western and northern Alaska. II. Biology and ecology. U.S. Dep. Commer., NOAA, OCSEAP Final Rep. 56(1988):221-357. [Northwest Alaska]; A Sampling occurred in June, a time when most Alaskan belugas are born. It is possible non-pregnant 8-9 GLGs belugas would have conceived before their 10-11 GLGs birth date. B. For some female belugas. This was a tentative conclusion based on high conception rates noted in some females between the ages of 12-13 GLGs and 44-45 GLGs. (4) Robeck, T. R., S. L. Monfort, P. P. Calle, J. L. Dunn, E. Jensen, J. R. Boehm, S. Young, and S. T. Clark. 2005. Reproduction, growth and development in captive beluga (Delphinapterus leucas). Zoo Biol. 24:29-49. [captive belugas]. (5) Suydam, R.S., 2009. Age, growth, reproduction, and movements of beluga whales (Delphinapterus leucas) from the eastern Chukchi Sea. Ph.D. dissert., Univ. Wash., 152 p. (6) Allen, K., and T. Smith. 1978. A note on the relation between pregnancy rate, age at maturity and adult and juvenile mortality rates. Rep. Int. Whal. Comm. 28:477-478. Reviewed in Braham, 1984. (7) Braham, H. W. 1984. Review of reproduction in the white whale, Delphinapterus leucas, narwhal, Monodon monoceros, and Irrawaddy dolphin, Orcaella brevirostris, with comments on stock assessment. Rep. Int. Whal. Comm. (Spec. Issue 6):81-89. A. analysis of data collected by Seaman and Burns (1981). B. based this assumption on data from Brodie (1971) and Sergeant (1973) that age at fi rst pregnancy is 6 years (12 GLGs) and last pregnancy is about 21 years (42 GLGs) resulting in a 14-15 year breeding period, which would allow only 6 calves rather than the 10 calves predicted by the authors if a female's reproductive cycle is three years. However, this calculation was based on 2 GLGs = 1 year, using 42-12 = a 30 year breeding period and a 3-year reproductive cycle would produce 10 calves. (8) Ohsumi, 1979. Interspecies relationships among some biological parameters in cetaceans and estimation of the natural mortality coeffi cient of the southern hemisphere minke whale. Rep. Int. Whal. Comm. 29:397-406. (9) Beland, P., S. De Guise, and R. Plante. 1992. Toxicologie et pathologie des mammiferes marins du Saint-Laurent. INELS, Montreal, QC for the Fond Mondial pour la Nature (Canada), Toronto. Canada] St. Lawrence population. (10) Lesage, V., and M. C. S. Kingsley. 1998. Updated status of the St. Lawrence River population of the beluga, Delphinapterus leucas. Can. Field-Nat. 112(1): 98-114. [Canada] St. Lawrence population. (11) Brodie, P. F., J. L. Parsons, and D. E. Sergeant. 1981. Present status of the white whale (Delphinapterus leucas) in Cumberland Sound, Baffi n Island. Rep. Int. Whal. Comm. 31:579- 582 [Canada] Cumberland Sound, Baffi n Island. (12) Ray, G. C., D. Wartzok, and G. Taylor. 1984. Productivity and behavior of bowheads, Balaena mysticetus, and white whale, Delphinapterus leucas, as determined from remote sensing. Rep. Int. Whal. Comm. (Spec. Issue 6):199-209. (13) Davis, R.A., and K.J.Finley.1979.Distribution, migrations, abundance and stock identity of eastern Arctic white whales.Unpubl. doc.Submitted to Int.Whal.Comm.(SC/31/SM10).[Paper available from http://www.iwcoffi ce.org/].[eastern Arctic]. (14) Davis, R.A., and C.R.Evans.1982.Offshore distribution and numbers of white whales in the eastern Beaufort Sea and Amundsen Gulf, summer 1981.Report for SOHIO Alaska Petroleum Co., Anchorage, and Dome Petroleum Ltd., Calgary, Alberta, by LGL Ltd.76 p.[eastern Beaufort Sea and Amundsen Gulf]. (15) Breton-Provencher, M.1981.Survey of the beluga whale population in the Poste-de-la-Baleine region.Unpubl.Doc.Submitted to Int.Whal. Comm.(SC/32/SM16) [Paper available from http://www.iwcoffi ce.org/]. [Poste-de-la-Baleine region]. Table 3.--Estimated model parameters and their prior distributions. Parameter Description [[lambda].sub.0] Annual growth multiplier (or [e.sup.(intrinsic rate of growth)]) when the population size approaches 0 [[lambda].sub.H0] Annual growth multiplier (or [e.sup.(intrinsic rate of growth)]) for healthy population size approaches 0 [[lambda].sub.k] Annual growth multiplier when the population is at K. [S.sub.0] Survival rate when population size is small (close to 0) and survival is unmodified. [S.sub.k] Survival rate when population size is at K and survival is unmodified. [a.sub.mat] Age of maturity or age at which a female could first give birth [N.sub.1979] Total population size in 1979 K High density population size (derived from aerial survey results in 1979: Calkins [text footnote 1], Hobbs and Shelden [text footnote 5] [K.sub.1999] High density population after 1999 in models M, N and O z Shape parameter resulting in a maximum rate of population increase at MNPL (c.f. Taylor and DeMaster, 1993) [H.sub.t] Total number of whales killed in the hunt in year t (includes landed and struck and lost whales) Pr(Hunt Immature) Probability that an animal taken in the hunt was an immature whale Pr(Hunt Male) Probability that a mature animal taken in the hunt was a male [C.sub.p] Constant mortality effect (killer whale predation) (can be thought of as the number of deaths per carcass) [P.sub.Me] Probability of a catastrophic mortality event occurring in a given year [M.sub.e] Individual probability of mortality during a catastrophic mortality event [P.sub.g] Probability of a group mortality event occurring in a given year [M.sub.g] Individual probability of mortality in the affected group during a group mortality event [G.sub.t] Size of the group in which a group mortality event occurs at time t Parameter Prior distribution [[lambda].sub.0] Healthy model scenario U[1.02, 1.031]; alternative model scenarios U[0.94, 1.031] [[lambda].sub.H0] Healthy model scenario [[lambda].sub.H0] = [[lambda].sub.0]; alternative model scenarios U[1.02, 1.031] or U[[[lambda].sub.0], 1.031] if [[lambda].sub.0] >1.02 [[lambda].sub.k] 1.0 [S.sub.0] U[0.962, 0.974] [S.sub.k] U[[s.sub.0] [[lambda].sub.k]/[[[lambda].sub.0], [s.sub.0]] [a.sub.mat] Fixed at 9 GLGs [N.sub.1979] U[800, K] K U[811,2056] [K.sub.1999] U[100, 800] z Iterative solution for z of [H.sub.t] [H.sub.79-93] U(10, 40) [H.sub.94-12] Reported values (see Table 1) Pr(Hunt Immature) Beta(8, 27) Pr(Hunt Male) Beta(17, 11) [C.sub.p] Fixed at 0, 1, or 2 (depending on model scenario) [P.sub.Me] Fixed at 0.0 or calculated as (1/ E([S.sub.e,t]))/[M.sub.e] (depending on model scenario) [M.sub.e] Fixed at 0.0 or U[0.10, 0.50] (depending on model scenario) [P.sub.g] Fixed at 0.0, or calculated as (1/ E([S.sub.g,t]))/(0.20 [M.sub.e]) (depending on model scenario) [M.sub.g] U[0.10, 1.0] or 0.0 (depending on model scenario) [G.sub.t] Observed distribution of group sizes (with upper limit = [N.sub.t]) Table 4.--Other parameters and variables and their derivation, if applicable. Parameter Description [N.sub.t] Population size at time t [f.sub.a,t], Number of females and males, respectively, of age [m.sub.a,t] a at beginning of year t [f.sub.mat,t], Number of mature females and mature males, [m.sub.mat,t] respectively, at beginning of year t [S.sub.a,t] Probability of an individual at age a in year t surviving to age a+1 year t+1 [b.sub.t] Probability of a mature female giving birth to a live offspring in year t [b.sub.o] Values for fecundity when population size approaches 0 [b.sub.k] Values for fecundity when population size is at K. [C.sub.s] Constant mortality effect annual deaths depends on survival rate for each case [H.sub.j,t] Landings of immature whales in the hunt in year t [H.sub.m,t] Landings of mature male whales in the hunt in year t [H.sub.f,t] Landings of mature female whales in the hunt in year t [NV.sub.t] Total number of immature animals in the vulnerable age classes in year t [S.sub.h,j,t] Survival rate of juveniles from hunting in year t [S.sub.h,f,t] Survival rate of adult females from hunting in year t [S.sub.h,m,t] Survival rate of adult males from hunting in year t [S.sub.t] Population survival rate in year t [S.sub.c] Survival of a constant per capita impact [S.sub.p,t] Survival of a constant mortality effect such as predation by killer whales in year t [S.sub.e,t] Survival of a catastrophe, a stochastic event in year t [S.sub.g,t] survival of a group mortality, a stochastic event affecting a social group in year t [S.sub.x] Product of survival of modifications partitioned in to survival for individual effects [S.sub.c] [S.sub.p,t] E([S.sub.e,t]) and E([S.sub.g,t]). [B.sub.c] Reduced fecundity of a constant per capita impact [[sigma].sub.b] A scaling factor for environmental variation of fecundity [[epsilon].sub.t] A stationary correlated random environmental deviation with mean = 0 and variance = [[sigma].sup.2] and correlation = [rho] [[omega].sub.t] A stationary uncorrelated random environmental deviation with mean = 0 and variance = [[sigma].sup.2] [[sigma].sup.2] Variance of environmental effect on annual survival [rho] Correlation of environmental variation Parameter Derivation [N.sub.t] [summation] [f.sub.a,t] + [m.sub.a,t] [f.sub.a,t], Eq. 1 [m.sub.a,t] [f.sub.mat,t], Eq. 1 [m.sub.mat,t] [S.sub.a,t] Eq 3 [b.sub.t] Eq. 2 [b.sub.o] From Eq. 2, and required to be within interval 0.05 to 0.20 [b.sub.k] From Eq. 2, and required to be within interval 0.00 to [b.sub.o] [C.sub.s] 0.068(1 - [s.sub.346])346 [H.sub.j,t] B([H.sub.t], Pr(Hunt Immature)) [H.sub.m,t] B([H.sub.t] - [H.sub.i,t], Pr(Hunt Male)) [H.sub.f,t] [H.sub.t]-[H.sub.i,t]-[H.sub.m,t] [NV.sub.t] [summation][f.sub.a,t] + [m.sub.a,t] for age 4 to [a.sub.mat] [S.sub.h,j,t] 1-[H.sub.j,t]/[NV.sub.t] [S.sub.h,f,t] 1-[H.sub.f,t]/[f.sub.mat,t] [S.sub.h,m,t] 1-[H.sub.m,t]/[m.sub.mat,t] [S.sub.t] 1.0 or partition of [S.sub.x] [S.sub.c] Eq. 6 [S.sub.p,t] Eq. 7 [S.sub.e,t] Eq. 8; E([S.sub.e,t]) partition of [S.sub.x] [S.sub.g,t] Eq. 9; E([S.sub.g,t]) partition of [S.sub.x] [S.sub.x] Eq. 11 solved iteratively for [s.sub.0][S.sub.x] with [[lambda].sub.o] and [b.sub.o] fixed. [B.sub.c] Eq. 2, [B.sub.c] = 1-[S.sub.0]-1 - [S.sub.0][S.sub.c] [[sigma].sub.b] Eq. 2 [[epsilon].sub.t] Eq. 2 [[omega].sub.t] Eq. 2 [[sigma].sup.2] Eq. 2 [rho] Eq. 2 Table 5.--Model scenarios (Models A-O) for sources of mortality affecting recovery of Cook Inlet beluga whales. "All" indicates all of the additional mortality is attributed to a single mechanism mortality model. "Partition" indicates a random partitioning between two or more mortality models (with an exception when the constant mortality effect is included, its effect is determined by the parameter Cp and the remaining mortality is partitioned). Model Annual growth Per capita Constant ID multiplier survival rate mortality ([lambda].sub.0) and birth rate effect effect ([S.sub.c], ([C.sub.p]) [B.sub.c]) A 1.02,1.031 B 0.94,1.031 All C 1.02, 1.031 1 D 1.02, 1.031 2 1.02, 1.031 1 (year<99) 2(year [greater than or equal to] 99) F 0.94, 1.031 Partition 1 G 0.94, 1.031 H 0.94, 1.031 I 0.94, 1.031 1 J 0.94, 1.031 1 K 0.94, 1.031 Partition 1 L 0.94, 1.031 Partition 1 M 1.02, 1.031 N 1.02, 1.031 1 O 0.94, 1.031 All 1 Model Catastrophic Group mortality Change in Hypothesis ID mortality event percentage carrying ([M.sub.e], ([M.sub.g], capacity [P.sub.Me]) [P.sub.g]) 1979-99 A H1-7 B H1-2 C H3 D H3 H4 F H7 G All H5 H All H6 I Partition H7 J Partition H7 K Partition H7 L Partition H7 M Yes H8 N Yes H8 O Yes H8 Table 6.--Probability of extinction and recovery for Cook Inlet beluga whales and expected growth over the next 20 years based on model scenarios (Models A-O) for sources of mortality affecting this population. Probability of declining is the probability that [N.sub.2H4] < [N.sub.2014]. Probability of Probability of extinction (%) recovery (%) Model by 2064 by 2114 by 2064 by 2114 ID (50 years) (100 years) (50 years) (100 years) A 0 0 81 94 B 0 0 7 14 C 0 2 64 88 D 1 1 50 83 E 0 0 52 85 F 0 6 12 25 G 0 4 26 37 H 9 27 22 34 I 2 13 24 35 J 18 38 21 33 K 2 14 18 30 L 10 29 15 27 M 0 0 0 1 N 0 1 0 1 O 0 3 0 0 Model Mean (SD) of Probability ID annual growth of declining % 2014-34 (%) A 2.1 (1.0) 1 B 0.1 (1.3) 53 C 1.8 (1.2) 2 D 1.6 (1.0) 2 E 1.5 (1.1) 1 F 0.3 (1.3) 42 G 0.5 (2.4) 37 H -0.8 (4.1) 49 I 0.3 (2.4) 41 J -1.1 (4.5) 51 K 0.2 (2.1) 44 L -0.7 (3.6) 54 M 0.1 (0.8) 59 N 0.0 (0.9) 60 O -0.2 (0.9) 71 Table 7.--Statistics for the posterior distributions of the Cook Inlet beluga whale expected growth rate and the Bayes factors for each model scenario compared to Model A (the healthy population scenario) or Model B (the per capita mortality model). PSRF is the potential scale reduction factor. Values of PSRF close to 1.0 indicate that each chain is a representative sample of the posterior distribution. A value of 2(ln(Bayes factor)) greater than 2.0 or less than -2.0 indicates positive evidence of one model over the other. MCMC = Metropolis-Hastings Markov chain Monte Carlo algorithm (King et al., 2010). Model Average PSRF Standard Bayes factor 2 x Ln Bayes ID MCMC jumps error of 2X odds ratio comparison per case Ln Bayes relative to to Model A factor Model A A 3.28 1.44 0.61 1.0E+00 0 B 3.25 1.35 0.87 8.4E+05 27 C 3.28 1.20 0.18 7.3E-01 -1 D 3.26 1.12 0.53 1.4E+03 14 E 3.26 1.24 0.28 4.1E+02 12 F 3.26 1.16 0.42 3.2E+04 21 G 3.27 1.26 0.52 2.7E-03 -12 H 3.28 1.13 0.45 1.1E+02 9 I 3.27 1.11 0.32 4.5E-02 -6 J 3.24 1.23 0.49 6.3E-01 -1 K 3.27 1.40 0.82 6.5E+02 13 L 3.25 1.06 0.22 1.2E+02 10 M 3.23 1.13 0.40 1.1E+09 42 N 3.26 1.26 0.25 5.1E+09 45 O 3.24 1.28 0.22 1.4E+09 42 Model 2 x Ln Bayes ID comparison to Model B A -27 B 0 C -28 D -13 E -15 F -7 G -39 H -18 I -33 J -28 K -14 L -18 M 14 N 17 O 15
Please note: Some tables or figures were omitted from this article.
|Printer friendly Cite/link Email Feedback|
|Author:||Hobbs, Roderick C.; Wade, Paul R.; Shelden, Kim E.W.|
|Publication:||Marine Fisheries Review|
|Date:||Mar 22, 2015|
|Previous Article:||Annual calf indices for Beluga Whales, Delphinapterus leucas, in Cook Inlet, Alaska, 2006-12.|
|Next Article:||Potential natural and anthropogenic impediments to the conservation and recovery of Cook Inlet Beluga Whales, Delphinapterus leucas.|