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.