An Alternative Approach to Studying Humpback Whale Song Changes Over Time Using Occupancy Models

Abstract

Male humpback whales (Megaptera novaeangliae) are known for producing complex and structured songs, primarily at their breeding grounds. These songs are hierarchical in structure: units combine to form phrases, phrases combine to form themes, and themes together make up a complete song (Payne and McVay 1971). Several studies have shown that song elements' presence and acoustic structure gradually change over time, likely influenced by the movement of whales as they diverge into different breeding areas at low latitudes and later converge at high latitude feeding grounds (Payne and Payne 1985; Rasmussen et al. 2007; Acevedo et al. 2017; Herman 2017). Occasional extensive changes in song structure over short periods of time (as little as 2 years), known as “song revolutions,” are thought to occur when individuals from one population encounter the song of another and incorporate the new song into the population (Garland et al. 2013; Schulze et al. 2022). This phenomenon is solely reported in the South Pacific, where the circumpolar cultural transmission of whale song travels from eastern Australia to French Polynesia (Garland et al. 2011, 2013) and to South America (Schulze et al. 2022; Gonçalves et al. 2023). These rapid changes in song structure may result from individuals shifting reproductive areas within or between breeding seasons (Schulze et al. 2022) or by shifting environmental conditions—such as reductions in the Antarctic Sea ice—that increase the likelihood of interpopulation encounters during migration or on shared feeding grounds (Andrews-Goff et al. 2018; Schall et al. 2022; Bedriñana-Romano et al. 2022; Schulze et al. 2022; Bengtson Nash et al. 2023).

Whether gradual or rapid, these changes highlight the complexity of humpback whale song structure and underscore the challenges involved in tracking and measuring song variation over time. Changes in humpback whale songs are traditionally measured tracking the presence of themes or phrases over time and performing network and similarity indices to assess the rate and direction of song change (Garland et al. 2011; Allen et al. 2019) and by using Markovian transition probabilities to identify subtle changes in song structure (Allen et al. 2019; Schulze et al. 2022). However, these approaches typically assume perfect detection of produced song elements; they do not account for song elements that were produced but undetected by researchers. Detection may be influenced by the environmental conditions at the time of recording, the duration of the recording itself and of the song element, and the frequency spectrum of the song element in question, among other variables (Harris et al. 2025). Failure to account for imperfect detection can mask our understanding of whale song evolution and hinder the identification of management strategies that minimize disruption to whale communication. The goal of this study was to evaluate the potential application of occupancy modeling as a complementary approach to study humpback whale song change analyses, allowing for determination of the most influential ecological and acoustic factors in determining song element presence and detectability.

Occupancy modeling has been widely used to study the distribution of single species across sites while accounting for “imperfect” detection, when a species might not be detected during a survey even if it is present (MacKenzie et al. 2017). In the single season model (Figure 1A; MacKenzie et al. 2002), each site must be surveyed multiple times to estimate the occupancy probability of a target species (Ψ, the probability of presence) and detection probability (p, the probability that the target species is detected given presence). The occupancy framework has been extended to investigate community dynamics (Dorazio et al. 2006; Guillera-Arroita et al. 2019) and can be used to explore dynamics through time, either directly through the multi-season model (MacKenzie et al. 2003) or indirectly by incorporating time period as a covariate in the single season model (e.g., Twining et al. 2024). Regardless of the approach, each of the occupancy parameters can be functions of covariates that influence the state of the species (present or absent) and the detection of that state (detected or not).

Here, as proof of concept, we demonstrate the potential application of the multi-species occupancy framework to estimate changes in humpback whale song elements through time by treating song units as “species” and tracking the detection of song units across surveys in 4 time periods, enabling inferences about humpback whale song change when detection probability is less than 1 (Figure 1B).

We applied the single-season, multi-species occupancy model framework to a database of coded song units from 17 audio recordings obtained between 2007 and 2023 (July and August) from breeding humpback whales of the Breeding Stock G (BSG) population in the Gulf of Chiriqui in Panama. This audio database consisted of opportunistic audio recordings made during boat surveys for humpback whales with recordings made primarily to detect singers (Figure 2). Wind and weather conditions were partially controlled, as recordings were only obtained during good sea conditions (clear weather, Beaufort of less than 3). The audio recordings were made using an M-Audio MicroTrack 24/96 recorder for recordings made between 2007 and 2010 and a Zoom H4n recorder for all other recordings, and various models of hydrophones (Cetacean Research Technology SQ26-08 with a 10 m cable, an effective sensitivity of −194 dB, re 1 V/uPa, and frequency responses of 30 Hz to 30 kHz for recordings made between 2007 and 2013 and 2 Hz to 50 kHz for recordings made in 2021; High Tech Inc. HTI-96-MN with a sensitivity of −164 dB, re 1 V/uPa, a frequency response of 2 Hz–30 kHz, and a 10 m cable for recordings made between 2014 and 2017; Cetacean Research Technology C57 cylindrical hydrophone with an integral 20 dB preamplifier, underwater connector, effective sensitivity of −167 dB re 1 V/uPa, and a 15 m cable for recordings made between 2018 and 2019), all at a sampling rate of 48 kHz and with the boat engine off. The hydrophones were held over the side of the boat at an approximate depth of 10 m. The length of the recording varied from 3 to 30 min. Recordings that included boat noise or the presence of multiple whales were excluded to minimize acoustic interference.

Acoustic files were visually inspected in RavenPro 1.5 build 37 (Conservation Bioacoustics at the Cornell Lab of Ornithology 2025), and a spectrogram was generated with a fast Fourier transform size of 2048 points, an overlap of 50%, and a 2048-sample Hann window to find and categorize song units. Individual units were categorized based on signal type, contour shape, presence of harmonics, duration, and peak frequency according to a modification of the criteria detailed in Djokic (2021), and every unit was identifiable. In total, we characterized 56 unit types across 17 recordings.

To test the multi-species occupancy modeling approach to examine how issuance of song units change through time (Figure 1B), recordings were grouped into periods (period 1 = 2007–2009; period 2 = 2013–2015; period 3 = 2017–2019; period 4 = 2023; Figure 3A). Further, to meet the requirement for survey replication, an essential component for occupancy modeling (MacKenzie et al. 2002), each recording was divided into two equal segments, hereby referred to as secondary periods, with each segment treated as a survey (Figure 1B). The detection/non-detection of each of the total 56 units was determined for each of the 17 recordings in each of the surveys, creating the required encounter history input for occupancy modeling (n = 17 * 56 = 952 rows in the encounter history).

We used the single-season, multi-species occupancy framework to evaluate 28 alternative models that best explained the observed encounter histories of song units through time, where each model estimated occupancy (psi, Ψ) for each period uniquely and estimated detection (p) for each survey and period uniquely. Additional covariates associated with psi included combinations of mean duration, tone type, and mean peak frequency of song units (Table 1). Additional covariates associated with detection included recorded sea surface temperature (SST), length of recording, day of the year (DOY) (Table 1), all of which have been shown to influence the communication range of humpback whales (Clark 1982; Perazio et al. 2018). Variables were scaled by dividing by a constant to improve numerical stability. We used AIC model selection to evaluate which model best explained the observed encounter histories of song units through time (Akaike 1998). Occupancy model analysis was conducted in RStudio using the RPresence package (MacKenzie and Hines 2024).

Of the 28 models evaluated, one model had overwhelming support with over 99% of the model weight (Table S1). The top model identified mean peak frequency, song unit duration, and tone type as important variables shaping the issuance of song units in each period (Table 2). Parametric bootstrap-based model fit (n = 500; p = 0.19) suggested the model fit the observed encounter histories.

In general, the longer song units with higher peak frequencies dominated the unit repertoire in periods 1–2 (2007–2015), but this pattern relaxed in period 3 (2017–2019), where longer units had a higher probability of being issued, especially tonal units (Figure 4). In period 4 (2023), only those units with the lowest peak frequencies (< 2.5 Hz/1000) were not issued; this period was characterized by a high unit richness, with both pulsed and tonal units being issued across different unit lengths (Figures 2B and 4).

In terms of detection probability, the probability of detecting a song unit given it was issued was strongly related to recording duration (Figure 5).

These results illustrate the potential of occupancy models to investigate the pace of change in humpback whale song, using song units as the basis for analysis and accommodating datasets with variable recording lengths and efforts. By applying occupancy modeling to acoustic recordings with uneven durations and temporal coverage, we show that this framework can reveal how environmental factors and intrinsic acoustic characteristics influence the issuance of song units and their probabilities of detection. Therefore, researchers can use this information to improve their own sampling designs.

Several caveats are in order. First, given the small sample size of this study, we hesitate to rely on the model to make inferences regarding song unit evolution and change in this population. Second, we used the single-season occupancy modeling approach to investigate changes in song units through time. In our approach, we identified 56 song units that could have been issued at any time throughout the study under the assumption that all song units were available throughout the study (e.g., no new units arose during the study). When these assumptions cannot be met, capture-mark-recapture analyses may be a possible alternative to occupancy modeling, treating each song unit as an “individual” to study song dynamics that include addition of new units, and permanent extinction of old ones, over time (e.g., Pradel 1996, Telenský et al. 2024). In addition to the standard single season occupancy model, we explored several other alternative approaches, including the multi-season occupancy model (MacKenzie et al. 2003), the correlated detections model (Hines et al. 2010), and the heterogeneous detections model (Royle 2006); the latter methods were an attempt to control for correlations among surveys that arose from splitting a recording into 2 surveys. All of these approaches yielded the same general result. The multi-season version is a powerful approach (MacKenzie et al. 2003), especially if repeated surveys are truly independent, allowing direct estimation of factors related to unit extinction (the probability that a previously issued song unit in time t is not issued in time t + 1) and unit colonization (the probability that a song unit that was not issued in time t is issued in time t + 1). Third, all covariates included in this analysis to study the probability of detection were associated with the recording process and their corresponding environmental conditions (sea surface temperature and day of year in which the recording was collected). However, we were unable to account for the potential effects of hydrophone settings due to missing metadata. In addition, all unit detection and categorization were carried out manually in Raven Pro (2025), which may introduce observer-related biases. Future studies should consider incorporating both hydrophone settings and observer effects as covariates when modeling detection probabilities.

A benefit of occupancy modeling framework is its adaptability and evolution. With the ability to collect a vast number of recordings, researchers are turning to automated detectors and classifiers for humpback songs based on deep convolutional neural networks (CNNs), which have demonstrated high performance for some populations but are more likely to introduce false positives (Miller et al. 2011; Allen et al. 2021; Kather et al. 2024). Relatively new versions of occupancy models are able to combine these automated detections with human-made labels to mitigate this challenge and facilitate the evaluation of song change through time with consideration paid to detection probabilities (e.g., Miller et al. 2011).

This study highlights the value of occupancy modeling as a powerful framework for examining temporal dynamics in humpback whale song structure. By explicitly modeling detection probability among song components, sampling limitations are accounted for and their effects are distinguishable from true changes in song structure over time. Occupancy modeling may be used in conjunction with similarity- or network-based analyses to identify factors impacting detection probabilities and weigh predicted results. For example, our findings indicate that environmental factors during data collection such as recording duration may play important roles in unit detectability; such information enables researchers to disentangle the contributions of environmental drivers and intrinsic acoustic traits to song structure, as well as understand any shifts in song structure over time as a function of detectability. Finally, this framework has important applications for assessing the influence of anthropogenic noise in humpback whale habitats, which is an increasingly prominent factor in studies of cultural transmission.

This study was funded by Panacetacea, K.R., The University of Vermont via the Kay, Klieman, Larabee Undergraduate Research Award and the John Wheeler Graduate Student Research and Development Award to F.O., and by the NSF Awards #2335991 and #2246323 to L.J.M.-C. Special thanks to Dr. Matt Leslie and Luis “Pulo” Bernal for their support during fieldwork and Divna Djokic for access to her unit characterization protocol and edits to this manuscript. We thank the anonymous reviewer who provided helpful comments that greatly improved the manuscript. Panamanian Research permits # ARBG-0074. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the US Government. The Vermont Cooperative Fish and Wildlife Research Unit is jointly supported by the US Geological Survey, University of Vermont, Vermont Department of Fish and Wildlife, US Fish and Wildlife Service, and Wildlife Management Institute.

Therese Donovan: conceptualization, methodology, writing – review and editing, supervision, formal analysis, visualization, resources. Franny Oppenheimer: conceptualization, investigation, funding acquisition, methodology, writing – original draft, visualization, writing – review and editing, resources, formal analysis. James Hines: methodology, visualization, writing – review and editing, formal analysis. Laura J. May-Collado: conceptualization, data curation, supervision, resources, investigation, funding acquisition, writing – original draft, writing – review and editing, methodology. Kristin Rasmussen: data curation, resources, writing – review and editing, funding acquisition.

This work was supported by Panacetacea; National Science Foundation, Award #2335991, Award #2246323; University of Vermont.

The authors declare no conflicts of interest.

Data and scripts materials are available from the first and corresponding author upon request.