The SRKW population is closed to immigration and emigration, every individual in the population is known, and the population has been censused annually for decades11. Individuals were identified by their unique fin shapes, saddle patches, and the presence of any nicks or scratches, and sexed using distinctive pigmentation patterns around the genital slits. Male and female offspring remain within the natal, matrilineal unit, although mating occurs within and between these pods. The term “resident” refers to their residency in inshore waters of southern British Columbia (Canada) and Washington state (USA) in the summer months, when they feed almost exclusively on Chinook salmon13,14,54,55. Given that there is no dispersal from the population56, mortality was recorded if an individual’s matriline was observed in the population within a year but the individual did not appear.

We used values of demographic parameters calculated from the census data to build the population model in the Vortex PVA program17,18. We included temporal variation in demographic rates (“environmental variation”), based on inter-annual variability in parameters observed since 1976, and we included individual variation in age of maturity and probability of reproductive success. The Vortex simulation model of possible future population trajectories includes demographic stochasticity (binomial variation in individual fates); random assignment of sex and a bi-sexual mating system, resulting in fluctuations in sex ratio and mate availability that can affect small populations; and projections of loss of genetic diversity, allowing for inclusion of inbreeding depression. We quantified population growth as the mean exponential rate of increase (r = ln[N t+1 /N t ]).

Modelling was conducted in stages. First, a “Baseline” model was developed to represent the population trajectories if demographic rates remain the same as have been observed in recent decades. We confirmed this Baseline model by comparing simulated dynamics with recent population trends. Secondly, we conducted sensitivity tests on uncertain demographic rates in the model to determine which parameters had large effects on the projected population growth. Thirdly, we used a set of models of Individual Threats that tested ranges of values for the primary threats identified in the recovery plans to determine which would have the greatest effects on population projections. Fourthly, we examined Cumulative Threats scenarios to project the fate of the population if further industrial development increases existing threats and adds new ones. A set of Demographic Management scenarios was then examined to determine the population growth that could be achieved by improvements in demographic rates. Finally, we explored Threat Management scenarios to assess the plausibility of reaching sustained annual population growth of 2.3% given various options for increasing salmon abundance, reducing ocean noise levels, or reducing contaminant levels. The following section describes key parameter estimates used in the model. More detailed description of the modelling methods is presented in Supplementary Information. The input files for the Vortex project are available at http://www.vortex10.org/SRKW.zip and from the Dryad Digital Repository at https://doi.org/10.5061/dryad.46vq7.

Baseline PVA

We started the simulations with the ages, sexes, and pod membership of the killer whales living in 2015. We specified the mother of each animal, where known (for 76 of 80 living animals)57. Based on previous genetic data on paternity58, we specified in the simulation that females would not mate with their father, a son, or a maternal half-sibling. What effect lower levels of inbreeding or the inevitable accumulated inbreeding in a closed population will have on any cetacean is unknown. We modelled inbreeding depression as being caused by recessive lethal alleles, with 6.29 “lethal equivalents” (the negative of the slope of log(recruitment) against the inbreeding coefficient), the mean combined effect of inbreeding on fecundity and first-year survival in a survey of impacts on wild species59.

Demographic rates were calculated from individual animal histories compiled by the Center for Whale Research57, using data collected from 1976 through 2014. The time series begins when the population was depleted by live-captures for display in aquariums60. The time series therefore includes periods of moderate population growth (1976 to 1993), subsequent decline, and approximate stability. Demographic rates were estimated for the age-class groupings used in recent models24,61, except that we set an upper limit for female breeding at 45 y rather than 50 y, because no females in the population have been documented to produce calves at older ages. Thus, we calculated survival and (for adult females) fecundity rates for calves (first year), juveniles (defined as from 1 y through 9 y of age), young mature females (10–30 y), older reproductive females (31–45 y), post-reproductive females (46 y and older), young mature males (10–21 y), and older males (22 y and older). Killer whales can survive many years after reproductive senescence, but estimating maximum longevity is difficult in such a long-lived species62. We set an upper limit of age to 90 y in our models, although only about 2% of females would be expected to reach this age, and only about 2% of males (with higher mortality) would be expected to exceed 50 y. Females stop breeding long before the maximum age, so the long-term population growth would not be affected by the upper age limit unless post-reproductive females benefit the pod in ways other than through their own reproduction.

Mortality for each age-sex class was averaged across the 39 years of data to obtain mean annual rates. We did not try to partition observed mortality into presumed causes of death. The use of these historic data for our Baseline model makes the implicit assumption that the frequency of deaths due to the various causes remains the same as has been observed across recent decades. The variation in mortality observed across years has two components: 1) environmental variation (fluctuations in the probability of survival), and 2) demographic stochasticity (binomial variation in individual fates). To determine how much of the observed variation was due to environmental variation, the variance due to demographic stochasticity can be calculated from the expectation for a binomial process, and then subtracted from the total variation across years. Calculated annual mortality rates (and environmental variation) ranged from a low of 0.97% (SD = 0) for young adult females to 17.48% (SD = 17.96) for calves. Although the lack of evidence for annual variation in the mortality adult females beyond that expected from random sampling of a constant probability might seem optimistic, for long-lived species a low level of annual variation in rates would have negligible effect on long-term population trajectories. We confirmed through sensitivity tests (Supplementary Information) that the environmental variation entered into the population model has no effect on our results.

The breeding system is polygamous, with some males able to obtain multiple mates, females mating with different males over their lifetimes, and mating between and within pods. Males become sexually mature (actively breeding, which may occur several years after they are physiologically capable of breeding) from 12 to 18 y of age. Thus, in the model, each male was assigned an age of sexual maturity by randomly selecting a value from 12 to 18. Variance in reproductive success among individual females and males will cause genetic diversity to be depleted faster and inbreeding to accumulate faster than would occur if mating was assumed random. Information is available on male mating success51, and we incorporated variation in male and female reproductive success in the model (Supplementary Information). Our models project an effective population size that is 37% of the total size, close to an estimate obtained from genetic data58.

Breeding rates, expressed as the proportion of the females of an age class that produce a calf each year, were calculated from annual census data. Rates ranged from 0% for post-reproductive females (age >45 y), to 7.88% (SD = 4.15) for older adult females (age 31–45 y), to 12.04% (SD = 3.54) for young adult females (age 10–30 y).

The upper limit on population size was set to 300, so that carrying capacity (K) would not restrict future population growth except under the best conditions tested. In the projections of current or expected conditions, the SRKW populations never reached this limiting size, and rarely exceeded 150 animals in any of the independent iterations of each simulation. Population recovery was assessed by the mean growth rate each year calculated before any carrying capacity truncation. Thus, the growth rate reflects the demographic potential and is not affected by the limit on population size in the model.

The SRKW population was projected for 100 years. For the initial exploration of parameter uncertainty, the simulation was repeated in 10,000 independent iterations to obtain high precision in mean and variance estimations. For comparisons among alternative management scenarios, less iteration is needed to obtain the relative influence of input values, and tests were run with 1,000 iterations. Sensitivity tests were conducted by varying each basic demographic rate (life table values for fecundity and mortality) over a range of ± 10% around the baseline value. For several model variables that describe other aspects of the population dynamics and are also very uncertain, a wider range of values was tested (see Supplementary Information).

Individual Threats

We explored the effects of three threats identified in the recovery strategies. For each of prey abundance, noise disturbance, and PCB contaminants, we scaled impacts such that the estimated current level of the threat resulted in the mean demographic rates reported over recent decades. Effects of prey limitation were modelled using published relationships linking inter-annual variability in Chinook salmon to inter-annual variability in calf and adult mortality63 and fecundity13,61. A prey index was calculated by dividing the total salmon abundance in each year by its average abundance over the 1979–2008 period63. The relationship of mortality to prey abundance was modelled with a multiplier of baseline mortality that is a linear function scaled to 1 when salmon abundance was at the mean observed level over period of observation: MortalityFactor = 3.0412 – 2.0412 * PreyIndex. The relationship of birth rate to prey was modelled with logistic functions, with the intercept scaled to yield the observed birth rates for young females (12.04%) and older females (7.88%) when PreyIndex = 1. For relationships of form BirthRate = exp(A + B*PreyIndex)/[1 + exp(A + B*PreyIndex)], the function parameters were A = −3.0 and B = 1.0 for young females, and A = −3.46 and B = 1.0 for older females. (See Supplementary Information for more details on these relationships.) To explore the impacts of prey abundance across a range of plausible values, we varied the prey index from approximately the lowest level (0.60) reported since 1978 to approximately the highest level (1.30).

Effects of noise on demography were modelled using the approach outlined in previous analyses of loss of acoustic communication space4,64. We used summertime observations to estimate the proportion of time boats were present (during daylight hours) while the whales were foraging and the reduction in foraging expected with that amount of acoustic disturbance. We calibrated the model of noise impacts so that the mean Baseline demographic rates are obtained at the reported level of disturbance. We then simulated the relative change in foraging time and consequently demographic rates across the spectrum from no noise impact at all, to the upper limit expected if boat disturbance increased from current, already high, levels to 100% of time. We do not have data on the amount of acoustic disturbance in the winter feeding areas, but the modelling based on observed summertime disturbance provides a means to project a range of population consequences if changes in disturbance overall mirror those that are possible in the summertime habitat. Land-based observations have shown that SRKWs reduce their time spent feeding in the presence of boats by 25%65. Vessels are present 85% of the daytime, and SRKWs are foraging in the presence of vessels an estimated 78% of that time. Thus, for the 85% current (baseline) exposure to vessels, feeding is expected to be reduced by 16.6% ( = 85% × 78% × 25%) due to disturbance by boats. To translate the reduction in feeding into its demographic consequences, we multiplied the prey index by a factor of (1 − 0.195 * Noise)/(1 – 0.166) to obtain the proportional availability of prey. This proportion is thus 1 in the current, baseline conditions (Noise = 0.85), 0.965 when vessels are always present (Noise = 1.00), and 1.20 assuming no disturbance from vessels. The noise-modified index of prey availability was then used to determine the consequent mortality and fecundity rates. We recognize that anthropogenic noise can also have less direct effects on wildlife, including disruption of social behaviours and even impeding responsiveness to other sensory modalities66.

Our model of accumulation, depuration, and impact on calf survival of PCBs was based on the approach described by Hall et al.67,68 with modifications in rates for SRKW69. Calves obtain their initial load of contaminants from their dams through gestation and lactation, and females producing calves thereby depurate an estimated 77% of their contaminants67. Otherwise, males and non-breeding females accumulate PCBs in the blubber of at a rate that we varied from 0 to 5 ppm/y in our tests. Few data are available on PCBs in the SRKW population with which to calibrate the model of PCB bioaccumulation, and the levels of PCBs reported in SRKW might have been dropping slowly in recent years. Reported levels in adult female SRKW range from 55 ± 19 ppm sampled in 1993–1996, 37 ± 42 ppm sampled in 2004–2007, and 30 ± 31 ppm sampled in 2008–201330. Our population model generates a mean 28, 55, and 81 ppm PCBs in adult females when bioaccumulation rate is 1, 2, and 3 ppm/y, respectively. Effects of maternal PCB load on calf mortality were modelled using a logistic response function (survival = exp(2.65 − 0.02 * PCB)/[1 + exp(2.65 − 0.02 * PCB)]), fitted to the two observed data points for SRKW (survival = 0.8252) and the nearby northern resident killer whales (survival = 0.9218)24, with the mean PCB levels (55.4 ppm and 9.3 ppm, respectively)70 reported from the time period in the middle of the span over which mortality rates were calculated. If we use the more recent, lower estimates of PCB loads in SRKW to estimate the impacts, our response function would have a steeper slope. There are not yet sufficient data on effects of PCBs on other demographic rates to allow inclusion of any other effects of PCBs (or other contaminants) in our PVA model.

Cumulative Threats

We modelled two scenarios to represent the cumulative impacts of possible increases in threats, based in part on a recent environmental impact assessment submitted to Canada National Energy Board45 evaluating effects of a proposed oil pipeline and associated tanker traffic. For the purposes of this PVA, projected increases in anthropogenic threats are not meant to mimic any one industrial development, but rather a general process of industrialization reflecting the number of port expansions, pipeline proposals, and liquefied natural gas terminal proposals pending for the BC coast4. For a low level scenario, we used the catastrophe option in Vortex to add the possibility of large (>16,500 m3) and smaller (>8,250 m3) oil spills. The frequencies of a big spill (0.21% chance per year) and a smaller spill (1.08%) were based on an industry projection of the likelihood of such spills caused by proposed increase in tanker traffic71. Based on the percent overlap of oil coverage and critical habitat, we estimate that if a large oil spill were to occur, about 50% of the SRKWs would be killed due to direct exposure to the oil. We estimate that 12.5% of the SRKWs would be killed by exposure to oil from a smaller spill. For a scenario with higher level impacts of development, we doubled the frequency of oil spills.

These energy development scenarios also included an increase in vessel noise and disturbance of feeding, with the current vessel presence of 85% of time increased to 92.5% in the low level scenario and to 100% in the high level scenario. We also included a probability of additional deaths of killer whales due to ship strikes, with one death per decade in the low level and two deaths per decade in the high level scenario. Although some persistent organic pollutants might increase under increased industrial activity in the SRKW habitat, PCBs have been phased out of production and are in decline in at least some fish species in low-development basins72. Lacking data on likely long-term trends in the contaminant loads of SRKW prey, we did not include any change in such pollutants in these scenarios.

Climate change is projected to cause a decline in Chinook abundance28, and we modelled this possibility with a projected 25% (low scenario) or 50% (high scenario) decrease in Chinook over the next 100 years.

Demographic Management and Threat Management scenarios

We used the PVA to simulate how much improvement in demographic parameters or how much reduction in anthropogenic threats, singly or in combination, would be required to reach a stated recovery objective of sustained annual population growth of 2.3% for 28 years11. In calculating the growth for these models, we started the tally 20 years into the simulation to avoid short-term demographic fluctuations as the age structure adjusts to new demographic rates, and growth was tallied over the subsequent 28 years. For the set of Demographic Management scenarios, we assessed the relationship between improved demography and population growth. Birth rate was incremented by 1.1x, 1.2x, 1.3x, 1.4x, and 1.5x, whereas calf mortality and adult mortality were decreased by 0.9x, 0.8x, 0.7x, 0.6x, and 0.5x. Next, in Threat Management scenarios, we modelled the effects of reduced threats, with the consequences resulting from the functional relationships to demography. We increased salmon abundance (up to the highest level of the Chinook index observed between 1979 and 2008, namely 1.3 times the long-term average). We simulated the improved demography if acoustic disturbance were reduced or eliminated. We considered the population consequences of improved calf survival resulting from reduction of PCBs, testing rates of future accumulation in SRKW from the estimated current 2 ppm/y to down to 0 ppm/y. Finally, we tested scenarios that both reduced acoustic disturbance by half and increased salmon abundance up to 1.3x.