Ecological Applications, 22(2), 2012, pp. 668–684 Ó 2012 by the Ecological Society of America
A Bayesian hierarchical model of Antarctic fur seal foraging and pup growth related to sea ice and prey abundance LISA M. HIRUKI-RARING,1 JAY M. VER HOEF, PETER L. BOVENG,
JOHN L. BENGTSON
National Marine Mammal Laboratory, Alaska Fisheries Science Center, National Marine Fisheries Service, NOAA, 7600 Sand Point Way NE, Seattle, Washington 98115-6349 USA
Abstract. We created a Bayesian hierarchical model (BHM) to investigate ecosystem relationships between the physical ecosystem (sea ice extent), a prey measure (krill density), predator behaviors (diving and foraging effort of female Antarctic fur seals, Arctocephalus gazella, with pups) and predator characteristics (mass of maternal fur seals and pups). We collected data on Antarctic fur seals from 1987/1988 to 1994/1995 at Seal Island, Antarctica. The BHM allowed us to link together predators and prey into a model that uses all the data efﬁciently and accounts for major sources of uncertainty. Based on the literature, we made hypotheses about the relationships in the model, which we compared with the model outcome after ﬁtting the BHM. For each BHM parameter, we calculated the mean of the posterior density and the 95% credible interval. Our model conﬁrmed others’ ﬁndings that increased sea ice was related to increased krill density. Higher krill density led to reduced dive intensity of maternal fur seals, as measured by dive depth and duration, and to less time spent foraging by maternal fur seals. Heavier maternal fur seals and lower maternal foraging effort resulted in heavier pups at 22 d. No relationship was found between krill density and maternal mass, or between maternal mass and foraging effort on pup growth rates between 22 and 85 days of age. Maternal mass may have reﬂected environmental conditions prior to the pup provisioning season, rather than summer prey densities. Maternal mass and foraging effort were not related to pup growth rates between 22 and 85 d, possibly indicating that food was not limiting, food sources other than krill were being used, or differences occurred before pups reached age 22 d. Key words: Antarctic fur seal; Arctocephalus gazella; Bayesian hierarchical model; diving; ecosystem; ecosystem model; Euphausia superba; foraging; krill; pup growth; sea ice.
INTRODUCTION The western Antarctic Peninsula area has been the focus of attention in recent years, as the marine ecosystem in the Antarctic Peninsula area is one of the most rapidly warming regions on Earth (Clarke et al. 2007, Ducklow et al. 2007). Several large-scale ecosystem analyses have been done around the Antarctic Peninsula and South Georgia areas (e.g., Ducklow et al. 2007, Murphy et al. 2007a) in addition to studies focusing primarily on predator response to climate change and prey availability (Reid and Croxall 2001, Fraser and Hofmann 2003, Reid et al. 2005, Clarke et al. 2007, Trathan et al. 2007), and physical process studies (Murphy et al. 2007b, Stammerjohn et al. 2008). There is an urgent need to assess the response of marine populations to the warming trend, but relationships between environmental measures, prey, and predators are not always straightforward, and studies based on few species over a short period of time may not show the complexity of the ecosystem response (Clarke et al. 2007). Manuscript received 19 January 2011; revised 11 July 2011; accepted 4 August 2011; ﬁnal version received 31 August 2011. Corresponding Editor: J. J. Barber. 1 E-mail: [email protected]
Ecosystem models aim to depict the main components and processes of ecosystems, both to gain a better understanding of interactions between ecosystem elements, and to predict how they change in response to inputs to the model. Ecosystem models often address two aspects: estimating parameters, and predicting missing observations. An estimate of a parameter is a mathematical construct, where the data are related mathematically to give a value (e.g., the slope of a regression line) that gives a better understanding of the relationship between the data points. A prediction takes measurable data (such as time series data) and predicts values for years in which data were not collected. The value of ecosystem models is in their holistic approach to the relationships within the model. Rather than single relationships driving an analysis, each association in turn affects the other interactions within the model. With the emphasis in the past 25 years on ecosystem management (Christensen et al. 1996) rather than management of single species, a single uniﬁed model relating several trophic levels and predator–prey relationships provides information both on how the independent elements relate to one another and on making useful predictions from the linked components of the model. Our goal was to explore ecosystem relationships using such a multilevel uniﬁed, but
FUR SEAL BAYESIAN HIERARCHICAL MODEL
somewhat simpliﬁed, ecosystem model composed of linked linear regressions, accounting for sources of uncertainty as described in Cressie et al. (2009). This approach allows us to test each component regression, explore relationships between components, make inference on relationships, and predict values with uncertainty intervals for years when data are unavailable. In the Antarctic ecosystem, Antarctic fur seals (Arctocephalus gazella) are predators primarily on krill (Euphausia superba, Croxall and Pilcher 1984), which are inﬂuenced by oceanography and sea ice, lower trophic level plankton dynamics, and intermediate predators (Murphy et al. 2007a). Variation in Antarctic fur seal abundance and behavior may be linked to variation in krill distribution and population dynamics and the physical environment, as has been shown for krill-eating penguins in this same system (Fraser and Hofmann 2003). Our goal was to investigate relationships between a physical ecosystem feature (sea ice extent), a prey measure (krill density), predator behaviors (foraging trip durations and diving behavior of female Antarctic fur seals with pups [maternal fur seals]) and predator characteristics (maternal mass and pup growth). To characterize these relationships, we built a single Bayesian hierarchical model, BHM (see overviews by Clark 2007, Royle and Dorazio 2008, Cressie et al. 2009). We chose this model because it accounts for multiple sources of uncertainty, including sampling, process, and estimation variances. Prior to ﬁtting the model, we surveyed the literature and made 10 hypotheses, which were evaluated after ﬁtting the model. Although some of the relationships are obvious, they are still useful to include in the model, as the information in the relationship is used to help with other inferences in the model, such as predicting parameter values or estimating other relationships. Sea ice and krill Hypothesis 1: Increased two-year average of sea ice coverage results in increased krill density.—Annual krill abundance has been shown to be positively correlated with sea ice cover in the previous winter (Brierley et al. 1999a, Hewitt et al. 2003). Long winters with extensive sea ice coverage encouraged early onset of female krill spawning (Siegel and Loeb 1995). Larval krill examined in winters with heavy ice conditions were in better physiological condition (higher lipid content, condition factor, and growth rate) than those examined in winters with light ice conditions (Ross and Quetin 1991, Quetin et al. 1994). Thus, summers following heavy winter ice conditions generally showed increased recruitment of krill via increased spawning success and greater larval survival (Siegel and Loeb 1995). Krill stock size was not always inﬂuenced by recruitment success or failure in one year class, however, because krill reproduce for several years, with females ﬁrst spawning in their third summer (Ross and Quetin 1991, Siegel and Loeb 1995). Ice conditions during the two previous winters might
have a more detectable effect on the regional krill stock size (Siegel and Loeb 1995). Thus, we hypothesized that the two-year average of sea ice extent is positively related to krill abundance. Hypothesis 2: Increased krill density results in increased abundance of gravid female krill.—We considered the proportion of gravid female krill available in the population in addition to general krill abundance. Measurements of prey must take into account both overall quantity and the quality of prey. The diet of several krill predators, including Antarctic fur seals, showed a consistent preference for female krill (Croxall and Pilcher 1984, Hill et al. 1996, Reid et al. 1996, Osman et al. 2004). Antarctic fur seals appeared to actively select female krill at Cape Shireff in the South Shetland Islands (Osman et al. 2004), although in other studies it was unclear whether fur seals preferred female krill or selected large individuals that happen to be female (Hill et al. 1996, Reid et al. 1999). Selection of female krill would provide an energetic advantage to predators, given that female krill had higher lipid content than males, especially gravid females whose ovaries can contain up to 60% of their body lipid (Clarke 1984). Because fur seals appeared to target female krill, a general measure of krill biomass may not accurately reﬂect the proportion of the prey population consumed by fur seals. Gravid krill are found offshore of the general adult population, in deeper water (Trathan et al. 1993, Nicol et al. 2000), and eggs are laid in December through March (Marr 1962). Under certain environmental conditions, female krill can produce continuous batches of eggs throughout the four months of summer (Ross and Quetin 1983). In the Elephant Island area, gravid krill are found in lower densities offshore, compared to higher densities of immature krill in inshore areas (Ichii et al. 1998, Siegel et al. 2002), although gravid krill densities have been linked to overall krill densities (Siegel et al. 2002). We made the simple hypothesis that increased krill abundance results in increased abundance of gravid krill. Krill and maternal fur seal diving intensity Hypothesis 3: Increased krill density results in decreased dive depth intensity; Hypothesis 4: Decreased dive depth intensity results in decreased dive duration intensity.—We used dive depth intensity and dive duration intensity (the sum of dive depths and dive durations, respectively, divided by the number of days with dive data) as measures of the diving effort of maternal fur seals. Simultaneous measurements of prey abundance and predator diving effort are difﬁcult to obtain, especially when prey species are patchy and mobile (Boveng 1996), as is the case with Antarctic krill, the primary prey of Antarctic fur seals (Croxall and Pilcher 1984, Watkins et al. 1986, 1992, Costa et al. 1989, Reid and Arnould 1996). In years of high krill abundance, krill remained at shallower depths during the day than in other years (Godlewska 1996). Dive data from Antarctic
Ecological Applications Vol. 22, No. 2
LISA M. HIRUKI-RARING ET AL.
fur seals at Kerguelen and Bird Islands indicated that when krill were plentiful, dives were shallower (Lea et al. 2006) and more frequent (McCafferty et al. 1998). Therefore, we hypothesized that increased abundance of krill results in a decreased dive depth intensity, which then results in a decreased dive duration intensity. Krill and maternal fur seal mass Hypothesis 5: Increased abundance of gravid krill results in heavier maternal fur seals.—Female Antarctic fur seals and their pups at Kerguelen and Bird Islands had lower mass in years when prey availability was low (McCafferty et al. 1998, Lea et al. 2006). Therefore, we hypothesized that increased abundance of gravid krill results in increased mass of maternal fur seals. Krill and maternal fur seal foraging time Hypothesis 6: Increased abundance of gravid krill results in decreased time spent foraging.—Foraging trip duration of maternal fur seals can vary with prey availability (e.g., Lunn et al. 1993). We decided to use an indirect measure of the mother’s ability to care for her pup, the percentage of time spent away from the pup while foraging, as our measure of maternal foraging effort. We hypothesized that increased gravid krill availability results in a decreased percentage of time spent foraging by maternal fur seals. Maternal fur seal mass and pup growth Hypothesis 7: Heavier females have heavier pups; Hypothesis 8: Heavier females have pups that grow faster.—We considered two measures of pup growth: pup mass early in the summer and pup growth over the austral summer (December through February). Maternal fur seal mass may indicate age of the seal (older seals are heavier; Lunn and Boyd 1993). Older females gave birth to heavier pups at Bird Island (Lunn and Boyd 1993, Lunn et al. 1994). Heavier females had faster growing pups only when food was scarce for both subantarctic fur seals, Arctocephalus tropicalis (Georges and Guinet 2000), and Antarctic fur seals (Lea et al. 2006). Adult female body mass did not affect pup growth when food was plentiful or close to the colony (McCafferty et al. 1998, Guinet et al. 1999, Lea et al. 2006). We hypothesized that heavier females produce heavier pups that grow faster. Although the relationship between maternal mass and pup mass has been clearly demonstrated, we included it in our model so that the effect of the relationship on other elements of the model can be estimated. The information calculated by the model can then be used in the prediction of other parameters. Maternal fur seal foraging time and pup growth Hypothesis 9: Females that forage for a shorter time have heavier pups; Hypothesis 10: Females that forage for a shorter time have pups that grow faster.—Female Antarctic fur seals typically have one pup, which they
raise alone. The duration of the pup-rearing period varies little from year to year (Bonner 1968, Doidge and Croxall 1989). During the four-month pup-rearing period, females alternate foraging trips away from the breeding colony with visits ashore to feed pups (Bonner 1968), providing a direct link between foraging effort and transfer of nutrients to their young. However, the relationship between maternal foraging trip duration, prey abundance, and pup growth is not always clear (e.g., Reid 2002). Seasonal patterns of predator foraging behavior must also be taken into account (Boveng et al. 1998a, Jansen et al. 2002, Biuw et al. 2009), as foraging trip length may be longer at different times during the season. When food availability is low, females spend more time foraging (Lunn et al. 1993, McCafferty et al. 1998, Lea et al. 2006), fewer females have pups, and pups are lighter (Lunn et al. 1994). Conversely, when food availability is high, shorter trip durations are related to heavier birth mass (Lunn et al. 1994). In general, we hypothesized that when maternal fur seals spend less time foraging, pups are heavier and grow faster. METHODS Study site Antarctic fur seals were studied at two adjacent sites, North Annex and North Cove, at Seal Island, South Shetland Islands, Antarctica (60859 0 S, 55823 0 W; Fig. 1) from 1988 to 1995 (ﬁeld seasons were designated by the second year of the split-year austral summer; i.e., 1988 represented the austral summer of 1987/1988). Both colonies were relatively small, with North Cove pup production averaging about 214 pups per year, and North Annex about 67 pups per year (Boveng et al. 1998b). Data collection Adult female fur seals (n ¼ 386) were captured opportunistically, using noose-pole and restraint techniques (Gentry and Holt 1982). In December, females were captured after parturition, but before departing on their ﬁrst foraging trip. Females that were captured in January and February were observed to ensure they were nursing pups before they were captured. Females were weighed and instruments were attached to the seals’ fur on the dorsal midline with quick-setting epoxy glue. Of these females, 165 received both a time–depth recorder (TDR, Mark II or Mark III, Wildlife Computers, Redmond, Washington, USA), and a radio transmitter (Advanced Telemetry Systems, Isanti, Minnesota, USA), while 211 females were instrumented with only a radio transmitter, and 10 were instrumented with a satellite-linked transmitter (Toyocom, Tokyo, Japan) for a separate study (Ichii et al. 2007). TDRs were programmed to record dive depth, dive duration, bottom time, surface interval between dives, temperature, and time. Females alternated foraging trips at sea with visits on land to nurse pups. Each TDR was
FUR SEAL BAYESIAN HIERARCHICAL MODEL
Map of Seal Island and the Antarctic Peninsula.
recovered from the instrumented female at the ﬁrst opportunity after she had completed at least six foraging trips, or after she had lost her pup. Radio-receiving data loggers located at the breeding colony recorded the attendance patterns of instrumented seals. A foraging trip was deﬁned as the period of time an instrumented seal was not being recorded on the data logger, and therefore not at the breeding colony (Walker and Boveng 1995). We conducted daily visual observations to conﬁrm that individuals recorded on the data logger were present at the colony, and that individuals absent from the data logger were not on the beach. Maternal fur seals showed a seasonal nonlinear pattern in foraging behavior: foraging trip duration increased to a peak in late December, decreased until early February, and then remained constant at about one-half the peak duration for the remainder of February (Boveng et al. 1998a, Ichii et al. 2007). Because of this seasonal pattern, dive data (depth and duration) and foraging data (time spent foraging at sea) used in this analysis were restricted to those collected in January. Pup mass was also measured. Each year from 1988 to 1995, fur seal pups (n ¼ 100) were captured and weighed at two-week intervals over a two-month period (late December to February), when they were 22–85 days old. In the BHM, we related a sea ice index, krill density, gravid krill density, diving behavior of maternal fur seals (dive depth and dive duration), maternal foraging time, maternal mass, and pup growth (from 22 to 85 d). However, each variable except the sea ice index was estimated from a sample design; hence, we used data models, or measurement error models (terminology and ideas, including lack of variance identiﬁability, are described in Cressie et al. 2009), which were generally based on classical sampling variances for each variable
prior to its inclusion in the BHM. A description of each variable follows. Sea ice index A sea ice index was created by using data from Hewitt (1997: Table 1), which reﬂected annual sea ice availability, incorporating both the temporal (seasonal duration) and spatial extent of ice cover (106 km2-months; determined by numerical integration of annual curves of monthly sea ice cover measured in 106km2 to get one value per year) on the west side of the Antarctic Peninsula region in the Bellingshausen Sea. We calculated a sea ice index to reﬂect the overwintering success of krill, a long-lived species that spawns for multiple years (Ross and Quetin 1991). For our sea ice index, we added Hewitt’s (1997) ice indices from the current and past year together. For example, our sea ice index for 1986/1987, denoted 1987 in this paper, reﬂected the winters of 1986 and 1985, and therefore we added the ice indices of 1986 (6.77) and 1985 (4.15) from Hewitt (1997) for a value of 10.92 (Table 1). We denote this sea ice index as Si for i ¼ 1, 2, . . . N years. Krill density Yearly measures of krill density (krill/km3) for the study area, reﬂecting the mean density of krill between mid-December and mid-March of the austral summer, were taken from Siegel et al. (1998: Table 1). Krill data were collected from random oblique net tows conducted between mid-December and mid-March in the area between latitudes 608 S and 62830 0 S and longitudes 508 W and 57830 0 W (Siegel et al. 1998), which included the foraging area of female fur seals from Seal Island (Boveng 1996, Ichii et al. 2007). The data were given with standard errors for the sampling estimates. Because
Ecological Applications Vol. 22, No. 2
LISA M. HIRUKI-RARING ET AL.
TABLE 1. Annual estimates and sampling variance (r2) of variables in a Bayesian hierarchical model of the Antarctic ecosystem.
ln(krill density) à
ln(gravid krill) §
Sea ice, Si
1981 1982 1983 1984 1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996
11.72 11.33 10.02 9.02 8.38 8.63 10.92 12.56 9.55 6.32 6.91 9.68 9.6 8.18 8.77 9.86
5.08 5.78 5.62 5.49 4.45
0.285 0.088 0.362 0.168 0.579
3.32 4.38 2.73 1.87 2.99 3.37 3.50 4.29 4.79
0.041 0.156 0.051 0.067 0.171 0.066 0.119 0.562 0.199
2.925 2.442 1.398 1.645 2.138 1.647 2.514 4.257 4.773
0.049 0.257 0.078 0.070 0.179 0.098 0.130 0.562 0.199
Dive depth intensity (m/d)
Dive duration intensity (min/d)
2443 2589 4862 5248 4320 2666 4312 3214
66 853 57 941 94 730 32 340 37 029 82 641 181 119 19 180
145.9 108.5 221.8 259.3 196.3 146.0 205.3 200.5
288.3 263.4 149.8 59.5 131.3 148.1 197.4 122.2
Notes: See Methods for a description of each variable and how it was calculated. Units for variables are indicated in separate footnotes. Intensity of dive depth and duration, adult female mass, and percentage of time in January spent foraging all refer to maternal fur seals. Sea ice index (106 km2-months). à Log-transformed krill density, originally measured as krill/km3. § Log-estimated gravid krill density, originally measured as gravid krill/km3. } Estimated growth intercepts (mass at age 22 d) for female and male pups from exploratory regression ﬁt.
the data were skewed, we used a log transformation. If we denote krill density as Xi and the log of krill density as Ki ¼ ln(Xi ), then using the delta method (Dorfman 1938) we obtain an approximate variance: r2Ki ¼ varðKi Þ ’ðXi Þ2 varðXi Þ which is given in Table 1. Gravid krill In addition to krill density, Loeb et al. (1997: Table 1) reported the proportion of mature female krill in advanced reproductive stages during January–February in a region near our study area, along with yearly sample sizes ni. Let this be denoted Yi for the ith year. No variances or standard errors were given. However, the variance of a mean proportion from a series of samples from a binomial distribution with proportion p is p(1 p)/n, so we used as an approximate variance var(Yi ) ¼ Yi(1 Yi )/ni. We computed the expected density of gravid krill as XiYi. An unbiased estimator of the variance of a product is given by Goodman (1960) as: varðXi Yi Þ ¼ Xi2 varðYi Þ þ Yi2 varðXi Þ varðXi ÞvarðYi Þ: The estimated density of gravid krill was skewed, with a few very large values. A log transformation was used to normalize this variable: Gi ¼ ln(XiYi ). The delta method (Dorfman 1938) was used to approximate this variance: 2
r2Gi ¼ varðGi Þ ’ðXi Yi Þ varðXi Yi Þ which is given in Table 1.
Maternal fur seal diving intensity For each maternal fur seal with a time–depth recorder in each year (total n ¼ 165), dive depth intensity and dive duration intensity were calculated for the month of January. Dive depth intensity and dive duration intensity were obtained for each maternal fur seal by summing all dive depths and dive durations in January, and dividing by the number of days with dive data. Calculating these values allowed us to assess the daily diving effort expended by the seals, taking into account the number of dives as well as the depth and duration of the dives, rather than simply calculating average dive depth and dive duration. As previously described, only data from January were used, for consistency with maternal fur seal foraging data. Using diving data from January also allowed us to relate diving behavior to krill indices from January–February. An annual estimate of the average dive depth intensity, Pi, was computed by averaging the individual dive depth intensities for all maternal fur seals with time–depth recorders in each year. We also computed the variance, r2P;i , for this mean (Table 1). Similarly, we computed an annual estimate of the average dive duration intensity, Di, by averaging the individual dive duration intensities for all maternal fur seals with time–depth recorders in each year, and computed the variance, r2D;i , of this mean (Table 1). Maternal fur seal foraging time For each maternal fur seal with a radio transmitter in each year (total n ¼ 376), the percentage of time spent foraging (for brevity, referred to as ‘‘Foraging (%)’’ in
FUR SEAL BAYESIAN HIERARCHICAL MODEL
TABLE 1. Extended.
Adult female mass (kg)
Pup growth intercept (kg) }
Female, gˆ i,1
Male, gˆ i,2
66.34 67.97 67.94 74.63 73.88 70.00 72.58 77.38 65.55
3.30 6.03 3.02 0.65 0.53 0.70 0.59 0.52 2.88
39.85 39.12 39.64 40.47 43.37 47.45 43.14 46.36 49.34
1.07 1.31 1.10 0.84 0.87 1.10 0.95 1.27 0.84
5.71 5.24 4.85 5.38 5.84 5.13 5.88 6.11
6.47 5.38 5.56 5.73 6.43 5.56 5.79 6.30
Table 1) was computed for the month of January. An annual estimate of the percentage of time seals spent foraging in January, Fi, was computed by averaging across individuals for each year. We also computed the variance, r2F;i , for that same mean (Table 1). In January, maternal fur seals at Seal Island ranged an average maximum distance of 52 km from the island, with an average trip length of 2–3 d (Ichii et al. 2007). Fur seals took longer trips in the early part of the breeding season, with an average trip length of 5 d (maximum distance of 99 km) in December (Ichii et al. 2007). Because the length of foraging trips typically reﬂects the distance traveled by the seals to forage, the percentage of time spent away from the pup should be a reasonable representation of the effort spent foraging by maternal fur seals, incorporating both travel and foraging time. The percentage of time spent away from the pup also allows for comparisons of foraging effort regardless of differences in the foraging trip length of individual females. A female that takes less frequent, longer foraging trips can be compared directly with another female that takes more frequent, shorter trips. Maternal fur seal mass Maternal fur seals captured for instrumentation were weighed after capture. The mass of females weighed in December did not differ signiﬁcantly from the mass of those weighed in January and February of the same years (1988, P ¼ 0.76; 1990, P ¼ 0.27; 1991, P ¼ 0.78); therefore, all were pooled. An annual estimate of maternal mass, Ai, was computed by averaging among individuals for each year. We also computed the variance, r2A;i , of the mean (Table 1). Exploratory analysis and model selection Information criteria methods (e.g., Akaike information criterion, AIC, and Bayesian information criterion,
BIC) are useful for simpler models but become more difﬁcult as models become more complex. For example, we use BIC in an exploratory way to identify a subcomponent model within the BHM. As models become more complex, it is often difﬁcult to compute the likelihood, difﬁcult to count parameters (what is ﬁxed and what is random?), and it may require a long time to ﬁt a model. Therefore, we used the background literature given in the Introduction and exploratory analysis to deﬁne our model prior to ﬁtting it. We had within-season pup masses across years, by location and by sex, and we needed to identify how the factors interacted to create a submodel of pup growth for the BHM. Fig. 2 plots daily average mass as a function of age for each sex with a lowess-smoothed curve for each year. Although growth was expected to show some logistic-type pattern, it appeared that our data were closely approximated by a linear model over the range of 22–85 d old. Hence, we used linear models to ﬁt pup growth over the range of 22–85 d. We explored factors that affect pup growth by ﬁtting linear models with effects for location, sex, and year, in addition to age. First, to assess the effect of location (colonies at North Cove and North Annex), we ﬁt a full three-way interaction for intercepts (sex 3 year 3 location) and four-way interaction for slope (sex 3 year 3 location 3 age) (Model 1, Table 2). We then removed location from both intercept and slope (Model 2, Table 2). BIC indicated that location was not an important factor. Next, we tried dropping various combinations of both sex and year for both intercepts and as interactions with age (Models 3–6, Table 2). The most parsimonious model was Model 5, with separate intercepts for each sex/year combination, and a separate slope on age for each sex (Table 2). Next, we tried to make Model 5 simpler by dropping some sex and age combinations (Models 7–9, Table 2), but Model 5 remained the bestﬁtting model, and it was this form that was used in the BHM. For the BHM, we wished to create a structure reﬂecting the dependence among all of the variables, and we took a graphical exploratory approach. We examined scatterplots between all annual variables in Table 1 (Fig. 3), in addition to the ﬁtted intercepts for each sex (i.e., male and female pup mass at 22 d) in Model 5 of Table 2. Fig. 3 shows that there could be several important relationships among variables. Based on the relationships established in the Introduction and the scatterplots, we wished to estimate relationships shown by the arrows among plots in Fig. 3. Notice that some variables were initially response variables in some relationships, but then became predictor variables in other relationships. This forms one of the characteristics of a hierarchical model. All dependencies are shown as a directed graph in Fig. 4, where the þ and symbols show the relationships that we might expect, given the introductory literature and scatterplots in Fig. 3. We did not attempt any formal modeling using linear regression.
Ecological Applications Vol. 22, No. 2
LISA M. HIRUKI-RARING ET AL.
FIG. 2. Lowess ﬁts of daily average Antarctic fur seal pup mass for male and female pups at Seal Island, Antarctica, beginning at 22 days of age. A separate curve (solid line) was ﬁt for each year, 1988–1995.
There were many sources of uncertainty. The usual regression techniques assume that the predictor variable was observed without error. Moreover, because the response variable in one regression model became the predictor variable in another, we wanted a single model that linked all of the variables and accounted for all sources of uncertainty. For this, we used a BHM.
½Gi j lG;i ; r2G;i ¼ N ðlG;i ; r2G;i Þ ½lG;i j bG0 ; bG1 ; lK;i ; d2G ¼ N ðbG0 þ bG1 lK;i ; d2G Þ
Based on our exploratory analysis and literature review, we set up a BHM following the ideas and notation of Cressie et al. (2009). In the following 14 equations used to create the BHM, [A j B] means ‘‘the distribution of A given B’’ and N(l, r2) is a normal distribution with mean l and variance r2. The notation for the observed annual data Si, Ki, Pi, Di, Gi, Fi, Ai was given earlier in Methods and is found in Table 1, along with estimated sampling variances r2K;i , r2P;i , . . . , r2A;i , which we now assume to be known and ﬁxed. Wijk is the mass for the kth pup of the jth sex in the ith year and xijk is the age of that pup in days: ½Ki j lK;i ; r2K;i ¼ N ðlK;i ; r2K;i Þ
½Fi j lF;i ; r2F;i ¼ N ðlF;i ; r2F;i Þ
¼ N ðbK0 þ bK1 Si ;
½lP;i j bP0 ; bP1 ; lK;i ;
¼ N ðbP0 þ bP1 lK;i ;
½Di j lD;i ; r2D;i ¼ N ðlD;i ; r2D;i Þ ½lD;i j bD0 ; bD1 ; lP;i ; d2D ¼ N ðbD0 þ bD1 lP;i ; d2D Þ
ð4Þ ð5Þ ð6Þ
½lA;i j bA0 ; bA1 ; lG;i ; d2A ¼ N ðbA0 þ bA1 lG;i ; d2A Þ ð12Þ ½Wijk j gij ; sj ; xijk ; r2W ¼ N ðgij þ sj xijk ; r2W Þ
TABLE 2. Bayesian Information Criterion (BIC) values for linear models of Antarctic fur seal pup mass with effects for location, sex, and year, in addition to age.
½Pi j lP;i ; r2P;i ¼ N ðlP;i ; r2P;i Þ d2P
½Ai j lA;i ; r2A;i ¼ N ðlA;i ; r2A;i Þ
½lK;i j bK0 ; bK1 ; Si ;
½lF;i j bF0 ; bF1 ; lG;i ; d2F ¼ N ðbF0 þ bF1 lG;i ; d2F Þ ð10Þ
Bayesian hierarchical model
Sex 3 year 3 loc þ sex 3 year 3 loc 3 age Sex 3 year þ sex 3 year 3 age Sex 3 year þ year 3 age Year þ sex 3 year 3 age Sex 3 year þ sex 3 age Sex þ sex 3 year 3 age Sex 3 year þ age Year þ sex 3 age Sex þ sex 3 age
Model No. number parameters
3 4 5 6 7 8 9
25 25 19 19 18 11 5
12 784.93 12 789.88 12 669.87 12 737.13 12 760.50 12 729.76 12 741.57
The lowest BIC value is associated with model 5; this was the model used in the Bayesian hierarchical model.
FUR SEAL BAYESIAN HIERARCHICAL MODEL
Fig. 3. All pairwise scatterplots of variables in Table 1. The arrows indicate the direction of assumed dependence. ‘‘Forage percentage’’ refers to the percentage of time spent foraging by maternal fur seals in January. Units for other plots are as follows: pup mass at age 22 d for females and males, kg; maternal mass, kg; gravid krill, ln(gravid krill density), originally measured as krill/ km3; dive duration, min/d; dive depth, m/d; krill density, ln(krill density), originally measured as krill/km3; sea ice index, 106 km2months. ½gij j c0; j ; cF; j ; cA; j ; lF;i ; lA;i d2g ¼ N ðc0; j þ cF; j lF;i þ cA; j lA;i ; d2g Þ:
Eqs. 1 and 2 form a measurement model between krill and sea ice, where krill is measured with error and depends linearly on sea ice. Both Eqs. 1 and 2 create a BHM that can be used to combine sampling and process variance (Morris 1983, Ver Hoef 1996). Now, suppose that there was no sampling variance; that is, we could observe the krill index exactly each year. We would not
expect a perfect linear relationship between sea ice and krill density. Thus, we use a second source of variation to ‘‘absorb’’ the lack of fit between the true krill density and sea ice. The fit is the linear model bK0 þ bK1Si in Eq. 2, and the lack-of-fit variance is denoted d2K . This lack of fit is often due to unmeasured ecological processes, so we call this a ‘‘process’’ variance. For example, krill density can be affected by its own food source, predators, disease, and many other factors. Obviously, a linear relationship to sea ice is not a complete model, and d2K replaces knowledge of other ecological processes,
LISA M. HIRUKI-RARING ET AL.
Ecological Applications Vol. 22, No. 2
FIG. 4. Structure of Bayesian hierarchical model relating sea ice, krill density indices, maternal Antarctic fur seal diving parameters (dive depth intensity, dive duration intensity), foraging effort, and mass, and Antarctic fur seal pup growth. A ‘‘plus’’ (þ) indicates a hypothesized positive relationship, and a ‘‘minus’’ () indicates a hypothesized negative relationship. ‘‘Forage percentage’’ refers to the percentage of time spent foraging by maternal fur seals in January.
which justifies calling it a process variance. Also recall that we assume r2K;i to be fixed, which allows d2K to be estimated (if r2K;i were free, then we could not separate the two sources of variability, known as an ‘‘identifiability’’ issue in statistics). Note that because r2K;i itself is estimated, d2K may also absorb sample error variation not captured in r2K;i . Eqs. 3 and 4 then form a measurement error model where dive depth intensity, again measured with error, depends hierarchically on krill density, which is also uncertain. The same discussion regarding process variance and identiﬁability for Eqs. 1 and 2 applies here, as well as for the following sets of equations. Eqs. 5 and 6 form a measurement error model where dive duration intensity, measured with error, depends hierarchically on dive depth intensity, which is also uncertain. Eqs. 7 and 8 form a measurement error model where the proportion of gravid krill, measured with error, depends hierarchically on krill density, which is also uncertain. Eqs. 9 and 10 form a measurement error model where the percentage of time spent foraging by maternal fur seals, measured with error, depends hierarchically on the proportion of gravid krill, and Eqs. 11 and 12 form a measurement error model where maternal mass, measured with error, also depends hierarchically on the proportion of gravid krill. Eq. 13 is a linear growth model for both sexes of pups, with different intercepts and growth rates for each, and then Eq. 14 makes the intercepts of Eq. 13 depend linearly on the percentage of time spent foraging and maternal mass. Hence, Eqs. 1–14 establish the hierarchical model shown in Fig. 4. We used the following vague priors: ½bK0 ¼ ½bK1 ¼ ½bP0 ¼ ½bP1 ¼ . . . ¼ ½c0;1 ¼ ½cF;1 ¼ . . . ¼ ½cA;2 ¼ ½s1 ¼ ½s2 ¼ N ð0; 1:0 3 106 Þ
and ½d2K ¼ ½d2P ¼ . . . ¼ ½d2A ¼ ½d2g ¼ ½r2W ¼ IGð0:001; 0:001Þ
where IG(, ) is an inverse gamma distribution. The joint distribution of all parameters and data, given sea ice, days that pup mass was measured, and the sampling variances, is obtained by taking the product of all distributions given in brackets in Eqs. 1–16: ½ j fSi ; xijk ; r2K;i ; :::; r2A;i g ¼ ½Ki j lK;i ; r2K;i ½lK;i j bK0 ; BK1 ; Si ; d2K :::½bK0 :::½r2W for i ¼ 1, . . . , N, j ¼ 1, 2, and k ¼ 1, . . . , nij, where nij is the number of pups weighed of the jth sex in the ith year. A Markov chain Monte Carlo (MCMC) routine was used to sample from the posterior distribution of the parameters, given the data (Link et al. 2002, Cressie et al. 2009): ½h j fSi ; Ki ; :::; Ai ; Wijk ; xijk ; r2K;i ; :::; r2A;i ; i ¼ 1; :::; N; j ¼ 1; 2; k ¼ 1; 2; :::; Kij g where h is the vector h ¼ ð lK;i ; bK0 ; bK1 ; lP;i ; bP0 ; bP1 ; :::; c0;1 ; cF;1 ; :::; cA;2 ; gij ; s1 ; s2 ; d2K ; :::; d2g ; r2W Þ for i ¼ 1, 2, . . . , N and j ¼ 1, 2. We used the software WinBUGS to estimate all parameters. We used a burnin of 10 000 iterations, and then used 100 000 samples to estimate the posterior distributions of all parameters and any functions of the parameters. We examined the chains of all parameters and judged them to be converged and mixing well. Because our primary focus was on the highest level of the hierarchy, pup growth, we considered one variation of the model above to examine whether the percentage of time spent foraging, Fi, and maternal mass, Ai, affected the pup growth slope parameter (as well as intercepts), even though our exploratory analysis indicated no relationship. For that model, Eq. 13 is modiﬁed to ½Wijk j gij ; sj ; xijk ; r2W ¼ N ðgij þ sij xijk ; r2W Þ
FUR SEAL BAYESIAN HIERARCHICAL MODEL
FIG. 5. Relationships between variables from Table 1. See Results for explanation of symbols, lines, curves, and ellipses. Posterior densities (unitless) for slope parameters showing the relationship between the variables in each graph are shown as insets in each panel. ‘‘Forage percentage’’ refers to the percentage of time spent foraging by maternal fur seals in January; ln-transformed densities of all krill and of gravid krill were originally measured as krill/km3.
and we added the following:
Bayesian hierarchical model
½sij j n0; j ; nF; j ; nA; j ; lF;i ; lA;i d2n ¼ N ðn0; j þ nF; j lF;i þ nA; j lA;i ; d2n Þ where the similarity to Eq. 14 should be clear.
In Fig. 5A, the estimated log-transformed krill density is shown by the open squares. The black vertical lines are the 95% conﬁdence intervals derived from the
Ecological Applications Vol. 22, No. 2
LISA M. HIRUKI-RARING ET AL.
TABLE 3. Ecosystem relationships examined in a Bayesian hierarchical model. Linear model
Krill density and sea ice Krill density and dive depth intensity Dive duration intensity and dive depth intensity Krill density and gravid krill abundance Gravid krill density and average percentage of time spent foraging in January Gravid krill density and annual average maternal mas Pup mass and age Female Male
Eq. Eq. Eq. Eq. Eq.
Fig. Fig. Fig. Fig. Fig.
Average percentage time spent foraging in January and pup mass at 22 d Female Male
Average maternal mass and pup mass at 22 d Female Male
Average percentage time spent foraging in January and pup growth rate, 22–85 d
2 4 6 8 10
Eq. 12 Eq. 13
Female Male Average maternal mass and pup growth rate, 22–85 d
95% credible interval
5A 5B 5C 5D 5E
bK1 bP1 bD1 bG1 bF1
0.4449 1712 0.0424 1.178 4.276
(0.0204, (2772, (0.0277, (0.756, (9.042,
Fig. 6A Fig. 6D
(0.0982, 0.1071) (0.1290, 0.1374)
Fig. 6B Fig. 6E
(0.1190, 0.0339) (0.1447, 0.0064)
Fig. 6C Fig. 6F
(0.007103, 0.1719) (0.032879, 0.1320)
nˆ F,1 nˆ F,2
(0.0039, 0.0043) (0.0065, 0.0054)
nˆ W,1 nˆ W,2
(0.0065, 0.0042) (0.0064, 0.0052)
0.7100) 627.9) 0.0601) 1.647) 0.581)
Notes: Numbered linear models refer to equations in Methods. The ﬁgure references show the model ﬁt and the posterior density of the parameter.
sampling error (the estimated values and variances are given as Ki and r2K;i , respectively, in Table 1). The ﬁtted linear model is shown in Fig. 5A as the dashed line, and the 95% credibility region for it is shown by the outer dashed curves. The inset in Fig. 5A shows the posterior density for the slope parameter bK1. The posterior mean for bK1 (0.4449 ln (krill/km3) per 106km2-months, 95% credible interval (CI) ¼ (0.0204, 0.7100), Table 3) implies that there is good evidence that krill are more abundant when there is more sea ice in the previous two years. In Eq. 1, the true krill density is denoted lK,i. In Fig. 5A, the posterior means for lK,i are shown as solid black circles, and the 95% CIs are shown with gray vertical bars. Note that in some years we had sea ice values but not krill density. The lK,i values can still be predicted, but the CIs are much wider, which is evident for two years in Fig. 5A. Fig. 5B shows part of the BHM model between krill density and annual dive depth intensity in January. The ﬁtted linear model is shown in Fig. 5B as the dashed line, and the 95% credibility region for it is shown by the outer dashed curves. The inset in Fig. 5B shows the posterior density for the slope parameter bP1 in Eq. 4. The posterior mean for bP1 (1712 m/d per ln (krill/ km3), 95% CI ¼ (2772, 627.9); Table 3) implies that there is strong evidence that maternal fur seals make shallower dives in January, when krill are more abundant. Note that there are just eight years with observations on both krill density and annual dive
depth, yet we can estimate values for all years, some of which are outside the graph. As in Fig. 5A, the estimates of the true dive depth intensities are taken from the mean of the posterior distribution for the lP,i values of Eqs. 3 and 4, which are shown as solid circles in Fig. 5B. The 95% CIs in Fig. 5A were given as gray vertical lines. In Fig. 5B there is uncertainty in the relationship between lK,i and lP,i on both the x- and y-axes. We formed 95% credibility ellipses (shown in gray) from a normal approximation to the posterior using the MCMC samples of bivariate pairs of estimates for lK,i and lP,i. As before, large ellipses are seen when predicting lP,i in years without observations, and the ellipses are smaller in years with dive depth data. The 95% credibility ellipses for the estimated values in Fig. 5B that are outside of the graph are visible as partial ellipses. The response variable in Fig. 5B becomes the predictor variable in Fig. 5C, which shows the model for dive duration intensity as a function of dive depth intensity. The features of Fig. 5C are similar to Fig. 5B. The ﬁtted linear model in Fig. 5C shows a positive relationship between dive depth and duration intensities. The posterior density for the slope parameter bD1 in Eq. 6 had a mean of 0.0424 min/d per m/d (95% CI ¼ (0.0277, 0.0601); Table 3), probably reﬂecting a tendency for maternal fur seals to make shorter dives when they make shallower dives in January.
FUR SEAL BAYESIAN HIERARCHICAL MODEL
FIG. 6. Relationships between Antarctic fur seal pup mass, percentage of time spent foraging by maternal fur seals in January, and maternal fur seal mass. Female and male Antarctic fur seal pup mass at age 22 d for 1981–1987 and 1996 are estimated from the Bayesian hierarchical model. See Results for explanation of symbols, lines, curves, and ellipses. Posterior densities for slope parameters showing the relationship between pup mass and maternal foraging, and between pup mass and maternal mass, are shown as insets in the respective panels.
Fig. 5D uses krill density as a predictor variable, as in Fig. 5B, but now the response variable is gravid krill. Of course, gravid krill was computed, in part, directly from krill density, so we expected this to be a positive relationship. The ﬁtted linear model in Fig. 5D shows a positive relationship between krill and gravid krill densities. The posterior density for the slope parameter bG1 in Eq. 8 had a mean of 1.178 ln (krill/km3) per ln (krill/km3) (95% CI ¼ (0.756, 1.647); Table 3). This implies that there is strong evidence that gravid krill density is higher when krill density is higher. The response variable in Fig. 5D becomes the predictor variable in Fig. 5E, which models the annual average percentage of time that maternal fur seals spent away from their pups in January (for brevity, referred to as ‘‘forage percentage’’ in Fig. 5E) as a function of gravid krill density. The ﬁtted linear model in Fig. 5E shows a negative relationship between the percentage of time spent foraging and gravid krill densities. The posterior density for the slope parameter bF1 in Eq. 10 had a mean of 4.276% time per ln (krill/km3) (95% CI ¼ (9.042, 0.581); Table 3). This implies that there is good evidence that fur seals spend less time foraging in January in years with abundant gravid krill.
The response variable in Fig. 5D again becomes the predictor variable in Fig. 5F, which models annual average maternal mass as a function of gravid krill density. The ﬁtted linear model in Fig. 5F shows no relationship between maternal mass and gravid krill densities. The posterior density for the slope parameter bA1 in Eq. 12 had a mean of 0.245kg per ln (krill/km3) (95% CI ¼ (5.213, 3.872); Table 3). Thus there is no evidence that maternal fur seals weigh more in years with abundant gravid krill. The ﬁnal link in the BHM related annual average percentage of time spent foraging and maternal mass to pup mass at 22 d (i.e., the intercepts (Eq. 14) of the pup growth model (Eq. 13)). The annual posterior means for gi, j, along with 95% CIs, are shown in Fig. 6A, D. There appears to be a trend of higher pup mass at 22 days of age that began to decline in the mid-1980s, and then a possible recovery in the mid-1990s. The prediction intervals for unsampled years are very high, however. In all years, male pup mass values are clearly higher than female pup mass values. The ﬁtted linear models in Fig. 6B, E show negative relationships between annual average percentage of time spent foraging and pup growth intercepts. The posterior density for the slope parameter cF,1 (females) in Eq. 14
Ecological Applications Vol. 22, No. 2
LISA M. HIRUKI-RARING ET AL.
TABLE 4. Hypothesized relationships between model parameters compared with Bayesian hierarchical model (BHM) results. Model parameters Sea ice and krill density
Krill density and maternal fur seal diving intensity
Hypothesized relationship 1) Increased sea ice results in increased krill density. 2) Increased krill density results in increased gravid krill. 3) Increased krill density results in decreased dive depth intensity. 4) Decreased dive depth intensity results in decreased dive duration intensity.
BHM result supports hypothesis?
Krill density and maternal mass
5) Increased gravid krill results in heavier adult females.
Krill density and time spent foraging by female fur seals
6) Increased gravid krill results in decreased foraging time.
Maternal fur seal mass and pup growth
7) Heavier females have heavier pups. 8) Heavier females have pups that grow faster.
yes (both male and female pups) no relationship between female mass and pup growth after 22 d
Time spent foraging by female fur seals and pup growth
9) Decreased female foraging time results in heavier pups. 10) Decreased female foraging time results in pups that grow faster.
yes (both male and female pups)
had a mean of 0.0420 kg per percentage of time (95% CI ¼ (0.1190, 0.0339)), and the posterior density for the slope parameter cF,2 (males) in Eq. 14 had a mean of 0.0656 kg per percentage of time (95% CI ¼ (0.1447, 0.0064); Table 3). Thus, there is reasonable evidence that pup mass at 22 days of age is higher in years when maternal fur seals spend less time foraging. The ﬁtted linear models in Fig. 6C, F show positive relationships between annual average maternal mass and pup growth intercepts. The posterior density for the slope parameter cA,1 (females) in Eq. 14 had a mean of 0.09639 kg per kg (95% CI ¼ (0.007103, 01719)), and the posterior density for the slope parameter cA,2 (males) in Eq. 14 had a mean of 0.05549 kg per kg (95% CI ¼ (0.032879, 0.1320); Table 3). Thus there is good evidence that pup mass at 22 days of age is higher in years when maternal fur seal mass is higher. The intercepts of the pup growth models varied by sex and year (Fig. 6A, D), whereas the slopes of the pup growth models varied only by sex (Eq. 13; Fig. 2, Table 2). The posterior density for the slope parameter s1 (females) in Eq. 13 had a mean of 0.1025 kg/d (95% CI ¼ (0.0982, 0.1071)), and the posterior density for the slope parameter s2 (males) in Eq. 13 had a mean of 0.1331 kg/ d (95% CI ¼ (0.1290, 0.1374); Table 3). The pup growth models looked similar to Fig. 2. Pup growth models showed that growth rate did not vary by year (Table 2, Fig. 2). However, to investigate whether there could be a relationship between pup growth, maternal foraging behavior, and mass, we ﬁt another BHM in which we allowed the slope of pup mass (i.e., growth rate) for both male and female pups to depend on maternal fur seal mass and the percentage of time spent foraging by maternal fur seals in January. The posterior distribution of the estimated relationships
no relationship between foraging time and pup growth after 22 d
indicated that pup growth rate between 22 and 85 days old did not have a relationship with the maternal percentage of time spent foraging (nˆ F,1 ¼ 0.00074 kg/d per percentage time for female pups, 95% CI ¼ (0.0039, 0.0043); nˆ F,2 ¼ 0.00048 kg/d per percentage time for male pups, 95% CI ¼ (0.0065, 0.0054)) or with maternal female mass (nˆ W,1 ¼ 0.00072 kg/d per kg for female pups, 95% CI ¼ (0.0065, 0.0042); nˆ W,2 ¼ 0.0014 kg/d per kg for male pups, 95% CI ¼ (0.0064, 0.0052); Table 3). The relationships estimated in the BHM generally matched the hypotheses from our literature survey; seven of 10 hypotheses were supported by the BHM results (Table 4). The three hypotheses that were not matched dealt with krill density and maternal fur seal mass, maternal fur seal mass and pup growth after 22 d, and percentage of time spent foraging and pup growth after 22 d. DISCUSSION The BHM framework allowed us to estimate relationships between ecosystem features, with each relationship contributing data to the next level of analysis. We established relationships between ecosystem features and fur seal behavior and physiology, while ﬁnding that some of the anticipated associations were not veriﬁed. The model allowed us to make predictions (for example, pup mass at 22 days; Fig. 6A, D), with uncertainty intervals, for years when data were not collected. In general, the relationships estimated by the model matched our hypotheses (Table 4). The three relationships that did not match our hypotheses involved either maternal fur seal mass or pup growth rates between 22 and 85 days old. Maternal mass had no relationship with gravid krill availability, or with pup growth rate after 22
FUR SEAL BAYESIAN HIERARCHICAL MODEL
days of age. Because most of the females were weighed when captured at the beginning of the austral summer (December), maternal fur seal mass may be more a function of environmental conditions or prey availability over the previous winter (Lunn and Boyd 1993), rather than a function of summer (January) prey densities. We did not have a measure of winter krill density to relate to maternal fur seal mass, so we could not verify this relationship. Because Antarctic fur seals feed during the nursing period and do not rely solely on energy stores to feed their pups, the mass of the mother may be less of a factor for pup growth rates after 22 days (Lunn et al. 1994). Maternal fur seal mass still conferred an advantage to pups, however, as heavier females did have heavier pups at 22 days (Fig. 6C, F). Our data support ﬁndings that heavier females give their energetic advantage to their young early on, either with higher birth mass (Lunn et al. 1994) or higher growth rates prior to 22 days of age. The lack of relationship between pup growth rate after 22 days, maternal fur seal mass, and time spent foraging by females (Table 4) could simply reﬂect that food was abundant at Seal Island. Smaller maternal size and longer foraging trip duration of maternal females negatively inﬂuenced pup growth rate during the nursing period when food was scarce (Lunn et al. 1993, Georges and Guinet 2000, Lea et al. 2006). When food was readily available, maternal size and foraging trip duration had no effect on pup growth (McCafferty et al. 1998, Guinet et al. 1999, Reid 2002, Lea et al. 2006). Prey ﬂuctuations of the magnitude seen at Bird Island (e.g., Lunn et al. 1993, Brierley and Watkins 1996, Brierley et al. 1999b) have not been documented at our study site, Seal Island. The Elephant Island area (containing Seal Island; Fig. 1) and South Georgia (containing Bird Island) had the same general trends in krill abundance (i.e., highs and lows in the same years at both locations; Brierley et al. 1999a), and length– frequency distribution (Reid et al. 2002). However, the mortality rate of krill at South Georgia is much higher than in the Elephant Island area, which may account for the greater variability in krill biomass at South Georgia (Reid et al. 2002). Pup growth measurement can also be affected by the methodology used to sample pups, as well as mortality of pups over the summer weighing period (Lunn et al. 1993, Reid 2002). Cross-sectional sampling of pups may mask the true growth rate of pups, because it is restricted to the population available to be sampled, which can be biased by different factors (e.g., geographical distribution of pups to be sampled or mortality of pups during the sampling period; Trites 1991, 1993, Reid 2002). In our study, virtually all pups in our small colonies were sampled (Boveng et al. 1998b); thus, spatial factors related to sampling probably did not affect the calculation of our growth rates. At Seal Island, mortality of pups over the study period was greater at North Cove (pup counts declined
by .50% of peak values each year) than at North Annex (counts declined by ,20% of peak values; Boveng et al. 1998b). North Cove pups were subject to leopard seal (Hydrurga leptonyx) predation, which was responsible for the removal of an average of 34.4% of pups annually from 1991 to 1995 (Boveng et al. 1998b, Hiruki et al. 1999). However, pup growth rates did not differ between the two colonies, implying that leopard seals did not prey preferentially on one size class of pups (i.e., larger, more active pups or smaller, slower-swimming pups) and thus did not affect the population of pups available for sampling. Data improvements Our model could be improved by modifying or adding a few key variables. A better measure of intensity of feeding effort might be dive bottom time (time spent below 85% of the dive’s maximum depth; Mattlin et al. 1998, Page et al. 2005). In gray seals (Halichoerus grypus), a study using satellite transmitters, TDRs, and stomach-temperature recorders found that accumulated bottom time was the best predictor of feeding frequency (Austin et al. 2006), and therefore the dive intensity based on dive depth and duration may not necessarily reﬂect the amount of food ingested. Although fur seals are feeding on krill in the pelagic zone, and thus a distinct ‘‘bottom of the dive’’ may not be obvious, calculating the time spent below 85% of the dive’s maximum depth may still prove useful if the maximum dive depth is close to the depth where prey are available. At Bird Island, differences in foraging behavior (i.e., dive behavior) in Antarctic fur seals were most likely caused by differences in the foraging location (Staniland et al. 2004), implying that a variable incorporating location of foraging could add useful information to the model. Fur seals were probably feeding on other sources of prey in addition to krill. Antarctic fur seals tracked from Seal Island did not target the highest krill densities, but rather focused on intermediate densities in offshore areas early in the season (December) and on the slope northwest of Seal Island later in the season (January; Ichii et al. 2007). Seals may have chosen areas where krill were larger or where alternate prey, such as myctophid ﬁsh, were more abundant (Ichii et al. 2007). Scat samples collected at Seal Island in January and February of 1987, 1988, and 1991 indicate that ﬁsh was a common prey for Antarctic fur seals (J. L. Bengtson, unpublished data). Samples from Seal Island in 1991 indicate that the percentage of fur seal scats containing myctophids increases through the season, from 0% of scats containing myctophids in December (n ¼ 13), to 61% in January (n ¼ 33), 81% in February (n ¼ 26), and 100% in March (n ¼ 6; Ichii et al. 2007). Myctophid ﬁsh such as Gymnoscopelus spp. and Electrona antarctica are documented as important prey for birds and seals around the South Shetland Islands (Daneri 1996, Casaux et al. 1998, Jansen et al. 1998,
Ecological Applications Vol. 22, No. 2
LISA M. HIRUKI-RARING ET AL.
Osman et al. 2004). Fish were the most frequent prey (74.5%) and contributed the most mass (54.4%) in the diet of nonbreeding male fur seals (Casaux et al. 1998) and were the second most frequent prey (after krill) for breeding fur seals (Osman et al. 2004). Lactating fur seals shift from preying on krill to ﬁsh at Livingston Island (Banks 2001). Fish is also an important part of female fur seal diet during the breeding season at Heard, Macquarie, Kerguelen, and Marion Islands (e.g., Green et al. 1989, 1990, Cherel et al. 1997, Goldsworthy et al. 1997, Klages and Bester 1998, Lea et al. 2002), and at Bird Island (Staniland et al. 2004). We could improve the model if a time series of myctophid abundance were added, as a measure of a food source independent of krill density. SUMMARY
With the BHM approach used in this study, we used a hierarchical structure to link together widely variable ecosystem elements, account for error and uncertainty, and predict values for years with missing data. We linked an environmental variable (sea ice index), prey indices (krill density and gravid krill abundance), maternal fur seal foraging behavior (diving intensity, time spent foraging), and physical measures of fur seals (maternal fur seal mass, pup mass, pup growth) into one model. Our hypotheses for individual relationships were mostly veriﬁed, with a few interesting exceptions. We found that maternal foraging effort and mass did not have the association that we anticipated with more rapid pup growth after 22 days of age. We were able to predict pup body mass at 22 days of age in years when data were not gathered, allowing us to estimate trends over a much longer time series (more than double the length of our original data set). If we were able to extend the data sets even farther, we could explore whether the rise and fall of pup mass in our model (where it appeared that pup mass at 22 days of age was higher in the mid-1980s, and declined in the early 1990s before starting to recover) is part of a larger pattern in pup mass over time. We also identiﬁed additional variables that could improve the model and address additional factors in the ecosystem. Extension of the BHM to larger systems is straightforward. Future models of the Antarctic ecosystem can include multiple predator species, such as penguins and other seals, multiple prey bases including various ﬁsh species, more oceanographic features such as sea surface temperature, and climate variables such as annual average temperature. The BHM’s ability to predict results for years in which data are missing is useful in situations where remote ﬁeld sites may not allow for consistent data collection (Biuw et al. 2009). BHMs provide a natural framework for relating multiple sources of ecological data from the polar regions, which are difﬁcult and costly to study comprehensively, and are anticipated to be the areas that will be impacted ﬁrst and most dramatically by climate warming.
ACKNOWLEDGMENTS We are grateful to the many people who assisted in Seal Island research and logistics, especially those who were most involved with capturing and instrumenting adult female fur seals and weighing pups: M. Cameron, D. Croll, M. Goebel, R. Holt, H. Huber, J. Jansen, R. Merrick, W. Meyer, S. Osmek, M. Schwartz, and B. Walker. Thanks to the many people who processed Seal Island diet samples: M. Goebel, T. Ha¨rko¨nen, J. Jansen, R. Merrick, and J. E. Sease. Comments by J. Jansen, M. Simpkins, J. Sterling, N. Friday, and three anonymous reviewers improved the paper. This research was supported by the National Oceanic and Atmospheric Administration’s National Marine Fisheries Service, as part of the Antarctic Marine Living Resources Program. Reference to trade names does not imply endorsement by the National Marine Fisheries Service, NOAA. LITERATURE CITED Austin, D., W. D. Bowen, J. I. McMillan, and S. J. Iverson. 2006. Linking movement, diving, and habitat to foraging success in a large marine predator. Ecology 87:3095–3108. Banks, A. R. 2001. Energy content of Antarctic fur seal prey: implications for foraging behavior. Thesis. University of California, Santa Cruz, California, USA. Biuw, M., B. A. Krafft, G. J. G. Hofmeyr, C. Lydersen, and K. M. Kovacs. 2009. Time budgets and at-sea behaviour of lactating female Antarctic fur seals Arctocephalus gazella at Bouvetøya. Marine Ecology Progress Series 385:271–284. Bonner, W. N. 1968. The fur seal of South Georgia. British Antarctic Survey Scientiﬁc Reports 56:1–81. Boveng, P. L. 1996. Space and time relationships between vertebrate predators and a commercially harvested prey species in the Antarctic marine ecosystem. Proceedings of the Biometrics Section of the American Statistical Association 1:83–89. Boveng, P. L., J. L. Bengtson, and L. M. Hiruki. 1998a. Seasonal and annual patterns in time spent foraging by Antarctic fur seals. New Zealand Natural Sciences, Supplement 23:17. Boveng, P. L., L. M. Hiruki, M. K. Schwartz, and J. L. Bengtson. 1998b. Population growth of Antarctic fur seals: Limitation by a top predator, the leopard seal? Ecology 79:2863–2877. Brierley, A. S., D. A. Demer, J. L. Watkins, and R. P. Hewitt. 1999a. Concordance of interannual ﬂuctuations in acoustically estimated densities of Antarctic krill around South Georgia and Elephant Island: biological evidence of sameyear teleconnections across the Scotia Sea. Marine Biology 134:675–681. Brierley, A. S., and J. L. Watkins. 1996. Acoustic targets at South Georgia and the South Orkney Islands during a season of krill scarcity. Marine Ecology Progress Series 138:51–61. Brierley, A. S., J. L. Watkins, C. Goss, M. T. Wilkinson, and I. Everson. 1999b. Acoustic estimates of krill density at South Georgia, 1981 to 1998. CCAMLR Science 6:47–57. Casaux, R., A. Baroni, and A. Carlini. 1998. The diet of the Antarctic fur seal Arctocephalus gazella at Harmony Point, Nelson Island, South Shetland Islands. Polar Biology 20:424–428. Cherel, Y., C. Guinet, and Y. Tremblay. 1997. Fish prey of Antarctic fur seals Arctocephalus gazella at Ile de Croy, Kerguelen. Polar Biology 17:87–90. Christensen, N. L., A. M. Bartuska, J. H. Brown, S. Carpenter, C. D’Antonio, R. Francis, J. F. Franklin, J. A. MacMahon, R. F. Noss, D. J. Parsons, C. H. Peterson, and R. G. Woodsmansee. 1996. The report of the Ecological Society of America Committee on the scientiﬁc basis for ecosystem management. Ecological Applications 6:665–691. Clark, J. A. 2007. Models for ecological data: an introduction. Princeton University Press, Princeton, New Jersey, USA.
FUR SEAL BAYESIAN HIERARCHICAL MODEL
Clarke, A. 1984. Lipid content and composition of Antarctic krill Euphausia superba Dana. Journal of Crustacean Biology 4:285–294. Clarke, A., E. J. Murphy, M. P. Meredith, J. C. King, L. S. Peck, D. K. A. Barnes, and R. C. Smith. 2007. Climate change and the marine ecosystem of the western Antarctic Peninsula. Philosophical Transactions of the Royal Society B 362:149–166. Costa, D. P., J. P. Croxall, and C. Duck. 1989. Foraging energetics of Antarctic fur seals, Arctocephalus gazella, in relation to changes in prey availability. Ecology 70:596–606. Cressie, N., C. A. Calder, J. S. Clark, J. M. Ver Hoef, and C. K. Wikle. 2009. Accounting for uncertainty in ecological analysis: the strengths and limitations of hierarchical statistical modeling. Ecological Applications 19:553–570. Croxall, J. P., and M. N. Pilcher. 1984. Characteristics of krill Euphausia superba eaten by Antarctic fur seals Arctocephalus gazella at South Georgia. British Antarctic Survey Bulletin 63:117–125. Daneri, G. A. 1996. Fish diet of the Antarctic fur seal, Arctocephalus gazella, in summer, at Stranger Point, King George Island, South Shetland Islands. Polar Biology 16:353–355. Doidge, D. W., and J. P. Croxall. 1989. Factors affecting weaning weight in Antarctic fur seals Arctocephalus gazella at South Georgia. Polar Biology 9:155–160. Dorfman, R. 1938. A note on the delta-method for ﬁnding variance formulae. Biometric Bulletin 1:129–137. Ducklow, H. W., K. Baker, D. G. Martinson, L. B. Quetin, R. M. Ross, R. C. Smith, S. E. Stammerjohn, M. Vernet, and W. Fraser. 2007. Marine pelagic ecosystems: the West Antarctic Peninsula. Philosophical Transactions of the Royal Society B 362:67–94. Fraser, W. R., and E. E. Hofmann. 2003. A predator’s perspective on causal links between climate change, physical forcing and ecosystem response. Marine Ecology Progress Series 265:1–15. Gentry, R. L., and J. R. Holt. 1982. Equipment and techniques for handling northern fur seals. U.S. Department of Commerce, NOAA Technical Report NMFS SSRF-758. Georges, J.-Y., and C. Guinet. 2000. Maternal care in the subantarctic fur seal on Amsterdam Island. Ecology 81:295– 308. Godlewska, M. 1996. Vertical migrations of krill (Euphausia superba Dana). Polskie Archiwum Hydrobiologii 43:9–63. Goldsworthy, S. D., M. A. Hindell, and H. M. Crowley. 1997. Diet and diving behavior of sympatric fur seals Arctocephalus gazella and A. tropicalis at Macquarie Island. Pages 151–163 in M. Hindell and C. Kemper, editors. Marine mammal research in the southern hemisphere. Volume I: Status, ecology and medicine. Surrey Beatty and Sons, Chipping Norton, Australia. Goodman, L. A. 1960. On the exact variance of products. Journal of the American Statistical Association 55:708–713. Green, K., H. R. Burton, and R. Williams. 1989. The diet of Antarctic fur seals Arctocephalus gazella (Peters) during the breeding season at Heard Island. Antarctic Science 1:317– 324. Green, K., R. Williams, K. A. Handasyde, H. R. Burton, and P. D. Shaughnessy. 1990. Interspeciﬁc and intraspeciﬁc differences in the diets of fur seals, Arctocephalus species (Pinnipedia: Otariidae), at Macquarie Island. Australian Mammalogy 13:193–200. Guinet, C., S. Goldsworthy, and S. Robinson. 1999. Sex differences in mass loss rate and growth efﬁciency in Antarctic fur seal (Arctocephalus gazella) pups at Macquarie Island. Behavioural Ecology and Sociobiology 46:157–163. Hewitt, R. 1997. Areal and seasonal extent of sea-ice cover off the northwestern side of the Antarctic Peninsula: 1979 to 1996. CCAMLR Science 4:65–73. Hewitt, R. P., D. A. Demer, and J. H. Emery. 2003. An 8-year cycle in krill biomass density inferred from acoustic surveys
conducted in the vicinity of the South Shetland Islands during the austral summers of 1991–1992 through 2001– 2002. Aquatic Living Resources 16:205–213. Hill, H. J., P. N. Trathan, J. P. Croxall, and J. L. Watkins. 1996. A comparison of Antarctic krill caught in nets and taken by macaroni penguins Eudyptes chrysolophus: evidence for selection? Marine Ecology Progress Series 140:1–11. Hiruki, L. M., M. K. Schwartz, and P. L. Boveng. 1999. Hunting and social behaviour of leopard seals (Hydrurga leptonyx) at Seal Island, South Shetland Islands, Antarctica. Journal of Zoology, London 249:97–109. Ichii, T., J. L. Bengtson, P. L. Boveng, Y. Takao, J. K. Jansen, L. M. Hiruki-Raring, M. F. Cameron, H. Okamura, T. Hayashi, and M. Naganobu. 2007. Provisioning strategies of Antarctic fur seals and chinstrap penguins produce different responses to distribution of common prey and habitat. Marine Ecology Progress Series 344:277–297. Ichii, T., K. Katayama, N. Obitsu, H. Ishii, and M. Naganobu. 1998. Occurrence of Antarctic krill (Euphausia superba) concentrations in the vicinity of the South Shetland Islands: relationship to environmental parameters. Deep Sea Research Part I: Oceanographic Research Papers 45:1235–1262. Jansen, J. K., P. L. Boveng, and J. L. Bengtson. 1998. Foraging modes of chinstrap penguins: contrasts between day and night. Marine Ecology Progress Series 165:161–172. Jansen, J. K., R. W. Russell, and W. R. Meyer. 2002. Seasonal shifts in the provisioning behavior of chinstrap penguins, Pygoscelis antarctica. Oecologia 131:306–318. Klages, N. T. W., and M. N. Bester. 1998. Fish prey of fur seals Arctocephalus spp. at subantarctic Marion Island. Marine Biology 131:559–566. Lea, M. A., Y. Cherel, C. Guinet, and P. D. Nichols. 2002. Antarctic fur seals foraging in the Polar Frontal Zone: interannual shifts in diet as shown from fecal and fatty acid analyses. Marine Ecology Progress Series 245:281–297. Lea, M. A., C. Guinet, Y. Cherel, G. Duhamel, L. Dubroca, P. Pruvost, and M. Hindell. 2006. Impacts of climatic anomalies on provisioning strategies of a Southern Ocean predator. Marine Ecology Progress Series 310:77–94. Link, W. A., E. Cam, J. D. Nichols, and E. G. Cooch. 2002. Of BUGS and birds: Markov chain Monte Carlo for hierarchical modeling in wildlife research. Journal of Wildlife Management 66:277–291. Loeb, V. J., V. Siegel, O. Holm-Hansen, R. Hewitt, W. Fraser, W. Trivelpiece, and S. Trivelpiece. 1997. Effects of sea-ice extent and krill or salp dominance on the Antarctic food web. Nature 387:897–900. Lunn, N. J., and I. L. Boyd. 1993. Effects of maternal age and condition on parturition and the perinatal period of Antarctic fur seals. Journal of Zoology 229:55–67. Lunn, N. J., I. L. Boyd, and J. P. Croxall. 1993. Factors affecting the growth rate and mass at weaning of Antarctic fur seal pups at Bird Island, South Georgia. Journal of Mammalogy 74:908–919. Lunn, N. J., I. L. Boyd, and J. P. Croxall. 1994. Reproductive performance of female Antarctic fur seals: the inﬂuence of age, breeding experience, environmental variation and individual quality. Journal of Animal Ecology 63:827–840. Marr, J. W. S. 1962. The natural history and geography of the Antarctic krill (Euphausia superba Dana). Discovery Reports 32:33–464. Mattlin, R. H., N. J. Gales, and D. P. Costa. 1998. Seasonal dive behaviour of lactating New Zealand fur seals (Arctocephalus forsteri ). Canadian Journal of Zoology 76:350–360. McCafferty, D. J., I. L. Boyd, T. R. Walker, and R. I. Taylor. 1998. Foraging responses of Antarctic fur seals to changes in the marine environment. Marine Ecology Progress Series 166:285–299. Morris, C. N. 1983. Parametric empirical Bayes inference: theory and applications. Journal of the American Statistical Association 78:47–55 and 55–65.
LISA M. HIRUKI-RARING ET AL.
Murphy, E. J., P. N. Trathan, J. L. Watkins, K. Reid, M. P. Meredith, J. Forcada, S. E. Thorpe, N. M. Johnston, and P. Rothery. 2007b. Climatically driven ﬂuctuations in Southern Ocean ecosystems. Proceedings of the Royal Society B 274:3057–3067. Murphy, E. J., et al. 2007a. Spatial and temporal operation of the Scotia Sea ecosystem: a review of large-scale links in a krill centred food web. Philosophical Transactions of the Royal Society B 362:113–148. Nicol, S., J. Kitchener, R. King, G. W. Hosie, and W. K. de la Mare. 2000. Population structure and condition of Antarctic krill (Euphausia superba) off East Antarctica (80–1508 E) during the Austral summer of 1995/1996. Deep-Sea Research II 47:2489–2517. Osman, L. P., R. Hucke-Gaete, C. A. Moreno, and D. Torres. 2004. Feeding ecology of Antarctic fur seals at Cape Shirreff, South Shetlands, Antarctica. Polar Biology 27:92–98. Page, B., J. McKenzie, and S. D. Goldsworthy. 2005. Intersexual differences in New Zealand fur seal diving behavior. Marine Ecology Progress Series 304:249–264. Quetin, L. B., R. M. Ross, and A. Clarke. 1994. Krill energetics: seasonal and environmental aspects of the physiology of Euphausia superba. Pages 165–184 in S. Z. ElSayed, editor. Southern Ocean ecology: the BIOMASS perspective. Cambridge University Press, Cambridge, UK. Reid, K. 2002. Growth rates of Antarctic fur seals as indices of environmental conditions. Marine Mammal Science 18:469– 482. Reid, K., and J. P. Y. Arnould. 1996. The diet of Antarctic fur seals Arctocephalus gazella during the breeding season at South Georgia. Polar Biology 16:105–114. Reid, K., and J. P. Croxall. 2001. Environmental response of upper trophic-level predator reveals a system change in an Antarctic marine ecosystem. Proceedings of the Royal Society B 268:377–384. Reid, K., J. P. Croxall, D. R. Briggs, and E. J. Murphy. 2005. Antarctic ecosystem monitoring: quantifying the response of ecosystem indicators to variability in Antarctic krill. ICES Journal of Marine Science 62:366–373. Reid, K., E. J. Murphy, V. Loeb, and R. P. Hewitt. 2002. Krill population dynamics in the Scotia Sea: variability in growth and mortality within a single population. Journal of Marine Systems 36:1–10. Reid, K., P. N. Trathan, J. P. Croxall, and H. J. Hill. 1996. Krill caught by predators and nets: differences between species and techniques. Marine Ecology Progress Series 140:13–20. Reid, K., J. L. Watkins, J. P. Croxall, and E. J. Murphy. 1999. Krill population dynamics at South Georgia 1991–1997, based on data from predators and nets. Marine Ecology Progress Series 177:103–114. Ross, R. M., and L. B. Quetin. 1983. Spawning frequency and fecundity of the Antarctic krill Euphausia superba. Marine Biology 77:201–205.
Ecological Applications Vol. 22, No. 2
Ross, R. M., and L. B. Quetin. 1991. Ecological physiology of larval euphausiids, Euphausia superba (Euphausiacea). Memoirs of the Queensland Museum 31:321–333. Royle, J. A., and R. M. Dorazio. 2008. Hierarchical modeling and inference in ecology: The analysis of data from populations, metapopulations, and communities. Academic Press, San Diego, California, USA. Siegel, V., B. Bergstro¨m, U. Mu¨hlenhardt-Siegel, and M. Thomasson. 2002. Demography of krill in the Elephant Island area during summer 2001 and its signiﬁcance for stock recruitment. Antarctic Science 14:162–170. Siegel, V., and V. Loeb. 1995. Recruitment of Antarctic krill Euphausia superba and possible causes for its variability. Marine Ecology Progress Series 123:45–56. Siegel, V., V. Loeb, and J. Groger. 1998. Krill (Euphausia superba) density, proportional and absolute recruitment and biomass in the Elephant Island region (Antarctic Peninsula) during the period 1977 to 1997. Polar Biology 19:393–398. Stammerjohn, S. E., D. G. Martinson, R. C. Smith, and R. A. Iannuzzi. 2008. Sea ice in the western Antarctic Peninsula region: Spatio-temporal variability from ecological and climate change perspectives. Deep-Sea Research II 55:2041– 2058. Staniland, I. J., K. Reid, and I. L. Boyd. 2004. Comparing individual and spatial inﬂuences on foraging behaviour in Antarctic fur seals Arctocephalus gazella. Marine Ecology Progress Series 275:263–274. Trathan, P. N., J. Forcada, and E. J. Murphy. 2007. Environmental forcing and Southern Ocean marine predator populations: effects of climate change and variability. Philosophical Transactions of the Royal Society B 362:2351–2365. Trathan, P. N., J. Priddle, J. L. Watkins, D. G. M. Miller, and A. W. A. Murray. 1993. Spatial variability of Antarctic krill in relation to mesoscale hydrography. Marine Ecology Progress Series 98:61–71. Trites, A. W. 1991. Does tagging and handling affect the growth of northern fur seal pups (Callorhinus ursinus)? Canadian Journal of Fisheries and Aquatic Sciences 48:2436–2442. Trites, A. W. 1993. Biased estimates of fur seal pup mass: origins and implications. Journal of Zoology 229:515–525. Ver Hoef, J. M. 1996. Parametric empirical Bayes methods for ecological applications. Ecological Applications 6:1047–1055. Walker, B. G., and P. L. Boveng. 1995. Effects of time–depth recorders on maternal foraging and attendance behavior of Antarctic fur seals (Arctocephalus gazella). Canadian Journal of Zoology 73:1538–1544. Watkins, J. L., F. Buchholz, J. Priddle, D. J. Morris, and C. Ricketts. 1992. Variation in reproductive status of Antarctic krill swarms; evidence for a size-related sorting mechanism? Marine Ecology Progress Series 82:163–174. Watkins, J. L., D. J. Morris, C. Ricketts, and J. Priddle. 1986. Differences between swarms of Antarctic krill and some implications for sampling populations. Marine Biology 93:137–146.