Current global population size, post-whaling trend and historical trajectory of sperm whales | Scientific Reports

Selection of surveys and extraction of data

We selected published surveys that produced estimates of sperm whale population size or density (see Supplementary Information for methodology; surveys listed in Table 1). We extracted: the type of survey (ship, aerial; acoustic, visual), the years of data collection; the coordinates of the boundary of the study area; the estimates of g(0) and CV (g(0)) used to correct for availability bias, if given; and an estimate of sperm whale population or density in study area with CV. From these we calculated for each survey the survey area with waters greater than 1000 m deep (typical shallow depth limit of sperm whales3). When no value of g(0) was used (8 ship visual surveys) we corrected the population/density estimate using an assumed generic value of g(0) and recalculated the CV to include uncertainty in g(0) (as in Eq. 1 of8). Three ship visual surveys did calculate a single g(0) estimate: 0.62 (CV 0.35)32; 0.57 (CV 0.28)35; 0.61 (CV 0.25)37. These are consistent and suggest a generic g(0) = 0.60 (CV 0.29), also agreeing with g(0) = 0.60 estimated from pooled surveys in the California Current10.

Global habitat of sperm whales

To extrapolate sperm whale densities from surveyed study areas to the sperm whales’ global habitat, we created a one-degree latitude by one-degree longitude grid. We removed the following grid points as not being prime sperm whale habitat1,3,40: points on land or with central depths less than 1000 m; largely ice-covered points in the Beaufort Sea, and the waters north of Svalbard and Russia; the Black Sea and Red Sea both of which have shallow entrances that appear not to be traversable by sperm whales.

Generally, food abundance is a good predictor of species distribution. However, this is not possible for sperm whales as we have no good measures of the abundance or distribution of most of their prey, deep-water squid57. Instead, oceanographic measures have been used to describe sperm whale distributions over various spatial scales with a moderate level of success13,14. We follow this approach. Measures that might predict sperm whale density were collected for each grid point, some at just the surface, others at the surface, 500 m depth, 1000 m depth or an average of the measures at the different depths (Supplementary Table S2). Water depth was the strongest predictor in Mediterranean encounters, when compared to slope and distance to shore13. Temperature and salinity have been used as predictors for the distribution of fish and larger marine animals, which could translate into prey availability and thus density for sperm whales58,59. Primary productivity and dissolved oxygen generally dictate the biomass of wildlife in an area, while nitrate and phosphate levels limit the amount of primary productivity in an area60. Eddy kinetic energy is a measure of the dynamism of physical oceanography which is becoming a commonly used predictor of cetacean habitat61. We did not use: latitude and longitude as these primarily describe the general geographic distribution of the study areas, and geographic aggregates of sperm whale catches62 as these proved to have no predictive power. The mean values of the 14 predictor measures were calculated over calendar months for each grid point, and then over the grid points in each study area.

To obtain predictors of the sperm whale density at each grid point, we then made quadratic regressions of the density of sperm whales in each study area (i), d(i), on the mean values of the predictor measures, weighting each study area by its surface area. Because the surveys were conducted over different time periods, the densities were corrected based on the estimated trajectory of global sperm whale populations by multiplying d(i) by the ratio of the global population in 1993 over that in the mid-year of the survey (as in Fig. 4). Predictor variables were selected using forward stepwise selection based upon reduction in AIC.

Sperm whale population size

The population of sperm whales globally, N, was then calculated as follows:

$$N=\sum_{k}d\left(k\right)\cdot a\left(k\right),$$

(1)

where a{k} are the parameters of the regression; the summation is over k, the grid points; d(k) is the estimated sperm whale density at grid point k from the habitat suitability model; and a(k) is the area of the 1° cell centred on grid point k. Population estimates for other ocean areas (North Atlantic, North Pacific, Southern Hemisphere) were calculated similarly.

The CVs of these population estimates were calculated following the methodology in8, (although there is an error in Eq. (3) of8 such that the squareroot symbol covers both the numerator and denominator rather than just the numerator). The error due to uncertain density estimates for the different surveys is:

$$CV\left({D}_{T}\right)=\frac{\sqrt{\sum_{i}{\left(CV({n}_{i})\cdot {n}_{i}\right)}^{2}}}{\sum_{i}{n}_{i}}.$$

(2)

This is combined with the uncertainty in the extrapolation process (output from the linear models), CV(extrap.), to give an overall CV for the population estimate:

$$CV\left(N\right)=\sqrt{{CV({D}_{T})}^{2}+{CV(\mathrm{extrap}.)}^{2}.}$$

(3)

Post-whaling trend in population size

We compiled a database of series of surveys producing population estimates of the same study area during the period 1978 (by which time most commercial sperm whaling had ceased) and 2022. Each series had to span at least 10 years, and all of the surveys in the series had to be comparable in terms of area covered throughout the time span. There also had to have been at least 3 surveys for a data set to be included.

The data consisted of the survey area, A, the estimated population in area A in year y (for multi-year surveys, y would be the midpoint of the data collection years), nE(A,y), and the provided CV of that estimate, CV(nE(A,y)). The data series used for these analyses are summarized in Table 3.

For each survey area, A, we calculated the trend in logarithmic population size, r(A), over time using weighted linear regression:

$${\text{Log}}\left( {n_{E} \left( {A,y} \right)} \right) \, \sim {\text{ constant}}\left( A \right) \, + r\left( A \right) \cdot y. \left[ {{\text{weight }} = { 1}/\left( {{1} + {\text{ CV}}\left( {n_{E} \left( {A,y} \right)} \right)} \right)^{{2}} } \right]$$

(4)

Table 3 also includes other published estimates of sperm whale population trends, from sighting rates or mark-recapture analyses of photoidentification data, with these estimates also having to span at least 10 years of data collection, and include data collected in three or more different years.

Population trajectory

To examine possible trajectories of the global sperm whale population following the start of commercial whaling in 1712, we used a variant of the theta-logistic, a population model that has been employed in other recent analyses of the population trajectories of large cetaceans45,63. The theta-logistic model is:

$$n\left(y+1\right)=n\left(y\right)+r\cdot n\left(y\right)\left(1-{\left(\frac{n\left(y\right)}{n\left(1711\right)}\right)}^{\theta }\right)-f\left(y\right)\cdot c\left(y\right).$$

(5)

Here, n(y) is the population of sperm whales in year y, r is the maximum potential rate of increase of a sperm whale population, and θ describes how the rate of increase varies with population size relative to its basal level before whaling in 1711, n(1711). The recorded catch in year y is c(y) and f(y) is a correction for bias in recorded catches.

Whaling reduced the proportion of large breeding males64, likely disrupted the social cohesion of the females3, and may have had other lingering effects which reduced pregnancy or survival, and thus the rate of increase. Poaching has been found to reduce the reproductive output of African elephants, Loxodonta Africana, which have a similar social system to the sperm whales3, and this effect lingers well beyond the effective cessation of poaching46. There is some evidence for these effects of what we call “social disruption” on sperm whale population dynamics20,46,65. We added a term to the theta-logistic to account for such effects:

$$n\left(y+1\right)=n\left(y\right)\left[1+r\cdot \left(1-{\left(\frac{n\left(y\right)}{n\left(1711\right)}\right)}^{\theta }\right)-q\cdot \frac{\sum_{t=y-T}^{y}f(t)\cdot c(t)}{n\left(y-T\right)}\right]-f(y)\cdot c(y).$$

(6)

Here, \(\frac{\sum_{t=y-T}^{y}f(t)\cdot c(t)}{n\left(y-T\right)}\) is the proportion of the population killed over the last T years, and q is the reduction in the rate of increase when almost all the whales have been killed. This reduction is modelled to fall linearly as the proportion killed declines to zero.

The global sperm whale population has some geographic structure18. Females appear to rarely move between ocean basins, and males seem to largely stay within one basin. Furthermore, sperm whaling was progressive, moving from ocean area to ocean area as numbers were depleted4. We model this by assuming K largely separate sperm whale subpopulations of equal size. Exploitation in 1712 starts in subpopulation 1 and moves to subpopulations 1 and 2 when the population 1 falls to α% of its initial value, and so on for the other ocean areas. The catch in each year in each area being exploited is pro-rated by the sizes of the different subpopulations being exploited. The population model for subpopulation k, which is one of the KE subpopulations being exploited in year y, is:

$$n\left(k,y+1\right)=n\left(k,y\right)\left[1+r\cdot \left(1-{\left(\frac{n\left(k,y\right)}{n\left(k,1711\right)}\right)}^{\theta }\right)-q\cdot \frac{\sum_{t=y-T}^{y}C(k,t)}{n\left(k,y-T\right)}\right]-C\left(k,y\right),$$

(7)

where the estimated catch in year y in subpopulation k is given by: \(C\left(k,y\right)=f(y)\cdot c(y)\cdot n(k,y)/\sum_{{k}^{\mathrm{^{\prime}}}=<{K}_{E}}n\left({k}^{\mathrm{^{\prime}}},y\right).\)

Global catch data for sperm whales have been compiled by year for the periods 1762–179916, 1800–189917 and 1900–19996 (Fig. 4). This catch record is not complete. There are no compiled data on sperm whaling between its inception in 17124 and 1761. We filled this in by extrapolating a linear trend between the catch in 1711 (0 whales) and that in 1762 (351 whales). Catch records for the 18th and 19th Century are based on extrapolations from oil production records, and are not corrected for whales struck and lost that died, which will produce a downward bias in estimates of catches of the order of a few percent16,17. The 20th Century catch records are also likely downward biased by an unknown, but likely quite small, amount, primarily due to falsification of records6.

For any choice of population and exploitation parameters (n(1711), r, θ, {f}, K, α, q, T) a population trajectory can be calculated, giving a population size in 1993, n(1993). We have an estimate of the population in 1993 from the first part of this paper, but not in 1711. Thus, given values of the other parameters, the catch history {c}, and a value of n(1993), we found the value of n(1711) (using the “fzero” function of MATLAB2021b) so that the trajectory reached n(1993) in year 1993.

The best estimates, assumed distributions, and range of acceptable values for each of the parameters of the population model are given and justified in Table 4.

Using these parameter estimates, we ran the model in three ways: with the best estimate of each parameter; in a sensitivity analysis, using a range of values for each parameter in turn (the range given in Table 4) with the best estimates of the other parameters; and using random but reasonable values for all parameters (see Table 4 for assumed distribution of each) in each of 1,000 runs. Runs producing invalid (infinite or unreal) population estimates in 1711 were discarded.

The primary outputs of the population model were the estimated trajectory from 1711 to 2022, the current population (n(2022)), the estimated population in 1711 before the commencement of commercial whaling, and n(2022)/n(1940). Assuming three generations equals about 82 years66, n(2022)/n(1940) gives the proportional change in population size over 3 generations (the IUCN Red List measure for population trends). The proportional decline from pre-whaling numbers (i.e. 1711) to 2022 was also calculated.