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

Ted CheesemanpersName, Jay BarlowpersName, Jo Marie AcebespersName, Katherina AudleypersName, Lars BejderpersName, Caitlin BirdsallpersName, Olga Solis BracamontespersName, Amanda L. BradfordpersName, Josie ByingtonpersName, John CalambokidispersName, Rachel CartwrightpersName, Jen CedarleafpersName, Andrea Jacqueline García ChavezpersName, Jens CurriepersName, Rouenne Camille De CastropersName, Joëlle De WeerdtpersName, Nicole DoepersName, Thomas Doniol-ValcrozepersName, Karina DracottpersName, Olga FilatovapersName, Rachel FinnpersName, Kiirsten R. FlynnpersName, John FordpersName, Astrid Frisch-JordánpersName, Chris GabrielepersName, Beth GoodwinpersName, Craig HayslippersName, Jackie HilderingpersName, Marie C. HillpersName, Jeff K. JacobsenpersName, M. Esther Jiménez-LópezpersName, Meagan JonespersName, Nozomi KobayashipersName, Marc LammerspersName, Edward LymanpersName, Mark MallesonpersName, Evgeny MamaevpersName, Pamela Martínez LoustalotpersName, Annie MastermanpersName, Craig O. MatkinpersName, Christie McMillanpersName, Jeff MoorepersName, John MoranpersName, Janet L. NeilsonpersName, Hayley NewellpersName, Haruna OkabepersName, Marilia OliopersName, Christian D. Ortega-OrtizpersName, Adam A. PackpersName, Daniel M. PalaciospersName, Heidi PearsonpersName, Ester Quintana-RizzopersName, Raul Ramírez BarragánpersName, Nicola RansomepersName, Hiram Rosales-NanducapersName, Fred SharpepersName, Tasli ShawpersName, Ken SoutherlandpersName, Stephanie StackpersName, Iain StanilandpersName, Janice StraleypersName, Andrew SzabopersName, Suzie TeerlinkpersName, Olga TitovapersName, Jorge Urban-RamirezpersName, Martin van AswegenpersName, Marcel ViniciuspersName, Olga von ZiegesarpersName, Briana WitteveenpersName, Janie WraypersName, Kymberly YanopersName, Igor YeginpersName, Denny ZwiefelhoferpersName, Phil ClaphampersName (2024date). “Bellwethers of change: population modelling of North Pacific humpback whales from 2002 through 2021 reveals shift from recovery to climate response.” Royal Society Open Sciencetitle·j.

Object of study — annual ocean-basin-wide abundance of North Pacificterm·key humpback whalesterm·key from 2002date through 2021date

Research question — whether North Pacificterm·key humpback whalesterm·key are still recovering from whaling or now respond to climate-driven conditions

Contribution — a bias-corrected 20-year mark-recaptureterm·mth abundance time series and a shift to continuous monitoring

Full text — tagged

Bellwethers of change: population modelling of North Pacificterm·key humpback whalesterm·key from 2002date through 2021date reveals shift from recovery to climate response Abstract For the 40 years after the end of commercial whaling in 1976, humpback whaleterm·key populations in the North Pacificterm·key Ocean exhibited a prolonged period of recovery. Using mark–recaptureterm·mth methods on the largest individual photo-identificationterm·key dataset ever assembled for a cetacean, we estimated annual ocean-basin-wide abundance for the species from 2002date through 2021date. Trends in annual estimates describe strong post-whaling era population recovery from 16 875 (± 5955) in 2002date to a peak abundance estimate of 33 488 (± 4455) in 2012. An apparent 20% decline from 2012 to 2021date, 33 488 (± 4455) to 26 662 (± 4192), suggests the population abruptly reached carrying capacityterm·dom due to loss of prey resources. This was particularly evident for humpback whalesterm·key wintering in Hawai‘i, where, by 2021date, estimated abundance had declined by 34% from a peak in 2013, down to abundance levels previously seen in 2006, and contrasted to an absence of decline in Mainland Mexico breeding humpbacks. The strongest marine heatwaveterm·key recorded globally to date during the 2014–2016 period appeared to have altered the course of species recovery, with enduring effects. Extending this time series will allow humpback whalesterm·key to serve as an indicator species for the ecosystem in the face of a changing climate.

Keywords: carrying capacityterm·dom, marine heatwaveterm·key, mark–recaptureterm·mth modelling, abundance estimationterm·key, climate changeterm·dom, environmental variables Introduction Population monitoring has become increasingly important in conservation biology as anthropogenic activities and climate changeterm·dom continue to impact marine ecosystems [ 1 – 5 ]. Abundance trends provide critical insights into the dynamics of animal populations, enabling a better understanding of the ecological interactions and underlying drivers that influence their distribution, abundance and life history [ 6 , 7 ]. Abundance estimationterm·key is the basis for assessing the current protection status for humpback whalesterm·key in the USA [ 8 ] and Canada [ 9 ] and has helped the understanding of anthropogenic effects such as ship strikes on blue whales [ 10 ]. For marine mammal populations exploited to near extinction, abundance estimationterm·key has served as an indispensable tool for decision-makers, supporting the development of effective conservation and management strategies aimed at protecting these animals and the ecosystems they inhabit [ 11 , 12 ]. Understanding population trends is essential for addressing the growing challenges faced by marine mammals in a rapidly changing world [ 13 ]; however, detecting population changes can be difficult due to the imprecision of most marine mammal abundance estimates and the infrequent intervals at which most populations are surveyed [ 14 ]. For large cetaceans, such as baleen whales, these issues are largely due to the costs and logistical challenges associated with surveying vast, often remote marine environments, as well as elusive and often cryptic animal behaviour [ 15 ]. Limited available data may result in poor model parameter estimation and reduced predictive power, ultimately hindering effective conservation and management efforts [ 16 ]. Furthermore, the infrequent generation of updated population estimates may fail to capture rapid changes in whale populations due to anthropogenic disturbances or environmental fluctuations [ 17 ]. To address these challenges, we have drawn on an innovative ocean-basin-wide collaboration [ 18 ] to cost-effectively generate a time series of relatively precise abundance estimates for a long-distance migrating species, the humpback whaleterm·key ( Megaptera novaeangliae ) in the North Pacificterm·key Ocean.

In the North Pacificterm·key, humpback whalesterm·key are known to breed in waters off Japan, The Philippines, the Mariana Islands, Hawai‘i, Mexico and Central America [ 19 – 22 ]. These whales migrate to feed in coastal to continental shelf waters off Russia, the Bering Sea, Alaska, British Columbia, Washington, Oregon and California [ 20 ]. As with all large cetaceans worldwide, humpback whalesterm·key in the North Pacificterm·key were the target of extensive commercial whaling until late in the twentieth century, with an estimated 31 785 whales taken from 1900 to 1976 [ 23 , 24 ]. Humpback whaleterm·key populations in the North Pacificterm·key were severely depleted [ 24 ], with the remnant abundance crudely estimated at a low of approximately 1200–1600 individuals [ 25 , 26 ] around the end of humpback whaleterm·key commercial catches in 1976 [ 24 ]. The end of commercial whaling created a possibility of full recovery to pre-whaling abundance; it is, however, difficult to define what full recovery means without precise estimates of pre-whaling abundance. Rice [ 27 ] suggested that the humpback population was ‘stable at about 15 000 [animals]' prior to 1905, and ‘was reduced to about 1000' by 1965. However, he gave no supporting evidence for the former figure, and the latter was estimated using catch statistics that were available at the time (either published or provided by the USSR); it is now known that the Soviet statistics were significantly under-reported [ 28 ]. As Rice stated, no one at the time had attempted to estimate the pristine population size prior to modern whaling, and the number of animals killed in various aboriginal operations prior to the twentieth century is largely unknown. Consequently, there is currently no reliable estimate of pre-whaling abundance, although calculating one is a primary aim of an ongoing Comprehensive Assessment by the Scientific Committee of the International Whaling Commission [ 29 ].

To understand population status, two previous abundance estimationterm·key studies for North Pacificterm·key humpback whalesterm·key were undertaken using mark–recaptureterm·mth (also known as capture–recaptureterm·mth) techniques with synoptic datasets [ 30 , 31 ]. These studies depended on visual matching of photographs by skilled technicians to determine marks and recaptures of individual humpback whalesterm·key based on the unique shape and pigmentation of the ventral surface of the tail (fluke) [ 32 , 33 ]. The first photographic mark–recaptureterm·mth abundance estimate for humpback whalesterm·key in the full North Pacificterm·key Ocean was based on a historical database of 3650 fluke images gathered from 1990 to 1993 [ 34 ]. Important western Pacific Ocean feeding areas of Russia, the Aleutian Islands and the Bering Sea were not surveyed during the study, and because it was retrospective, sampling did not follow a systematic design. Using a geographically stratified Darroch method, that study estimated an abundance of approximately 6000 (s.e. = 474) humpback whalesterm·key in 1992 from three major wintering areas (Mexico, Hawai‘i and Japan), but the authors' appraisal of likely biases suggested an adjusted estimate of approximately 8000 [ 34 ]. Assuming the same coefficient of variation (CV) as the uncorrected estimate, a minimum estimate of the standard error for the bias-corrected estimate is 632. Petersen estimates comparing wintering and summer areas suggest that this may have been an underestimate, but a more accurate Petersen model approach could not be applied uniformly due to the missing data from the unsampled summer areas.

Some of the geographical biases in the 1990–1993 sampling were corrected in the 2004–2006 project Structure of Populations, Levels of Abundance, and Status of Humpback whalesterm·key (SPLASH) [ 20 ], when coordinated efforts by over 40 different research groups across the North Pacificterm·key obtained fluke identification photographs from all known wintering and summer areas in the North Pacificterm·key during three winters and two summers. The experimental design sought to apportion sampling effort in proportion to whale density to representatively sample whales in all areas. The SPLASH project's North Pacificterm·key abundance estimate was based on a Chapman–Petersen approach that compared all wintering fluke photographs with all summer fluke photographs [ 30 ]. In that study, simulations were used to correct biases associated with births and deaths (a common violation of the closed population assumption), not sampling calves and missed matches. Many of these biases were found to be offsetting, and the net estimate of resulting bias was only +3.5%. The bias-corrected SPLASH abundance estimate was 21 063 (CV = 0.04) humpback whalesterm·key in 2005.

Based on evidence of ongoing recovery indicated by the SPLASH study and other studies, the species was partially delisted from the United States Endangered Species Act (ESA) in 2016 [ 8 ], and was reassessed from Threatened to Special Concern by the Committee on the Status of Endangered Wildlife in Canada in 2011 [ 35 ]. However, multiple current anthropogenic stressors including ship strikes [ 36 ], entanglements [ 37 ] and climate changeterm·dom [ 5 , 38 – 40 ] may limit or reverse full recovery. Our 2001–2022 study period included seasons before, during and after periods of a severe marine heatwaveterm·key [ 41 ], with demonstrated effects on regional humpback whaleterm·key populations [ 39 ], and signs of compounding climate changeterm·dom and habitat use conflict such as fisheries interactions leading to exacerbation of these existing threats [ 42 , 43 ].

Our current study builds on past photo-identificationterm·key (photo-ID) efforts mentioned above, incorporating the historical data gathered in those studies, combined with long-term region-specific photo-ID efforts, and making broader use of opportunistic sources, especially whale watch voyage naturalists. The sheer scale of data now available is due to advances in automated image recognition fluke photo-ID matching [ 44 ], data management and research collaboration [ 18 ]. These advances have shifted the abundance estimationterm·key paradigm from data scarcity and periodic study to continuous and accessible tracking of the ocean-basin-wide population through time. Given that no single research group's data archive or effort could scale and maintain what would be necessary for such an ongoing study, we combined research-based and community science data in a broad collaboration facilitated by the data management Web platform Happywhale.com [ 18 ]. Where data deficiency previously imposed limits on understanding, this study used a near-real-time population modelling infrastructure and with this, we aimed to evaluate population change with a robust time series of 20 annual abundance estimates from 2002date to 2021date. We used the same Chapman–Petersen mark–recaptureterm·mth approach for abundance estimationterm·key as was used for the 2004–2006 SPLASH study [ 30 ], based on a comparison of three wintering samples with two summer samples for each calculated annual estimation. To address geographically biased sampling in pre- and post-SPLASH years, we used a sub-sampling scheme of the SPLASH years to obtain a bias correctionterm·key. We examined whether the results of this study might be biased by long-term changes in migratory patterns and in the relatively large proportions of the entire population that winter in Mainland Mexico and Hawai‘i, respectively. Finally, we discussed findings of a shift in population trajectories from recovery from industrial whaling to responding to current ecosystem conditions.

Methods Sample collection and matching methods For this study, all available North Pacificterm·key humpback whaleterm·key 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 CheesemanpersName 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.

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 correctionterm·key 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 whalesterm·key 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 Pacificterm·key 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 estimationterm·key 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.

Sample size of unique photo-identificationsterm·key of humpback whaleterm·key identifications in each of seven summer regions for each season-year. Regions, west to east, include the Kamchatka Peninsula (Kamchatka), the western Bering Sea (WBerSea), the eastern and northern Bering Sea and southern Chukchi Sea (E&NBerSea), the Gulf of Alaska (GulfOfAK), southeast Alaska and northern British Columbia (SEAK&NorthBC), southern British Columbia and Washington (SouthBC&WA), and California and Oregon (CA&OR).

Sample size of unique photo-identificationsterm·key of humpback whaleterm·key identifications in each of six winter regions for each season-year. Regions, west to east, include the western Pacific (WPac), the Hawaiian Islands (Hawai‘i), Baja California, Mexico (MexBajaCal), the Revillagigedo Archipelago, Mexico (MexIsl), Mainland Mexico (MexMld) and Central America and southern Mexico (CenAm&SMex).

Locations of photo-identificationterm·key samples from humpback whalesterm·key in the North Pacificterm·key 2001–2021date colour coded by geographical strata. Summer feeding area strata north of 32°N, west to east are: Kamchatka Peninsula, Russia (black), west Bering Sea (cyan), north and western Bering Sea (green), Gulf of Alaska (yellow), southeast Alaska and northern British Columbia, Canada (blue and purple), southern British Columbia, Canada and Washington, United States (salmon), and California and Oregon, United States (orange). Wintering area strata south of 32°N, west to east are West Pacific (purple), Hawai‘i, United States (blue), Baja California, Mexico (black), Revillagigedo Islands, Mexico (cyan), Mainland Mexico (green), southern Mexico and Central America (yellow).

Absolute abundance estimationterm·key methods: full North Pacificterm·key 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, 2002date through 2021date 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 estimationterm·key 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 correctionterm·key 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 (2002date–2021date) 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 2010 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 ′ ), to obtain corrected estimates, i.e. N t = N t ′ × 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 ′ , 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–2021date, 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 abundanceterm·mth estimates over the study period to detect shifts in differential population growth rates between the two most sampled wintering areas.

Relative abundanceterm·mth estimation methods: Hawai‘i and Mexico To add resolution to understanding humpback whaleterm·key 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 whaleterm·key abundance in the North Pacificterm·key, 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 whalesterm·key 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 abundanceterm·mth.

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 Pacificterm·key and for the estimates of relative abundanceterm·mth 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 correctionterm·key) 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, x i ^ was the i th jackknife estimate, and x ¯ 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 correctionterm·key 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.

Results Sampling distribution The full 2001–2022 dataset consisted of 30 484 individual humpback whalesterm·key 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 ].

Sample size of unique summer and winter area whales (three and two seasons pooled, respectively), number of between-season matches, and resulting North Pacificterm·key abundance estimates without and with bias correctionsterm·key for non-representative geographical sampling. Standard errors are estimated from jackknife samples ( n = 20). Between season match rates are (matches between season)/(total unique whales in winter and summer). ΔYoY is the year-over-year rate of change in abundance, presented as 3-year moving average values for 2003–2020.

Recapture rates and bias correctionterm·key 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 Pacificterm·key humpback whaleterm·key 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 correctionterm·key factors ( F t ) averaged 1.11 (0.94–1.22, s.d. 0.07, table 3 ), adjusting abundance estimationsterm·key upward, in most cases, from uncorrected abundance estimates biased by non-representative sampling efforts.

Abundance estimationterm·key and trends Bias-corrected estimates of humpback whalesterm·key in the full North Pacificterm·key increased from 16 875 (s.e. = 5955) in sample year 2002date 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 2021date ( figure 2 and table 3 ). Growth from 2002date to 2013 was consistently positive apart from one anomalous year, 2010, averaging 5.9% yr −1 . From 2014 to 2021date, estimated abundance consistently declined, averaging −3.0% yr −1 . Hawai‘i region relative abundanceterm·mth 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 −1 from 2002date to 2013, then declined at an average of −4.3% yr −1 from 2014 to 2021date. By contrast, relative abundanceterm·mth for the Mainland Mexico region grew at an average of 7.1% yr −1 from 2002date to 2015 before appearing to stabilize with an average of 0.9% yr −1 growth from 2016 to 2021date ( table 4 and figure 3 ).

Bias-corrected estimates of humpback whaleterm·key abundance for the full North Pacificterm·key Ocean. The red line represents a three-year moving average of abundance estimates. Vertical lines represent confidence intervals based on ± twice the standard errors (s.e.). S.e. include both the uncertainty in the abundance for the SPLASH years (CV = 0.04, mid-sample year = 2005 [ 30 ]) and the uncertainty in the bias-correction process. Inset: North Pacificterm·key humpback whaleterm·key abundance in the context of post-whaling population estimates [ 25 , 26 ] and 1993 estimate [ 31 ].

Estimates of relative abundanceterm·mth for the Hawai‘i and Mainland Mexico regions and the percentage of their total abundance represented by Hawai‘i. Values and standard errors (s.e.) are the same as in figure 3 . ΔYoY is the year-over-year rate of change in estimated abundance, presented as 3-year moving average values for 2003–2020, with 2002date and 2021date omitted due to lack of multiple years to average creating excessive variabiliaty. Year over year growth rates above 11.8% exceed estimated maximum plausible population growth rates [ 52 ] and thus are expected to be estimation errors.

Relative estimates of humpback whaleterm·key abundance for the Hawai‘i (red) and Mainland Mexico (black) regions based on mark–recaptureterm·mth comparisons with northern British Columbia and southeast Alaska, and California and Oregon, respectively. The lines represent 3-year moving averages of abundance estimates. Vertical lines represent confidence intervals based on ± twice the standard error (s.e.).

Migration patterns and relative abundanceterm·mth 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–2021date found 591 and 698 individuals in wintering areas, respectively, with no significant difference in migratory destinations between the two time periods ( χ 2 = 0.263, p = 0.877; table 5 ). The proportionate relative abundanceterm·mth of whales between Hawai‘i and Mexico in 2002date–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 abundancesterm·mth were 65% versus 35% respectively by the 2017–2021date sampling period (ANOVA p < 0.001).

Migratory destinations of 2004–2006 and 2019–2021date southeast Alaska and northern British Columbia humpback whalesterm·key show no change over the study period.

Proportion of Hawai‘i and Mainland Mexico relative abundanceterm·mth of humpback whalesterm·key during study period of 3-year pooled samples 2002date through 2021date, with a slope of −0.8% per year.

Discussion This study estimated a total abundance of 26 662 (s.e. = 4192) humpback whalesterm·key in the North Pacificterm·key as of 2021date, growing from 2002date to 2021date at an average rate of 3% yr −1 , with a 20% decline from a peak abundance of 33 488 (s.e. = 4455) in 2012. Across the North Pacificterm·key Ocean, we described a distinct inflection from post-whaling era rapid population recovery to a state where population dynamics may now be more constrained by environmental factors, including variability in resources induced by climate changeterm·dom. From 2002date through 2013, growth averaging 5.9% yr −1 agreed with the 5.2–8.6% yr −1 95% confidence interval growth rate described by the SPLASH study [ 30 ] and was well within an estimated maximum biological potential growth rate of 11.8% yr −1 [ 52 ]. The timing of the inflection from long-term growth, most notably evident in our Hawai‘i relative population estimation ( figure 3 ), is consistent with observed changes in oceanographic conditions in feeding areas, specifically the severe northeast Pacific Marine Heatwaveterm·dom (PMHterm·dom) of 2014–2016 [ 53 ], which depressed biological productivity and thus prey availability in northeast Pacific Ocean whale feeding areas. Two other known sources of humpback whaleterm·key mortality, ship strikes and entanglements, may have contributed to the apparent decline; both have exceeded estimated potential biological removal rates [ 54 , 55 ] and entanglement rates surged on the US West Coast in correlation with the PMHterm·dom apparently due to habitat compression [ 43 , 56 , 57 ] but cannot account for an estimated population decline approaching 7000 individuals across the North Pacificterm·key in just 9 years (2012–2021date).

The PMHterm·dom, the most extreme in recorded history [ 41 ], had profound environmental impacts including elevated sea surface temperatures (by as much as 3 s.d. in some areas with maximum sea surface temperature anomalies sometimes exceeding 3–6°C to depths ranging from 50 to 200 m), a decline in sea surface winds, reduced upwellingterm·dom and decreases in nutrient-rich water. This led to reduced phytoplankton biomass and restructured zooplankton communities in favour of lower-calorie species [ 58 , 59 ]. These changes, in concert with increased metabolically driven food demands of ectothermic forage fish resulting in quality and quantity declines in forage fish [ 60 ], and an increased demand for forage fish by large ectothermic predatory fish (e.g. groundfish and salmon), ultimately increased competition beyond the demands already imposed by top predators such as marine mammals and sea birds. The result created what Piatt et al . [ 61 ] proposed as an ‘ectothermic vise' on the diminished forage resource. High latitude waters within the PMHterm·dom witnessed a decline in local abundance of humpback whalesterm·key, where they traditionally feed on krill and small schooling fish such as herring, capelin, sand lance and juvenile salmon [ 39 ]. Reduced prey biomass of zooplankton combined with increased competition for diminished forage fish resulted in severely diminished presence of calves, prevalence of ‘skinny' whales, increased strandings [ 62 ], and the absence of individual whales with high site fidelity to specific locales. The impacts were also witnessed in humpback whaleterm·key breeding grounds including significantly fewer whales and calves of the year (e.g. [ 38 , 63 – 65 ]). Co-occurring with the negative effects on humpback whalesterm·key in their feeding grounds, several mass mortalities of fish-eating seabirds were documented suggesting far reaching impacts of the PMHterm·dom. For example, over 60 000 dead or dying emaciated common murres ( Uria aalge ) were documented in 2015–2016 from California to Alaska, apparently from starvation, with high rates of reproductive failure [ 61 ]. The eastern Bearing Sea witnessed shifts in zooplankton community composition and forage fish distribution with a mass mortality of emaciated tufted puffins ( Fratercula cirrhata ) during the onset of molt, apparently associated with starvation [ 66 ]. Other species, whose food resources partially overlap with those diminished in association with the PMHterm·dom, and which also experienced declines in body condition and unusual mortality rates included: various Alaskan large predatory groundfish, California sea lions ( Zaolphus californianus ) and Guadalupe fur seals ( Arctocephalus townsendi ) (summarized in [ 61 ]). Given all of these co-occurring impacts with the PMHterm·dom, and because humpback whalesterm·key feed on both krill and small schooling fish, and may be flexible in their ability to shift their dominant food resource between the two when faced with variations in oceanographic and ecological conditions affecting their food resources [ 67 ], they may serve as a valuable bellwether indicator species of ocean-basin ecosystem health. Their lack of recovery by 2021date is notable.

While conservation management has focused on populations defined by breeding areas [ 8 ], feeding area resource limitations are likely the greatest determinants of carrying capacityterm·dom. Regional evidence in Alaska shows variable recovery in feeding areas where local abundance was severely reduced by the PMHterm·dom, still depressed in Prince William Sound [ 40 ] compared to partial recovery in Glacier Bay National Park [ 39 ]. The documented downward inflection in the population growth curve of North Pacificterm·key humpback whalesterm·key may indicate that population impacts of ecological stressors like the PMHterm·dom may be broad, and supports evidence of persistence for years after oceanographic conditions return closer to the long-term mean [ 53 ]. Similar findings of the lingering effects of the PMHterm·dom have been documented in seabirds in Cook Inlet, Alaska [ 68 ].

Bias correctionterm·key 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 whalesterm·key [ 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 estimationterm·key at the scale of the full North Pacificterm·key 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 Pacificterm·key humpback whaleterm·key 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 estimationterm·key 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 Pacificterm·key abundance estimate, compared to narrower ranges in Hawai‘i and Mexico relative abundanceterm·mth 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 estimationterm·key 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 whalesterm·key 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 abundanceterm·mth 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 abundanceterm·mth between Hawai‘i and Mexico, and possible undetected shifts in relative abundanceterm·mth between other study area regions, doespersName violate this assumption and thus may have introduced error in our geographical bias correctionterm·key factor affecting absolute abundance estimates. We realize that the 13% relative abundanceterm·mth shift is a violation of the assumption of consistent geographical distribution underlying our sampling bias correctionterm·key 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 Pacificterm·key humpback population was captured in most 3-year samples and (ii) the bias correctionsterm·key themselves are relatively small. The overall trends without the bias correctionterm·key 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 whalesterm·key in the North Pacificterm·key should include absolute abundance estimationterm·key by region.

Relative abundanceterm·mth 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 estimationsterm·key 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 abundanceterm·mth 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 whaleterm·key relative abundanceterm·mth off west Maui (which traditionally has hosted the greatest concentrations of humpback whalesterm·key including calves of the year, e.g. [ 77 – 79 ]) was reported after the 2014–2016 PMHterm·dom, 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 abundanceterm·mth 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 PMHterm·dom. Fluctuating but resilient humpback whaleterm·key crude birth rates in Hawai‘i [ 38 , 80 ] suggest the measured −4.3% yr −1 decline from 2014 to 2021date 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 capacityterm·dom triggered by rapid climate changeterm·dom.

Conclusion The end of the industrial whaling era left oceans largely empty of great whales [ 23 , 24 ]. Yet with no species hunted to global extinction during this period, full recovery became a possibility. For humpback whalesterm·key, steady population growth was the 40-year trend in the North Pacificterm·key, and it is to be celebrated that humpback whalesterm·key appeared sufficiently recovered to qualify for partial delisting from legal protected status regimes in the USA and Canada. Ironically, the timing of legal conservation status changes for the species was coincident with the dramatic ocean-warming-induced decline documented in this study. Some populations of humpback whalesterm·key may no longer be a priority for endangered species conservation funding [ 81 ], but now offer high value as an indicator species of ocean-basin-wide ecosystem health [ 82 ] in a world where we can expect increased frequency and severity of marine heatwavesterm·key. Past population monitoring efforts offered periodic abundance estimates interspersed with many years with very coarse assessment of population trajectory ( figure 2 ). This study creates a shift to continuous population monitoring with rapid feedback, enabled by an advance in AI-assisted photo-ID [ 44 ] and an ocean basin-wide collaborative effort [ 18 ]. If this cost-effective collaborative effort is maintained with funding to gather robust sample sizes from most wintering and summer areas of the species, we can expect to be able to detect shifts in abundance in near real time through ongoing updates of the model with data from each passing year.

This study establishes a beginning, and represents just one of many possible applications of a rich and continuous collaborative dataset to monitor humpback whaleterm·key abundance. Aggregation of historical data together with a coordinated future study could enable regional abundance estimationsterm·key, fine-scale understanding of migratory patterns, and time-varying survivorship and mortality rates. Our method of estimating population size was deliberately based on the simplest possible approach (pooling all wintering and summer areas) in order to be comparable to the most recent previous estimate of humpback whaleterm·key abundance for the North Pacificterm·key [ 30 ]. There were substantial limitations to this approach: it did not provide region-specific estimates of abundance or trends, it was a closed-population model and did not allow for estimation of mortality or birth rates, and individual estimates for a given year did not use information available in the longer time-series to improve the estimate for that year. It is our hope and expectation that this is just the first of many explorations of this rich and living dataset. Continued collaborative efforts across North Pacificterm·key wintering and feeding areas have great potential for population monitoring and as an indicator of ecosystem health; some of this is already underway. The In-depth Assessment process of the International Whaling Commission Scientific Committee is currently developing a population model seeking to include a changing carrying capacityterm·dom term that accounts for environmental variability [ 29 ]. And a current initiative in Hawai‘i to fit an integrated population model to these data may refine understanding of regional abundance in the context of post-delisting monitoring under the ESA [ 83 ].

Acknowledgements The Happywhale dataset would not exist without the dedication, support, and trust of a very large community; see electronic supplementary material, II, for full acknowledgements.

Ethics This work did not require ethical approval from a human subject or animal welfare committee.

Data accessibility Data and code are accessible here: https://github.com/tedcheese/RSOS-NPAC-abundance .

Supplementary material is available online [ 84 ].

Declaration of AI use We have used AI-assisted technologies in creating this article.

Authors' contributions T.C.: conceptualization, data curation, formal analysis, funding acquisition, investigation, methodology, project administration, resources, software, supervision, validation, visualization, writing—original draft, writing—review and editing; J.B.: conceptualization, data curation, formal analysis, methodology, software, validation, writing—original draft, writing—review and editing; J.M.A.: data curation, writing—review and editing; K.A.: data curation, funding acquisition, writing—review and editing; L.B.: data curation, writing—review and editing; C.B.: data curation; O.S.B.: data curation; A.L.B.: data curation, writing—review and editing; J.K.B.: data curation, writing—review and editing; J.C.: conceptualization, data curation, funding acquisition, writing—review and editing; R.C.: data curation, writing—review and editing; J.C.: data curation; A.J.G.C.: data curation; J.C.: data curation, funding acquisition, writing—review and editing; R.C.D.C.: data curation; J.D.W.: data curation, writing—review and editing; N.D.: data curation; T.D.-V.: data curation, writing—review and editing; K.D.: data curation; O.A.F.: data curation; R.F.: data curation; K.R.F.: data curation, writing—review and editing; J.K.F.: data curation; A.F.-J.: data curation; C.G.: data curation, writing—review and editing; B.G.: data curation; C.H.: data curation; J.H.: data curation; M.C.H.: data curation, writing—review and editing; J.K.J.: data curation; M.E.J.-L.: data curation; M.J.: data curation; N.K.: data curation; M.L.: data curation, writing—review and editing; E.G.L.: data curation, writing—review and editing; M.M.: data curation; E.G.M.: data curation; P.M.L.: data curation; A.M.: data curation; C.O.M.: data curation; C.M.: data curation, writing—review and editing; J.M.: data curation, funding acquisition, writing—original draft, writing—review and editing; J.M.: data curation; J.L.N.: data curation, writing—review and editing; H.N.: data curation; H.O.: data curation; M.O.: data curation; C.D.O.-O.: data curation; A.A.P.: data curation, writing—review and editing; D.M.P.: data curation, writing—original draft, writing—review and editing; H.C.P.: data curation; E.Q.-R.: data curation, writing—review and editing; R.R.B.: data curation; N.R.: data curation, writing—review and editing; H.R.-N.: data curation; F.S.: data curation; T.J.H.S.: data curation; K.S.: data curation, project administration, software, writing—review and editing; S.S.: data curation, writing—review and editing; I.S.: data curation; J.M.S.: conceptualization, data curation, funding acquisition, writing—review and editing; A.S.: data curation; S.T.: data curation; O.T.: data curation; J.U.-R.: data curation; M.v.A.: data curation; M.V.: data curation; O.v.Z.: data curation; B.W.: data curation; J.W.: data curation; K.M.Y.: data curation, writing—review and editing; I.Y.: formal analysis, software, visualization; D.Z.: data curation; P.C.: conceptualization, supervision, writing—review and editing.

All authors gave final approval for publication and agreed to be held accountable for the work performed therein.

Conflict of interest declaration We declare we have no competing interests.

Funding Direct funding for this work was provided by NOAA Fisheries, the Eric and Wendy Schmidt Family Foundation, CheesemanspersName’ Ecology Safaris, Cascadia Research Collective, the Marine Mammal Center, the Hurtigruten Foundation, Defenders of Wildlife, and many individual supporters of Happywhale.

References Associated Data Data Citations Data Availability Statement Data and code are accessible here: https://github.com/tedcheese/RSOS-NPAC-abundance .

Supplementary material is available online [ 84 ].