The Center for Whale Research (www.whaleresearch.com) collected the
long-term dataset for Southern residents. With all alive members of this
population being observed each year of the study, this provides a
near-complete life history for the population including information on
survival for males and females and reproduction for females. Fisheries
and Oceans Canada collected the long-term data on both Bigg’s and
Northern residents, which includes the survival for most individuals and
reproduction for females. However, a multi-year lack of encounters with
sime matrilineal groups generated uncertainty for year of death of some
individuals (Table 2). We used the photo-identification data as
capture-mark-recapture data in the format of a sighting matrix, i.e.
absence/presence for each individual during each year.
In all three datasets, year of birth was only included for animals that
were born after the start of the given study, otherwise year of birth
was zero, indicating birth year as unknown (Table 2). In Southern
residents, newborn individuals were observed within their first year of
life. In Northern residents and Bigg’s, if an individual was not
observed within the first year of life, year of birth was estimated
based on body length of the individual when it was first
observed41. For all populations, year of death of an
individual was determined either directly from strandings or based on an
individual missing from its matrilineal group on a single (individual
less than 3 years old) or on several (individual more than 3 years old)
occasions41. If year of birth and/or death were
unknown for a given individual, these values were assigned a zero and
would be estimated by the model (see Table 2 for number of individuals
with known birth and death year).
Estimating age-specific survival
As killer whales are a long-lived species, a common feature of the
observation data from all three populations is that some individuals
were born before the studies started (i.e. left truncated) and some were
still alive at the end of the study (i.e. right-censored). This gives
rise to uncertainties regarding some birth and death years in the
datasets. Further, some individuals have gaps in their sighting history
of more than 1 year likely due to their preference for waters beyond the
core study areas. This introduces uncertainty in the recapture
probability and year of death of these individuals. To account for the
uncertainties and missing observations in the data we used a Bayesian
hierarchical framework to estimate the age-specific survival and
mortality for all three populations42. This framework
estimates the unknowns and uncertainties as latent variables (i.e.
variables to be estimated) and combines this with flexible parametric
mortality functions to predict the age-specific survival of the three
killer whale populations49. Although the data on
Southern residents is near-complete, we have used the same approach on
all datasets to ensure the results are directly comparable. Given that
males and females likely have different mortality
trajectories30 the sex of individuals was included as
a covariate in the analyses.
We fitted ten different mortality and survival models to the each of the
three datasets using the BaSTA package49 in R version
3.6.250. The different functions tested either
describes a constant mortality that is independent of age (Exponential),
an exponential increase in mortality with age
(Gompertz)51, an increase or decrease in mortality as
a power function of age (Weibull)52 or an initial
exponential increase in mortality that plateaus after a given age
(Logistic)53. Adding a shape term to the mortality
function allows the mortality to have an initial decline from birth
(bathtub shape)54 or an added constant mortality rate
that is independent of age (Makeham)55. We ran the
models for four Markov chain Monte Carlo chains (MCMC) to evaluate
convergence of the model. The final model specifications, where
convergence for all model types had been reached, were: four MCMC
chains, 1,000,000 iterations, 200,001 burn-in, and 10,001
thinning56. We visually inspected convergence in the
parameter traces for the four chains, to ensure that all chains had
mixed properly. Further, convergence was assessed formally for each
model on the basis of the scale reduction (\(\hat{R}\)), where\(\hat{R}\) < 1.1 indicates
convergence49,57.
Testing the fit of the model
We used the deviance information criterion (DIC) to evaluate the model
fit and predictive power of the different models, a measure analogous to
Aikaiki’s information criterion (AIC)58 (but see
Spiegelhalter et al (2014)59). The importance of
including the categorical covariate of sex was investigated using the
Kullback-Leibler discrepancy, which is informed by the overlap of the
posterior distribution of parameter estimates60–62.
For data including only individuals of known sex, the Gompertz model
with a bathtub shape described the mortality and survival trajectory
well for all three populations, as the second best fit for Southern and
Northern residents and the best fit for Bigg’s (Table 3).
Mathematically, the Gompertz bathtub model consists of three elements
\begin{equation}
{\ \mu\left(x\right)=\ e}^{a_{0}-a_{1}x}+c+exp\left(b_{0}+b_{1}x\right)\nonumber \\
\end{equation}where a0 and a1 defines
the initial decline in mortality from birth, c defines the
mortality through the adult stage and b0 andb1 defines the exponential increase in mortality
at the senescent stage. This model therefore offers great flexibility in
the age of onset of aging as well as changes in the rate of aging
throughout life42,49,54 and we use this model for
quantifying the post-reproductive lifespan for all three populations to
best ensure direct comparison between the three populations. Given the
lack of data on older whales and the related uncertainty around ages at
death, we included weakly informative prior distributions that allowed
the models to explore the parameter space while reflecting a plausible
range of resulting values for survival. These prior distributions were
informed using known life history traits for the
populations29. For the final model (Gompertz with
bathtub shape) we ran the model with the following prior means (and
standard deviations): a0 = -3 (σ = 1),a1 = 0.2 (σ = 0.05), c = 0.001 (σ =
0.001), b0 = -4 (σ = 1), b1 = 0.05 (σ = 0.01). All prior distributions were normally distributed
and truncated at 0, except for a0 and
b0 . Further, recapture probabilities for all three
populations were assumed to be time-dependent as there are variations in
observation efforts across the study periods. This time-dependency of
recapture probabilities was therefore included in the model, as it
allows for variation of the parameter estimates for each time period
(10-year period).
Permutations of individuals with unknown sex
The sex of a substantial portion of individuals were not determined in
the Northern resident and Bigg’s populations (Table 2). These
individuals were all under 15 years old. In mammals, mortality is often
higher in early life, and by excluding juvenile individuals from the
analyses, we are likely underestimating such early-life mortality.
Instead of including them as a third category of “unknown sex”, which
would – in essence – calculate the mortality trajectory of an
artificial short-lived sex, we used a permutation approach to assign a
sex to these individals randomly. This way we are able to include the
early life mortality into the full age-specific mortality trajectory. We
implemented the same Bayesian hierarchical model that was the best
overall fit for individuals of known sex, the Gompertz bathtub model,
which was confirmed by a trial run on a random dataset as the best fit
for the full data. Although only 21 individuals were of unknown sex in
the Southern residents, we also ran the same number of permutations on
this population for comparability. We ran 1000 permutations, where sex
was randomly assigned to individuals of unknown sex for each
permutation. Both populations have a female-biased adult sex ratio
(Table 2), and we assumed a sex ratio of 1:1 at birth. We used the
permuted output to calculate the variation in post-reproductive
representation in all three populations.