Bellwethers of change: population modelling of North Pacific humpback whales from 2002 through 2021 reveals shift from recovery to climate response

2.1. Sample collection and matching methods

For this study, all available North Pacific humpback whale photo-ID data from a research collaboration of 46 organizations and 4292 community science contributors were aggregated within a study period of 2001 through 2022, and reconciled in a single dataset through the research collaboration and community science Web platform Happywhale.com. This dataset is fully described by Cheeseman et al. [18]. Photo-ID data for the purpose of this study consisted of encounters of individually identified whales, with associated dates and geographical locations. Multiple encounters of the same individual within the same season (summer or winter of a given year) were condensed to one capture occasion (capture), using the observation date closest to the midpoint of the respective season.

2.2. Geographical stratification

Each capture was assigned to a geographical stratum based on location, to segregate samples between summer (feeding) areas and winter (breeding) areas and as a method to allocate samples for a bias correction estimation described below. Geographical strata include seven summer areas and six wintering areas (tables 1 and 2; figure 1). Stratification generally followed that used in the SPLASH study [20], with strata based on similarities in migratory destinations. For example, the California and Oregon stratum was separated from the Washington and southern British Columbia stratum because whales from the former migrate almost exclusively to the coasts of Mexico and Central America, whereas a portion of the latter also migrate to the Mexican offshore islands and Hawaiian Island regions. Our large sample size allowed us to discern some migration patterns that were not apparent in the SPLASH study, which resulted in the following changes in the stratification scheme: humpback whales in southern coastal Mexico had similar distribution and migratory destinations to whales in Central America [45,46], so we separated southern Mexico from the rest of Mainland Mexico and pooled it with Central America. The ‘Russia' stratum in the SPLASH study included whales from the Kamchatka Peninsula, the Commander Islands and the Gulf of Anadyr. Based on winter migratory destinations, we revised the stratification of Russia as follows: we included the Kamchatka Peninsula as a separate stratum (Kamchatka whales migrated almost entirely to islands in the western Pacific), we combined the Commander Islands with the other western Aleutian Islands to be a new western Bering Sea stratum, we combined northern Russia coastal whales with a stratum that includes the northern and eastern Bering Sea, and samples from the Mariana Islands (which was not sampled in the SPLASH study) were pooled with the western Pacific stratum. Whale distributions in the North Pacific are continuous, with blurred boundaries across most regions (e.g. [47]), with regional delineation used here only to illustrate varied sampling over time and to develop bias-correction factors that address non-random sampling, as described below. We do not attempt to estimate abundance for each stratum and our primary abundance estimation method depends on all summer (feeding) areas pooled, and all winter (presumed breeding) areas pooled. Consequently, these geographical stratification changes have minimal impact on actual abundance estimates, serving only to align sample distributions between SPLASH study years and non-SPLASH study years. See electronic supplementary material, I, for stratum boundaries.

• Download figure
• Open in new tab
• Download PowerPoint

2.3. Absolute abundance estimation methods: full North Pacific Ocean

Prior to estimating abundance, whale identification data from Happywhale were pre-processed with a custom R script [48]. For each capture, sample season was assigned as either summer or winter. As with the SPLASH study, midpoints for each capture period were defined as 1 March (for winter) and 1 August (for summer), hereafter ‘season-year'. In the rare cases of an individual encountered in multiple geographical strata in a given season-year, the individual was assigned to only one geographical stratum, based on the closest capture occasion to the seasonal midpoint.

To be comparable to the 2004–2006 SPLASH abundance estimates [30], we used the Chapman bias-corrected version of Petersen estimator [49] with wintering area samples (pooled over three consecutive years) as one capture occasion and with summer area samples (pooled over the two summers between the three years of winter samples) as the second capture occasion. Modelling has shown that this estimator is robust to sample bias created by heterogeneity in capture probabilities if the factors affecting differential availability are different in the two capture occasions [30]. We relied on the different geographical distribution and behaviour of whales between summer and winter to reduce sample bias between mark and recapture samples. Recaptures included those individuals photographed in both the pooled winter and pooled summer samples, i.e. during both ‘capture occasions'. We repeated this estimation process for 20 overlapping 3-year periods, 2002 through 2021 with data from 2001 to 2022, to produce a time series of 20 annual abundance estimates. Abundance estimates are labelled based on the year of the first summer in each 3-year period.

Although the simple abundance estimation approach outlined above should work well if the geographical areas were as representatively sampled throughout the time series as during the 2004–2006 SPLASH study, the achieved geographical distribution of samples was uneven in many years (table 1). In particular, the Revillagigedo Archipelago (MxIsl), Baja California, Mexico (MxBaja), Kamchatka, Russia (Kamch) and the eastern and western Bering Sea (BerE and BerW) were intermittently and generally under-sampled. This was expected to bias the abundance estimates outside of the SPLASH years. We therefore used a novel geographical bias correction approach, constructed as follows.

We assumed that the bias-corrected SPLASH abundance estimate of 21 063 whales [30] was an unbiased estimate for study years (2004–2006). We denoted this reference estimate as N_SPLASH. Let the photo-ID data underpinning N_SPLASH be denoted Y_SPLASH. Let N_t denote the time series (2002–2021) of abundance values estimated in our analysis and let Y_t be the associated photo-ID data. Note, each Y_t value represented 3 years of data centred on year t (e.g. Y₂₀₁₀ are the data from 2009, 2010, 2011). For each t, we subsampled Y_SPLASH (call this subsample Y_{SPLASH, sub(t)}) so that it had the same geographical distribution as Y_t (see next paragraph). Y_{SPLASH, sub(t)} was then used to obtain an alternative SPLASH abundance estimate (N_{SPLASH, alt(t)}). We then calculated the ratios between N_SPLASH and each N_{SPLASH, alt(t)}. These ratios (F_t = N_SPLASH/N_{SPLASH, alt(t)}) were interpreted as the amount of bias associated with the geographically biased subsamples and were applied as correction factors to the uncorrected annual Chapman estimates ($N_t^{\mathrm{\prime }}$), to obtain corrected estimates, i.e. $N_t=N_t^{\mathrm{\prime }}\times F_t$.

To generate subsamples from the SPLASH data (Y_{SPLASH, sub(t)}) with the same geographical distributions as the Y_t, we used a largest common sample size (LCSS) approach. We calculated the number of unique whale identifications in all summer and wintering geographical strata for both the SPLASH period (2004–2006) and for the estimation period (t − 1 to t + 1). For each summer and wintering stratum we then calculated the largest sample size that could be taken from both Y_SPLASH and Y_t. We then subsampled both datasets by randomly sampling this LCSS from each (without replacement) to obtain Y_{SPLASH, sub(t)} and Y_t(sub). The Chapman–Petersen estimator was used with these data subsets to obtain N_{SPLASH, alt(t)} and $N_t^{\mathrm{\prime }}$, which were used to find F_t and N_t (see previous paragraph).

The primary assumptions of this bias-correction approach are that (i) the SPLASH abundance estimate itself is not biased by non-representative geographical sampling, (ii) migration patterns have not changed since the SPLASH period, and (iii) the geographical distribution of whales has not changed during the course of the time series, i.e. that they have remained proportionally distributed throughout the study area as they were during the SPLASH period. The assumption of geographical sampling bias in the SPLASH study was addressed as part of that study [20,30]. The second assumption, of migration pattern consistency, was evaluated by comparing migratory destinations of southeast Alaska and northern British Columbia from two samples from different time periods, 2004–2006 and 2019–2021, three seasons each. The third assumption, of consistency of geographical distribution of whales, was evaluated through an ANOVA test of Hawai‘i versus Mainland Mexico annual mean relative abundance estimates over the study period to detect shifts in differential population growth rates between the two most sampled wintering areas.

2.4. Relative abundance estimation methods: Hawai‘i and Mexico

To add resolution to understanding humpback whale abundance over the study period, we estimated a time series of abundances for two wintering areas which were relatively well sampled in most years (Hawai‘i and Mainland Mexico, table 2) with consistent research effort over the study period, and which comprised a high proportion of total humpback whale abundance in the North Pacific, with migratory contributions from multiple feeding areas. Chapman–Petersen abundance estimates were made for both wintering areas based on photo-ID samples in consistently well-sampled feeding areas. Feeding areas were selected based on (i) being a substantial migratory contribution to the given wintering area [20] and (ii) having a consistently high sample size of identified whales. The Hawai‘i estimate was based on using that wintering area as one sampling occasion and the southeast Alaska/northern British Columbia feeding area stratum as the second sampling occasion. The Mainland Mexico estimate was based on using that wintering area as one sampling occasion and the California/Oregon feeding area stratum as the second sampling occasion. Humpback whales in both of these paired wintering and feeding areas exhibited elements of migratory herd units [50] with wintering areas shared with whales from multiple feeding areas, and feeding areas populated with whales almost entirely from their respective paired wintering areas. As with the ocean-basin-wide model, estimates were based on running 3-year periods with two feeding area samples and three wintering area samples. To interpret these abundance estimates as absolute abundance would have required that the wintering area sample was a random sample of the migratory whales from all feeding areas that contributed to it. Because that assumption was not evaluated in the scope of this study, we treated our estimates as an index of relative abundance.

2.5. Standard error estimation

To estimate standard errors (s.e.) in abundance estimates, we used a jackknife approach [51]. For both the bias-corrected estimates for the entire North Pacific and for the estimates of relative abundance for the Hawaiian Island and Mainland Mexico strata, jackknife samples were taken by randomly excluding 5% of the fluke identifications (without replacement, after season-year duplicates were excluded for each geographical stratum). The estimation of abundance (including the bias correction) was repeated for each of the 20 jackknife samples, and s.e. were estimated using the standard formula:

where n was the number of jackknife estimates, $\widehat{x_i}$ was the ith jackknife estimate, and ${\overline{x}}_{\mathrm{jack}}$ was the mean of the jackknife samples. The 2004–2006 SPLASH abundance estimate was used as a reference for correcting the bias for all other sample periods. To include the uncertainty in the SPLASH bias correction estimate (CV = 0.04 [30]), the CV of the bias-corrected abundance was estimated as the square root of the sum of the squared CV of the uncorrected abundance estimate plus the squared CV of the SPLASH estimate.

3.1. Sampling distribution

The full 2001–2022 dataset consisted of 30 484 individual humpback whales in 132 684 unique season-year encounters, distilled from 192 869 unique encounters. The number of unique individuals that were photographically sampled each season-year varied from 1940 to 5668, and from 1252 to 3236 for the 21 winters and summers, respectively (table 3). Photo-ID samples were post-stratified into seven feeding areas (table 1) and six wintering areas (table 2). Sample sizes varied among seasons and regions due to varying effort; 20 of 294 season-regions had no effort and zero identifications. Kamchatka and the Mexican Islands regions were consistently under-sampled, with an average of 14 and 49 seasonal identifications, respectively. Hawai‘i (833), California/Oregon (648) and southeast Alaska/northern British Columbia (772) averaged the highest seasonal identifications, resulting from both high abundance and high effort. Across all seasons and regions, we found a mean of 448 unique identifications per season, and a maximum of 1798 unique identifications in one season. Sample sizes increased on average fourfold over the years of the study period, with the highest number of unique identifications from 2017 onward (table 3). Further description of this sample distribution is found in [18].

3.2. Recapture rates and bias correction

Between-season matches ranged from 4.95% to 11.43% of individuals per season (table 3). Over the full study period, 63% of all individuals were captured in multiple seasons and/or regions, suggesting that a high proportion of the North Pacific humpback whale population was included in the sample [18]. While sample sizes were very limited for some regions, no region stood out with a disproportionately low recapture rate. Recapture rates ranged from 56.8% of individuals captured in multiple seasons and/or regions in the Mariana Islands to 87.2% in southern Mexico and Central America. Following the photo-ID methodology of [44], we assumed that 97–99% of potential matches were found.

Bias correction factors (F_t) averaged 1.11 (0.94–1.22, s.d. 0.07, table 3), adjusting abundance estimations upward, in most cases, from uncorrected abundance estimates biased by non-representative sampling efforts.

3.3. Abundance estimation and trends

Bias-corrected estimates of humpback whales in the full North Pacific increased from 16 875 (s.e. = 5955) in sample year 2002 to a peak of 33 488 (s.e. = 4455) in 2012 before an inflection with a decline to 26 662 individuals (s.e. = 4192) by sample year 2021 (figure 2 and table 3). Growth from 2002 to 2013 was consistently positive apart from one anomalous year, 2010, averaging 5.9% yr⁻¹. From 2014 to 2021, estimated abundance consistently declined, averaging −3.0% yr⁻¹. Hawai‘i region relative abundance estimates showed a peak in 2013, then declined to abundance levels equal to the mid-2000s from 2017 onward (table 4 and figure 3). The Hawai‘i region growth phase was steeper than the North Pacific-wide aggregated dataset, averaging 6.9% yr⁻¹ from 2002 to 2013, then declined at an average of −4.3% yr⁻¹ from 2014 to 2021. By contrast, relative abundance for the Mainland Mexico region grew at an average of 7.1% yr⁻¹ from 2002 to 2015 before appearing to stabilize with an average of 0.9% yr⁻¹ growth from 2016 to 2021 (table 4 and figure 3).

• Download figure
• Open in new tab
• Download PowerPoint

• Download figure
• Open in new tab
• Download PowerPoint

3.4. Migration patterns and relative abundance

Analysis of migratory destinations of 2022 individuals found in southeast Alaska and northern British Columbia in 2004–2006 and 2025 individuals in the same area in 2019–2021 found 591 and 698 individuals in wintering areas, respectively, with no significant difference in migratory destinations between the two time periods (χ² = 0.263, p = 0.877; table 5). The proportionate relative abundance of whales between Hawai‘i and Mexico in 2002–2006 was 78% versus 22%, respectively (table 4 and figure 4). The Hawai‘i percentage declined by an average of −0.8% annually over the study period, so that the relative abundances were 65% versus 35% respectively by the 2017–2021 sampling period (ANOVA p < 0.001).

• Download figure
• Open in new tab
• Download PowerPoint

4.1. Bias correction validity and assumptions

The bias-correction approach allowed us to account for variable sampling effort over time and space, necessitating three assumptions. First, the assumption that the SPLASH abundance estimate is unbiased by non-representative geographical sampling is supported by an experimental design for the SPLASH study that pursued representative sampling of all wintering and summer areas known at the time [20,30]. The SPLASH study had no effort in the Mariana Islands, which has since been found to be a breeding area for humpback whales [21]; however, after seven seasons of effort there with 37 individuals identified [18], we are confident that this omission from SPLASH sampling effort would not significantly change the abundance estimation at the scale of the full North Pacific Ocean. The northwest Hawaiian Islands are now known to be a wintering area [69] with as yet unknown abundance and relationship to the main Hawaiian Islands. Of 37 individuals identified in two northwest Hawaiian Islands voyages of limited scope and duration (2007 and 2019), 24 individuals were resighted in Hawai‘i, interchange that suggests a substantial degree of mixing. We believe this is the only North Pacific humpback whale wintering area that may host an under-sampled, site-faithful population, though with a majority of identified individuals also encountered in the main Hawaiian Islands. Further study may thus reveal minimal impact on total abundance estimation statistics. Russian feeding areas were represented by only 102 individuals during SPLASH years [20]; in contrast, this study benefited from data from 2296 individuals [18]. This heterogeneity in sampling effort contributed to larger s.e. ranges derived from our jackknife approach in the full North Pacific abundance estimate, compared to narrower ranges in Hawai‘i and Mexico relative abundance estimates.

A further consideration of potential sampling bias is of equal detection probability. High accuracy of the AI image recognition tool used for matching [44] meant any fluking whale could be confidently identified irrespective of distinctiveness. Females with calves, however, are known to fluke less often in wintering areas [70] though Craig and Herman demonstrated that extended focal follows of non-fluking individuals would greatly reduce sampling bias [71]. The same study found lower wintering resight rates of females. Further, not all humpbacks migrate every year [72]. We believe sampling bias is effectively minimized by using multi-year pooled samples for each annual abundance estimation and the methodology of drawing from a dataset of nearly 200 000 identified encounters (192 869 unique encounters distilled to 132 684 unique season-year encounters) that may mimic the high detection rates of extended focal follows [71].

The second assumption, that migratory patterns have not changed, was conclusively found to be valid based on no difference between wintering destinations of over 2000 whales each in time periods near the beginning and end of our study period. This assumption is further supported by demonstrated strong maternal fidelity and natal philopatry in humpback whales with sufficient stability over time to show genetic differentiation [73–75].

Validity of the third assumption would require that population trends were uniform between and across all summer and wintering areas during the study period. We found that relative abundance between Hawai‘i and Mexico shifted from 78% versus 22% to 65% versus 35% over the course of the study (figure 4) apparently because of different population growth rates (figure 3). This 13% shift in relative abundance between Hawai‘i and Mexico, and possible undetected shifts in relative abundance between other study area regions, does violate this assumption and thus may have introduced error in our geographical bias correction factor affecting absolute abundance estimates. We realize that the 13% relative abundance shift is a violation of the assumption of consistent geographical distribution underlying our sampling bias correction methodology; however, we believe its effects on observed abundance trends (as documented by the study) are likely to be minimal because (i) the large samples created by multi-year pooled effort where approximately one-third to half of the total estimated North Pacific humpback population was captured in most 3-year samples and (ii) the bias corrections themselves are relatively small. The overall trends without the bias correction are similar to the bias-corrected trends (table 2), although the values are lower in all years except one. We cannot fully discount this as a potential source of bias, particularly after 2013 when the proportion of the population in Hawai‘i decreased substantially (figure 4). To address these sources of uncertainty, future study for humpback whales in the North Pacific should include absolute abundance estimation by region.

4.2. Relative abundance estimation

We sought to refine the ocean-basin-wide abundance estimate with regional estimates based on the two best-sampled subset wintering areas, Hawai‘i and the Mainland Mexico. These subset models could not be considered absolute abundance estimations because each 3-year sample was not assumed to be a closed population or controlled for sample bias such as, in the case of Mainland Mexico, an overlap of Central American whales. These whales are known to co-occur using the region as a migratory corridor [46,76], and therefore the subset abundance estimate sampled a population with undefined boundaries. In Hawai‘i, sample effort was concentrated in the Maui-Nui channel and to lesser extents off the islands of Hawai‘i, Oahu and Kauai, without a clear definition of whether this is a comprehensive sample of an equally mixed Hawaiian island population. The value of the relative abundance estimations was in describing different population growth trajectories with far lower s.e. estimations than the North Pacific-wide model. An apparent sharp decline in humpback whale relative abundance off west Maui (which traditionally has hosted the greatest concentrations of humpback whales including calves of the year, e.g. [77–79]) was reported after the 2014–2016 PMH, using boat-based transect surveys [38] and passive acoustic monitoring [65], as well as off the Kohala coast of Hawai‘i Island using shore-based systematic scans [64]. The relative abundance trend measured agrees with data from density estimates by fixed-wing aerial surveys of the entire Maui Nui region for 1993–2003 and 2019–2020 [80], though these surveys did not occur during the period of apparent steep increase then decrease, over the 2005–2016 period. Our results support the conclusions of Frankel et al. [64] that environmental variables likely caused Hawaiian humpback populations to decline due to the PMH. Fluctuating but resilient humpback whale crude birth rates in Hawai‘i [38,80] suggest the measured −4.3% yr⁻¹ decline from 2014 to 2021 may be temporary, stabilizing if ocean conditions in feeding areas return to the long-term mean. These trends point to an apparent ecological shift from a half century of recovery response following near-extirpation by industrial whaling to carrying capacity-limited abundance. It is notable that the end of the shift did not appear to occur from population growth but rather from a decrease in carrying capacity triggered by rapid climate change.