Male Humpback Whale Chorusing in Hawai‘i and Its Relationship With Whale Abundance and Density

Kügler, Anke; Lammers, Marc O.; Zang, Eden J.; Pack, Adam A.

doi:10.3389/fmars.2021.735664

ORIGINAL RESEARCH article

Front. Mar. Sci., 23 December 2021

Sec. Ocean Observation

Volume 8 - 2021 | https://doi.org/10.3389/fmars.2021.735664

This article is part of the Research Topic Innovation and Discoveries in Marine Soundscape Research View all 20 articles

Male Humpback Whale Chorusing in Hawai‘i and Its Relationship With Whale Abundance and Density

$\r\nAnke Kügler,*$ Anke Kügler^1,2*

Marc O. Lammers^2,3

Eden J. Zang^2,3 $Adam A. Pack,\r\n$ Adam A. Pack^4,5

¹Hawai‘i Institute of Marine Biology, University of Hawai‘i at Mānoa, Kāne‘ohe, HI, United States
²Oceanwide Science Institute, Honolulu, HI, United States
³Hawaiian Islands Humpback Whale National Marine Sanctuary, Kīhei, HI, United States
⁴Department of Psychology and Biology, University of Hawai‘i at Hilo, Hilo, HI, United States
⁵The Dolphin Institute, Hilo, HI, United States

Passive acoustic monitoring (PAM) with autonomous bottom-moored recorders is widely used to study cetacean occurrence, distribution and behaviors, as it is less affected by factors that limit other observation methods (e.g., vessel, land and aerial-based surveys) such as inclement weather, sighting conditions, or remoteness of study sites. During the winter months in Hawai‘i, humpback whale male song chorusing becomes the predominant contributor to the local soundscape and previous studies showed a strong seasonal pattern, suggesting a correlation with relative whale abundance. However, the relationship between chorusing levels and abundance, including non-singing whales, is still poorly understood. To investigate how accurately acoustic monitoring of singing humpback whales tracks their abundance, and therefore is a viable tool for studying whale ecology and population trends, we collected long-term PAM data from three bottom-moored Ecological Acoustic Recorders off west Maui, Hawaii during the winter and spring months of 2016–2021. We calculated daily medians of root-mean-square sound pressure levels (RMS SPL) of the low frequency acoustic energy (0–1.5 kHz) as a measure of cumulative chorusing intensity. In addition, between December and April we conducted a total of 26 vessel-based line-transect surveys during the 2018/19 through 2020/21 seasons and weekly visual surveys (n = 74) from a land-based station between 2016 and 2020, in which the location of sighted whale pods was determined with a theodolite. Combining the visual and acoustic data, we found a strong positive second-order polynomial correlation between SPLs and abundance (land: 0.72 ≤ R² ≤ 0.75, vessel: 0.81 ≤ R² ≤ 0.85 for three different PAM locations; Generalized Linear Model: p_land ≪ 0.001, p_vessel ≪ 0.001) that was independent from recording location (p_land = 0.23, p_vessel = 0.9880). Our findings demonstrate that PAM is a relatively low-cost, robust complement and alternative for studying and monitoring humpback whales in their breeding grounds that is able to capture small-scale fluctuations during the season and can inform managers about population trends in a timely manner. It also has the potential to be adapted for use in other regions that have previously presented challenges due to their remoteness or other limitations for conducting traditional surveys.

Introduction

Many animals across numerous taxa produce sounds (e.g., invertebrates: Henry and Wells, 2010; Bohnenstiehl et al., 2016; amphibians: Tobias et al., 2010; Dapper et al., 2011; reptiles: Hibbitts et al., 2006; Sicuro et al., 2013; birds: Catchpole and Slater, 2008; Phongkangsananan et al., 2014; fish: Maruska et al., 2012; Tricas and Boyle, 2014; terrestrial mammals: Boughman, 1997; De La Torre and Snowdon, 2009; Van Belle et al., 2014; Pitcher et al., 2015; marine mammals: Smolker and Pepper, 1999; Van Parijs et al., 2000; Rogers and Cato, 2002; Belikov and Bel’kovich, 2007; Sanvito et al., 2007), that can be recorded and analyzed to study their behaviors and ecology. Passive acoustic monitoring (PAM) using autonomous recorders has become increasingly widespread in wildlife research to study animals in terrestrial, aerial, and marine environments (Sugai et al., 2019; Mooney et al., 2020). In the study of cetacean occurrence, distribution, and behaviors, PAM has become an essential tool as it less influenced by factors that often limit visual observation methods (e.g., vessel, land, and aerial-based surveys) such as inclement weather, sighting conditions, the accessibility of study sites, and funding (Zimmer, 2011). PAM projects on various whale and dolphin species range widely from simple detection (e.g., presence/absence in an area) and distribution studies (Van Parijs et al., 2002; Oswald et al., 2011; Soldevilla et al., 2011; Schaffeld et al., 2016) to addressing questions about exposure to potentially disturbing anthropogenic sounds (Lucke et al., 2009; Dunlop et al., 2018). Automated detection and species classification algorithms have become particularly useful for monitoring certain populations (Oswald et al., 2007; Baumgartner and Mussoline, 2011; Baumann-Pickering et al., 2013; Allen et al., 2021), while other species remain problematic. Progress has also been made using PAM to infer density estimates and population trends using animals’ acoustic vocalization rates (Marques et al., 2013; Küsel et al., 2016), for example in beaked whales (Marques et al., 2009; Barlow et al., 2021), pygmy (Kogia breviceps) and dwarf (Kogia sima) sperm whales (Hildebrand et al., 2019), gray whales (Eschrichtius robustus) (Burnham and Duffus, 2020), minke whales (Balaenoptera acutorostrata) (Martin et al., 2013), fin whales (Balaenoptera physalus) (McDonald and Fox, 1999; Harris et al., 2018), and various delphinid species (Van Parijs et al., 2002; Frasier, 2015; Martin et al., 2020). However, using individual’s vocalization rates to estimate densities becomes challenging when an abundance of animals vocalize so often that individual sounds overlap and cannot be easily distinguished in recordings, as is often the case in large aggregations including on the breeding grounds of certain marine mammals, such as the humpback whale (Megaptera novaeangliae) (Au et al., 2000; Seger et al., 2016).

Humpback whales are medium-sized baleen whales that occur globally and migrate annually between high-latitude summer feeding and low-latitude winter calving and mating regions (Clapham, 2018). In the North Pacific, the Hawaiian archipelago is the species’ principal breeding ground, with more than half of the total stock—termed the Hawai‘i Distinct Population Segment (DPS, NOAA, 2016)—wintering in the islands (Herman and Antinoja, 1977; Calambokidis et al., 2008) while the other DPSs migrate to Mexico, Central America, southern Japan, the Philippines, and the Mariana Archipelago (Baker et al., 1986; Calambokidis et al., 2008; NOAA, 2016; Clapham, 2018; Hill et al., 2020).

Humpback whales produce a variety of vocalizations (e.g., Stimpert et al., 2007, 2011; Au and Hastings, 2008). So-called “social sounds” are produced by all ages and sexes and occur on the feeding grounds, along migratory routes, and on the breeding grounds (Silber, 1986; Dunlop et al., 2008; Zoidis et al., 2008; Darling, 2015). In contrast, “song,” known by sailors for centuries, but first formally described by Payne and McVay (1971), is a complex acoustic display produced by mature male humpback whales and some non-calf subadults (Herman et al., 2013). While most prevalent in breeding areas, song may be heard along migratory routes and to a lesser extent in high-latitude feeding areas (e.g., Stimpert et al., 2012; Vu et al., 2012). Song is one of the most recognized and studied features of the species, and because of its prominence during the breeding season, is hypothesized to play an important function in the humpback whale mating system, possibly signaling to females and/or mediating male-male interactions (Herman and Tavolga, 1980; Darling et al., 2006; Smith et al., 2008; Cholewiak et al., 2018; see Herman, 2017 for a comprehensive review on the function of song).

The cacophony of multiple males singing in an ongoing asynchronous chorus becomes the dominant source of low frequency ambient noise in the marine soundscape during the winter and spring months on many breeding grounds (Au et al., 2006; Seger et al., 2016) (Figure 1). In contrast to singing males, females and mother-calf pairs are acoustically cryptic on the breeding grounds (Videsen et al., 2017), male calves and yearlings do not yet sing (Herman et al., 2013), and social sounds in other types of assemblages such as competitive groups tend to be relatively infrequent compared to song (Silber, 1986). Consequently, although non-song vocalizations are also present, PAM studies conducted on humpback whale breeding grounds will predominantly capture the acoustic activity of singing males. In Hawai‘i, seasonal increases and decreases in the amplitude of song chorusing (Au et al., 2000; Chen, 2017; Kügler et al., 2020) mirror the whales’ historically documented seasonal changes in relative abundance (e.g., Baker and Herman, 1981; Mobley et al., 1999). Specifically, (Au et al., 2000) found that the peak in chorusing from singing males corresponds with the peak in abundance obtained from aerial surveys at the same time (Mobley et al., 1999, 2001). These aerial surveys covered all age and sex classes and included areas that are heavily used by mother-calf pairs and other non-singing whales (Mobley et al., 2001). However, the relationship between chorusing levels and overall whale numbers, including non-singing whales, remains poorly defined. In Australia, Noad et al. (2017) found a positive correlation between the number of individual singers and the number of whales sighted on a migration route. Similarly, Seger et al. (2016) showed that daily noise levels correlated positively with counts of whales from visual observations on a Mexican breeding ground. However, whale densities along the Australian migration route and in Mexico differ markedly from the Hawaiian breeding ground, making comparisons difficult. Furthermore, although the recent development in automated detections of song has shown to be a valuable tool to monitor humpback whale occurrence in low-density and transient regions as well as to document seasonal patterns (Allen et al., 2021), it is currently of limited use for estimating abundance until better models are available, particularly in high-density areas such as the four-island region of Maui, Moloka‘i, Lāna‘i, and Kaho‘olawe (also called “Maui Nui”), which, along with Penguin Bank (southeast of Moloka‘i), is considered the core winter humpback whale habitat in the Main Hawaiian Islands (Herman and Antinoja, 1977; Mobley et al., 1999).

FIGURE 1

Figure 1. (A) Long-term spectral average spectrogram between 0 and 5 kHz of 24 h recorded at Olowalu (deep) on 02/15/2020, (B) spectrogram of 30 s in the 0–6 kHz band, illustrating the contribution of humpback whale chorusing to the local marine soundscape. Yellow and red color hues indicate higher levels of acoustic energy.

The objective of the current study was to evaluate the applicability of PAM to establish and monitor humpback whale abundance as a non-invasive and comparatively low-cost alternative or complement to traditional survey methods in the high-density environment off Maui. We sought to define the quantitative relationship between male song chorusing levels and overall whale numbers, as well as to identify drivers impacting this relationship, through the synthesis of acoustic and visual survey data collected over four breeding seasons. Being able to use male chorusing as a proxy for whale abundance can provide a tool that allows the study of a wide range of humpback whale’s ecological questions, ranging from long-term monitoring of population trends to micro-scale habitat use patterns.

Materials and Methods

Data Collection

Acoustic Data

Acoustic data were collected with two deep-water and one shallow-water bottom-moored, autonomous Ecological Acoustic Recorders (EARs) using a Sensor Technology SQ26-01 hydrophone with a sensitivity of −193.5 dB (Lammers et al., 2008) during the humpback whale breeding season (ca. December through April) at three locations off Maui, Hawai‘i (Figure 2). Deployment locations, depths, and recording periods are summarized in Table 1. All sites were within or near the 200 m isobath, a depth within which most humpback whales are traditionally found in Hawai‘i (e.g., Herman and Antinoja, 1977; Mobley et al., 1999) and located on sand (Maui6, Olowalu—deep) or near a coral reef (Olowalu—shallow). The shallow EAR was anchored to a concrete block on the bottom and was diver deployed. The two deep EARs were deployed off a small vessel and were each coupled with a synthetic foam float, two SubSea Sonics AR-60 acoustic releases, and approximately 75 kg of sandbags to anchor the mooring, positioning the EAR into the water column approximately 4 m off the bottom. The EARs recorded at a sample rate of either 25 kHz (Olowalu—shallow) or 50 kHz (Maui6, Olowalu—deep) and a 10% duty cycle, recording for 30 s every 5 min. Acoustic data were analyzed for the period between December 1st and April 30th of each year.

FIGURE 2

Figure 2. The locations of two deep-water Ecological Acoustic Recorders, EARs (dark blue circles) and one shallow-water EAR (light-gray circle) that were deployed between 2016/17 and 2020/21 and the location of the Olowalu land observation station (triangle). (A) The hatched area shows the 15 km field of view from the Olowalu land station. (B) The line shows the transect legs for the vessel surveys in the 2018/19 through 2020/21 seasons where the hatched area is the approximate area covered during the surveys including a 3 km buffer.

TABLE 1

Table 1. Summary of EAR site, location, depth, and bottom type for continuous recording periods during the Hawaiian whale season (approx. December through April) 2016/17 through 2020/21 for two shallow and one deep EARs off Maui, Hawai‘i.

Land-Based Visual Surveys

Throughout each whale season, visual land-based scan samples (Martin and Bateson, 1986) were conducted one to three times each week from a hill near Olowalu on west Maui (20.829, −156.6315, 83.7 m elevation; Figure 2) one to three times each week. The total field of view was scanned for whales once per scan from one end across to the other at a slow and steady rate by the same observer during all surveys, using a pair of Bushnell 7 × 50 marine binoculars. Scans took place for 30 min every hour, up to seven times starting at 08:00 local standard time, weather permitting. The scanning direction was alternated during the first year and randomized during all subsequent years to reduce bias by inadvertently spending unequal amounts of time surveying an area segment at the beginning or the end of each scan. Any whale or group of humpback whales, regardless of number of individuals, was called a pod. When a pod was sighted, the location in vertical and horizontal angles was fixed by the observer with a Lietz DT5 digital surveyor’s theodolite and digitally transferred to a computer upon prompt by the computer operator whilst the observer continued their scan. The geo-referenced program Pythagoras (Gailey and Ortega-Ortiz, 2000) automatically converted the time-stamped angles into latitude and longitude coordinates in real-time and calculated the distance of the fix from the land station. Additional metadata that were entered into Pythagoras by the computer operator with each fix included the immediately observable pod composition (whale, dyad, competitive pod = comp pod, mother-calf pair = M/c, mother-calf-escort = M/c/E) and behavior (blow, surfacing, fluke up dive, surface active, resting/logging) as called out by the observer. Because these were immediate assessments and no prolonged time was spent observing pods to refine their composition and behavior before continuing with a scan, recorded pod composition and behavior were considered conservative estimates, with the presence of a calf likely underreported because of its relatively small size and profile at the surface. In terms of defining pod compositions, all pods of the “whale” category consisted of a single observed blow or tail fluke, but true solitary status could not be unequivocally confirmed. Although a humpback whale competitive group traditionally consists of two or more males competing for access to a single female, for our purposes a comp pod was specified simply if the pod contained at least three adult whales regardless of whether competitive behavior was observed or not. If possible, the true minimum number of whales in a comp pod was noted, but we assume an underestimation of the actual number of whales as this often was not feasible given the survey design and protocols. Dyads consisted of two non-calf whales. Mother-calf pairs were designated based on the sighting of an adult-sized whale and a whale in close proximity approximately 1/3 the size of the adult. If a third whale of adult size was present in this group, the pod was designated as a mother-calf-escort.

If a pod could not be fixed with the theodolite, compass bearing and distance in reticles were determined with the binoculars and noted as a time-stamped “non-fix” with the same composition and behavior information as fixes. During each scan, no pod was recorded twice to avoid pseudo-replication. Scans were conducted in a non-focal-follow manner and regarded as discrete; thus any pod was recorded as a new pod in subsequent scans, even if previously fixed. Environmental data (glare, visibility/haze, Beaufort sea state, wind speed, wind direction) were collected prior to each scan or when conditions changed during a scan. Scans were conducted in conditions of Beaufort sea states of 3 or less (to maximize the sightability of whales including calves; cf. methods described in e.g., Craig et al., 2014), wind speeds at the land site under 15 knots, and dry conditions.

Vessel-Based Visual Surveys

During the 2018/19, 2019/20, and 2020/21 whale seasons, vessel-based surveys following a line-transect design (Figure 2) were conducted using the Hawaiian Islands Humpback Whale National Marine Sanctuary’s R/V Koholā, a 11.6 m long, two-level rigid-hulled inflatable boat. The vessel travelled along the transect lines at a nominal speed of 10 knots (±2 knots) during surveys. Two observers positioned on each side of the flying bridge located 3.5 m above sea level, respectively, scanned the starboard and port sides. The observers actively scanned using unaided eyes from 90° abeam of their side of the vessel to 10° beyond the midline of the bow, resulting in a 100° viewing arc each. Upon sighting of a humpback whale pod, an observer measured the angle between the sighting and the track line using a TrekGuideLITE™ digital compass, obtained the minimum, maximum and best estimated number of animals, and estimated the distance in meters between 0 m and 2,000 m and in 2,000–2,500 m, 2,500–3,000 m, and >3,000 m bins, respectively, for distances greater than 2,000 m. These metrics were reported to a data recorder while the other observer continued to scan and the vessel maintained travel speed and distance. All observers underwent intensive on-water training in distance estimation from the vessel prior to formal surveys. In addition, each survey day, the observers were calibrated by obtaining blind estimates for six randomly chosen distances from a fixed reference buoy whose true distance was logged by the vessel captain. These estimates were subsequently used to obtain calibration curves for each observer and accordingly correct all sighting distance estimates obtained while on effort (Zang and Lammers, 2021). All sighting information as well as the GPS position of the vessel and the time of the sighting were logged into a custom digital form created in the program Poimapper™ on an Apple iPad mini^® by the data recorder. Additional environmental data (sea state, swell size in m, cloud cover, wind speed, glare) were collected at the beginning of each transect leg or when conditions changed.

Port observer, starboard observer, and data recorder rotated every two transect legs to reduce fatigue. Survey effort was conducted in conditions that ranged between Beaufort sea state 0 and 4. If the sea state exceeded 4, effort was paused while the vessel continued to travel along the track line until conditions improved or until the survey was aborted due to unworkable conditions.

Data Processing and Analysis

Acoustic Data

All acoustic recordings were down-sampled to a 3 kHz sample rate using Matlab (MATLAB and Statistics Toolbox Release 2015b, The MathWorks, Inc., Natick, MA, United States), resulting in a 0–1.5 kHz analysis bandwidth, which corresponds to the range where the majority of humpback whale song units have their peak frequencies (Payne and McVay, 1971; Au et al., 2006) and reduces influences from snapping shrimps (Au and Banks, 1998) and environmental sounds such as wind and waves (Urick, 1983). A custom Matlab program converted the 30 s recording from sample values to pressure, calculated the root-mean-square sound pressure levels (RMS SPL), and converted the result into dB re 1 μPa² (following standards described in Ainslie et al., 2018), as

RMS SPL = 20 \log \sqrt{\frac{1}{T} \int_{0}^{T} p^{2} (t) d t},

where T is the duration of each file (30 s) and p(t) is the pressure p re 1 μPa² at time t. Humpback whale chorusing is the dominant source of low-frequency acoustic energy during the whale season (Au et al., 2000; Au and Green, 2000), therefore RMS SPL is a representative metric of the cumulative amount of singing happening during any given 30 s recording period.

All subsequent data processing and statistical analyses were conducted in R 4.0.1 (R Core Team, 2016). The median RMS SPL in dB re 1 μPa² was calculated for each day and recording site. This reduces auto-correlation within the time-series, as recordings are not independent from each other when a song continues over several files. It also reduces the impact of outliers from intermittent high-intensity noise (such as transiting small vessels) and from individual whales singing close to the recorders, temporarily increasing the SPL levels of individual recordings. Daily medians were selected that coincided with concurrent land-based and/or vessel-based survey days (cf. Table 2).

TABLE 2

Table 2. Summary of survey effort for land-based visual (land) and vessel-based line-transect (vessel) surveys off west Maui between 2017 and 2021 and whether acoustic data were obtained from the three acoustic recorders on the same day.

Land- and Vessel-Based Whale-Sighting Data

Compass bearings and reticle distances of all non-fix observations from the land station were converted into latitude and longitude coordinates with an adapted version of the whalemap() function of the bangarang package (Keen, 2016); the function also calculates distance to the observer. All pods with distances greater than 15 km were excluded, to create a field of view that is approximately equal to the area covered by the vessel surveys and extends to the offshore Maui6 EAR (Figure 2). Three days during which only one scan was completed due to poor weather were excluded from further analysis (cf. Table 2).

For all pods, the number of individual whales was calculated depending on their composition as: “whale”: n = 1, “dyad”: n = 2, “M/c”: n = 2, “M/c/E”: n = 3, “comp pod”: n = 3 or whatever number of whales above three was detected upon sighting. The total number of pods and total number of whales were determined for each completed scan per survey day. As noted earlier, the same pod was potentially fixed repeatedly during consecutive 30 min scans of the same area. Therefore, calculating the total number of whales per day could result in pseudo-replication and misrepresentation of the actual number of whales present during the survey. Instead, the daily median number of whales per scan was determined, thereby also standardizing for differences in the number of completed scans among days due to factors such as weather. The proportion of pods with a calf and the proportion of competitive pods were calculated from the total number of pods per day and the daily total number of the respective two pod types.

For each day of vessel surveys, absolute abundance and density estimates were calculated using the Distance package (Miller et al., 2019). A full discussion of the analysis resulting in those estimates can be found in Zang and Lammers (2021). Applying the same protocol as for the land-based data, any pod observed during transects was considered a comp pod when the best estimate for group size was at least three adult whales, and the proportion of competitive pods as well as the proportion of pods with a calf were calculated from the total number of observed pods per survey. If a comp pod contained a calf, it was included in both categories.

Statistical Analysis

Multivariate models were fit to the data to determine the quantitative and qualitative relationship between abundance of whales and chorusing. To select the model that best describes the data following initial visual exploration, non-linear least squares models using the nls() function from the base stats package in R for SPL ∼ whales were fit for a linear, quadratic, logistic, and asymptotic relationship for each recorder site and survey type (cf. Figure 3 for visualizations from Maui6 as a representative example), where “whales” was the median number of whales per 30 min scan for land-based surveys and density in whales per square-kilometer for vessel-based surveys. Bias-corrected second-order Akaike Information Criteria (AICc) for each model per site were calculated using the AICc() function from the MuMIn package (Bartoń, 2020). Pairwise differences between the AICc with the lowest value and all other model AICcs (ΔAICc) were also calculated for each site and survey type. Different models were considered to support the data the same if their ΔAICc was not greater than 2. For all three recorder locations and the vessel-based data as well as Maui6, Olowalu—shallow—and the land-based data, the quadratic model had the lowest AICc values, followed by the logistic and the asymptotic models (Table 3). For the land-station and the deep Olowalu location, the logistic model had the lowest AICc, but the quadratic, logistic, and asymptotic models were all equally supported by the data. The linear model had the least support overall. A second-degree polynomial model was therefore chosen for subsequent analysis.

FIGURE 3

Figure 3. Comparison of four different types of regression fit with a non-linear least squares model for median daily root-mean-squared sound pressure levels (RMS SPL) in dB re 1 μPa² in the 0–1.5 kHz frequency band recorded at the Maui6 EAR site off Maui, Hawai‘i and median number of whales per 30 min scan from a land station.

TABLE 3

Table 3. Evaluation of four different models for the recorder location-specific relationship between daily median RMS SPL (dB re 1 μPa²) and: (a) median number of whales per 30 min scan from land, (b): density estimates (whales/km²) from vessel-based line-transect surveys.

To describe the quantitative relationship between (relative) abundance of whales (whales) and average cumulative chorusing levels (SPL), a generalized linear model (GLM) was fit with the glm() function from the stats package using the equation,

S P L = ß {}_{0}+ ß_{1} \cdot whales + ß_{2} \cdot {whales}^{2} + ß_{3} \cdot whales \cdot Site + ß_{4} \cdot {whales}^{2} \cdot Site + Site + ε,

where ß₀ is the intercept, ß₁ and ß₂ are the respective parameter coefficients, Site is a factor for the three EAR locations, and ε is the error term. Two separate GLMs were fit, one for the land-based data and one for the vessel-based data. In both models, there were no significant interactions between whales and recorder location (land: χ² = 5.65, DF = 4, p = 0.23; vessel: χ² = 0.48, DF = 4, p = 0.98), so a Site interaction was excluded from subsequent modeling.

To explore other factors impacting chorusing levels as well as the relationship between whale numbers and chorusing levels, a generalized additive model (GAM) was fit with the mgcv package (Wood, 2006) with the following equation,

SPL = ß_{0} + f_{1} (whales) + f_{2} (% calves) + f_{3} (% compPods) + f_{4} (whales \cdot day . season) + f_{5} (whales \cdot % calves) + f_{6} (whales \cdot % compPod) + Site + ε,

where %calves is the percentage of pods with a calf per day, %compPods is the percentage of competitive pods per day, and day.season is the count of days into the season for the survey date, with the first day of the season defined as December 1st. For GLMs, Type II likelihood ratios were determined with the Anova() function from the car package (Fox and Weisberg, 2019) and a statistical significance of α = 0.05 was used. For the GAM, likelihood ratios were determined with the summary() function from the mgcv package and using α = 0.05.

While full models including Season main effect and interaction terms suggested that the relationship between whale numbers and chorusing levels significantly differed across the different whale seasons (but not sites), we argue that this is an artifact of survey bias from small sample sizes rather than a true underlying biological difference between years. By pooling our data across years, we increase our overall sample size and therefore the statistical power of our models, and are thus better able to achieve our goal of describing a generalized relationship between whale abundance and RMS SPL.

Results

Recordings during the breeding season were obtained for five consecutive years between 2016/17 and 2020/21 for two recording sites (Maui6 and Olowalu—shallow) and for 4 years at the deep Olowalu EAR site (2017, 2019, 2020, 2021). A total of 74 visual surveys from the Olowalu land station (2017: n = 19, 2018: n = 26, 2019: n = 16, 2020: n = 12) and 26 vessel-based line-transect surveys (2019: n = 8, 2020: n = 8, 2021: n = 10) were conducted. Three land survey days (one in 2018 and two in 2020) were excluded because only one scan was completed, resulting in 71 days with data that were analyzed. The overlapping survey effort for the different methods is summarized in Table 2.

Overall, daily median RMS SPL levels ranged from 88.67 to 114.78 dB re 1 μPa² at Maui6 (Median ± SD in 2017: 104.58 ± 5.20, 2018: 100.80 ± 5.58, 2019: 103.75 ± 6.04, 2020: 105.22 ± 6.75, 2021: 106.63 ± 6.10), from 90.18 to 121.85 dB re 1 μPa² at Olowalu (deep) (Median ± SD in 2017: 109.22 ± 6.43, 2019: 110.48 ± 7.32, 2020: 115.43 ± 8.23, 2021: 113.17 ± 5.93), and from 94.95 to 118.99 dB re 1 μPa² at Olowalu (shallow) (Median ± SD in 2017: 109.86 ± 5.39, 2018: 107.06 ± 4.98, 2019: 109.35 ± 7.10, 2020: 109.48 ± 7.1, 2021: 110.67 ± 6.07). RMS SPL daily medians showed a strong seasonal pattern with an increase in December through January, a peak between February and March, and a decreasing trend through April (Figure 4).

FIGURE 4

Figure 4. Median daily root-mean-squared sound pressure levels (RMS SPLs) in dB re 1 μPa² in the 0–1.5 kHz frequency band, between 1 December and 30 April of 2016 through 2021 at two deep (Maui6, Olowalu—deep) and one shallow (Olowalu—shallow) EAR sites off Maui, Hawai‘i. Green triangles indicate dates of visual surveys from the Olowalu land station; red asterisks indicate dates of vessel-based line-transect surveys.

The median number of whales per 30 min scan from the land station ranged from 2 to 37 (Median ± SD in 2017: 13.50 ± 6.57, 2018: 11.25 ± 7.47, 2019: 14.50 ± 8.59, 2020: 20.25 ± 9.86). Density estimates from the vessel surveys ranged from 0.12 to 1.85 whales/km² (Median ± SD in 2019: 0.53 ± 0.61 whales/km², 2020: 1.07 ± 0.47 whales/km², 2021: 0.62 ± 0.53 whales/km²). Whale numbers from both visual survey methods showed the same seasonal patterns as described for the acoustic data.

Quantitative Relationship Between Chorusing Levels and Whale Numbers

The GLMs showed that the number of whales (land: χ² = 454.27, df = 2, p ≪ 0.001; vessel: χ² = 312.54, df = 2, p ≪ 0.001) and recording location (land: χ² = 134.18, df = 2, p ≪ 0.001; vessel: χ² = 89.59, df = 2, p ≪ 0.001) are significant predictors of RMS SPL (Figure 5). Both the linear and the quadratic term of number of whales are significant (land: χ²_linear = 177.542, p_linear ≪ 0.001, χ²_quadratic = 62.304, p_quadratic ≪ 0.001, vessel: χ²_linear = 138.03, p_linear ≪ 0.001, χ²_quadratic = 67.84, p_quadratic ≪ 0.001), supporting the second-degree polynomial relationship between whale numbers and chorusing levels. Coefficient of determinations, R² ranged from 0.71 to 0.75 for the land-based data (Maui6: R² = 0.71, Olowalu—deep: R² = 0.76, Olowalu—shallow: R² = 0.71) and from 0.81 to 0.85 for the vessel-based data (Maui6: R² = 0.82, Olowalu—deep: R² = 0.81, Olowalu—shallow: R² = 0.85). All parameter estimates from the model and predicted maximum RMS SPL (SPL_max) values are summarized in Table 4.

FIGURE 5

Figure 5. Median daily root-mean-squared sound pressure levels (RMS SPL) in dB re 1 μPa² in the 0–1.5 kHz frequency band recorded at three EAR locations off Maui, Hawai‘i and whale numbers; (A) median number of whales per 30 min scan from a land station (n = 71) between 2017 and 2020, (B) density estimates from vessel-based line-transect surveys (n = 26) between 2019 and 2021. No acoustic monitoring was conducted at Olowalu (deep) in 2018. Lines show the regressions from a second-degree generalized linear model for each EAR site; shaded areas are the respective 95% confidence intervals. Equations indicate parameter estimates and R² values from the models.

TABLE 4

Table 4. Coefficient estimates and their standard errors (SE) from a generalized linear model (GLM) fitting a second-degree polynomial relationship between daily median RMS SPL (dB re 1 μPa²) and (a) median number of whales per 30 min scan observed from a land station and (b) abundance estimates from vessel-based line-transect surveys for three recorder locations.

Whale densities (whales/km²) off West Maui can be estimated for a known RMS SPL value up to ∼SPL_max as:

(1)

N_{M a u i 6} = 1.36 - \sqrt{\frac{S P L - 91.23}{- 10.97} + 1.85} (1)

N_{O l o w a l u (d e e p)} = 1.38 - \sqrt{\frac{S P L - 97.35}{- 13.00} + 1.69} (2)

N_{O l o w a l u (s h a l l o w)} = 1.40 - \sqrt{\frac{S P L - 95.16}{- 11.74} + 1.84} (3)

where SPL_max is the SPL value at the maximum of the quadratic regression curve. Alternatively, density estimates that include a standard error as a measure of uncertainty can be obtained from the model fit with glm() using the invest() function from the investr package in R (Greenwell and Schubert Kabban, 2014).

Qualitative Relationship Between Chorusing Levels and Whale Numbers

When adding additional covariates in a GAM framework for the land-based data, number of whales remains a significant predictor of SPL (F = 106.529, edf = 1.969, p ≪ 0.001), but this relationship depends on the time of the season (F = 4.968, edf = 15.13, p ≪ 0.001). Low whale numbers at the beginning of the season produced comparatively lower SPLs than toward the end, while during the peak of the season the higher whale densities resulted in slightly lowered chorusing levels (Figure 6). The proportion of pods with a calf is not significantly related to SPL levels (F = 1.795, edf = 4.786, p = 0.14), but chorusing levels appear to increase with increasing proportion of calves after a certain point. The proportion of calves does, however, significantly influence the relationship between SPL and whale numbers (F = 4.155, edf = 5.909, p = 0.0003). At low whale numbers, a high proportion of present calf pods is related to a decrease in chorusing. In contrast, the proportion of comp pods is significantly related with chorusing levels (F = 5.716, edf = 2.645, p = 0.0008). SPLs initially decreased with increasing proportion of competitive pods while increasing again after a midpoint with further increasing competitive pods. This relationship is independent of whale numbers (F = 0.910, edf = 2.0, p = 0.42). However, the heatmap (Figure 6) suggests that at very low whale numbers a higher proportion of competitive pods is positively related to RMS SPL levels. Overall variation explained by the model is 91.5%. Visualizations of all partial effects are summarized in Figure 6.

FIGURE 6

Figure 6. Plots of the partial effects of predictors fit with a generalized additive model (GAM) for daily median root-mean-squared sound pressure levels (RMS SPL) in dB re 1 μPa² at three EAR locations off Maui, Hawai‘i, and median number of whales per 30 min from the land station, the proportion of pods with a calf per day, the proportion of competitive pods per day, and the interactions of relative whale numbers and day of the season (where day₁ = December 1), relative whale numbers and proportion of pods with a calf, and relative whale numbers and proportion of competitive pods. Shaded areas indicate 95% confidence intervals. The colors in the heatmaps indicate the relative response of RMS SPL as a result of the partial effect of the respective two predictor variables with darker colors corresponding with a higher response. The isobaths specify areas of 1-unit difference in the response. White areas in the heatmaps reflect specific parameter combinations for which no data was available.

Results from the GAM for the vessel-based data are similar and explain 94.2% of deviance; density (F = 46.022, edf = 1.97, p ≪ 0.001), proportion of pods with a calf (F = 10.568, edf = 1, p = 0.002), but not proportion of comp pods (F = 3.212, edf = 1.854, p = 0.064) are significant predictors of RMS SPL and there are significant interactions between density and day of the season (F = 7.184, edf = 2.866, p ≪ 0.001), and density and proportion of comp pods (F = 3.644, edf = 1.802, p = 0.046), but not density and proportion of calf pods (F = 2.898, edf = 1.824, p = 0.102). However, we urge caution with the vessel-based results due to a relatively large number of terms included in the model compared to the small sample size, and interpretation of interaction effects is limited.

Discussion

Only adult and some juvenile male humpback whales sing (Herman et al., 2013), but their increasing and decreasing contribution to the marine soundscape in Hawai‘i during the breeding season mirrors the bell-shaped abundance curve resulting from the species’ staggered, age-class and sex-segregated migration pattern (Baker and Herman, 1981; Mobley et al., 1999; Au et al., 2000; Craig et al., 2001, 2003; Chen, 2017; Kügler et al., 2020). Our findings confirm that male chorusing intensity is predictive of overall whale numbers including non-singing animals off west Maui, Hawai‘i.

Both visual observation methods used to obtain whale numbers (whale counts from a shore station and density estimates from line-transect vessel surveys) had the same second-order polynomial correlation with acoustic recordings obtained from the three PAM sites within the study area. Observed differences among recording sites are likely the result of sound propagation effects at different recorder depths that shifted the SPL baseline, but not the general qualitative relationship between chorusing levels and whale numbers. These findings are in agreement with those reported along the southbound migration route off East Australia that also showed a positive polynomial correlation of both the number of singers and song minutes with the number of whales sighted (Noad et al., 2017). They also agree with a study conducted on a breeding ground off Mexico that showed a positive correlation between the number of whales and daily sound levels (Seger et al., 2016).

While some males may sing while escorting a mother-calf pair, or as part of a dyad or larger group (Baker and Herman, 1984; Darling and Bérubé, 2001; Smith et al., 2008; Chen, 2017), most singers are lone males (Darling et al., 2006; Herman et al., 2013) that can stay submerged for up to 20 min at a time while singing and may repeat their song in sessions that last many hours. It was therefore difficult to unambiguously identify singers with our visual observation methods, so the proportion of singers in our population is unknown. While Noad et al. (2017) showed a strong positive correlation between the number of singers and the overall number of whales observed along an Australian migratory route, they found that the proportion of singers decreased as the population increased over the course of almost 20 years. This contrasts with results from Cato et al. (2002), who showed that the proportion of singers remained consistent between 1981 and 1994 in that same population. Noad and colleagues considered that the observed decrease in the proportion of singers may have been caused by singing whales moving further offshore and outside of their in-shore study area, but concluded this was unlikely based on their observational data. They also reject the possibility of a decrease in numbers of mature whales resulting from a potential change in age structure. Instead, they proposed that the increase in population size during their recovery from commercial whaling likely changed whale singing behavior, resulting in proportionally fewer males singing. However, singing behavior during migration may differ markedly from that occurring on the breeding grounds. In fact, Seger et al. (2016) found that the proportion of singers in the population remained constant on the Mexican breeding ground. We have no reason to suspect that an abundance-dependency in the proportion of singers as observed on the Australian migration route exists on the Hawaiian breeding ground. In a recent study on population fluctuations off Maui inferred from chorusing levels as a proxy for relative abundance, Kügler et al. (2020) reported a decrease in RMS SPL levels during the month of February, the peak of the whale season, of 6.1 dB between 2014/15 and 2017/18 at the shallow Olowalu EAR site. A −6 dB change in RMS SPLs represents a 50% decrease in acoustic energy. Applying the regression presented here, we calculate that a decrease occurred from 1.21 whales/km² in 2015 to 0.62 whales/km² in 2018, or a change of −48.85%. Cartwright et al. (2019) reported a decline in visually observed mother-calf pairs of 76.5% and of 39% of adult groups between 2013 and 2018 in the same study area off west Maui. In addition, researchers off Hawai‘i Island reported a 50% decline in the mean number of whales counted from a shore station between 2014 and 2018 (NOAA, 2019; Frankel et al., 2021). The similarity in magnitudes of the trends described among these studies support our conclusion that variations in chorusing levels are the result of changes in whale abundance and not changing proportions of singers. Thus, SPL appears to be a reliable indicator of whale numbers over time, although further examination of the relationship between whale numbers and chorusing levels could help to define the relationship and possible limitations resulting from social dynamics more precisely.

Chorusing levels do not increase indefinitely with increasing whale numbers, but rather reach a plateau (logistic and asymptotic regressions) or even start decreasing (quadratic regression) after whale densities reach a threshold. At peak densities in February and early March, opportunities to find potential mates are abundant, which often involves physical contest competition with other males and may result in a lower proportion of males singing (Figure 6). Further, multiple past studies point to song playing an important role as a spacing mechanism among singers (Winn and Winn, 1978; Tyack, 1981; Frankel et al., 1995; Seger et al., 2016). As local densities become high, singers may not be able to maintain their minimum distance from each other, and individual males may choose to switch their behavioral strategy, for example by joining mother-calf pairs (Herman and Antinoja, 1977; Craig et al., 2002) or competitive pods (Au et al., 2000; Darling and Bérubé, 2001; Craig et al., 2002; Pack et al., 2012) rather than sing. At higher densities, singers may also be joined more frequently by other males, which often leads to cessation of singing (Darling et al., 2006).

In general, as the proportion of competitive pods increases, chorusing levels decrease (Figure 6). However, toward the end of the season when there is a high prevalence of competitive pods as the availability of receptive females becomes low and males compete over the last opportunities to mate (Craig et al., 2002), the relative amount of chorusing increases again. This also may indicate a change in male behavioral strategy. Humpback whale song has been proposed to provide information about a male’s competitive fitness (Lammers et al., 2017) and during high competitive activity males may choose to invest additional time separately advertising their fitness to competitors and/or females. Similarly, our data suggest that a high relative presence of pods with a calf is tied to increased chorusing levels, except at low whale densities (Figure 6). The number of calf pods relative to the total whale abundance increases as the whale season progresses and is highest in late March through April (Craig et al., 2003). As females without a calf become rare, the increase in chorusing with the higher proportion of calf pods may reflect an increase in the number of escorts singing in mother/calf groups (e.g., Herman and Tavolga, 1980; Chen, 2017). On the other hand, the decrease in SPL at low whale densities, but high prevalence of calf pods probably results from differences in migratory timing, as mother-calf pairs are among the last to depart the islands for their northward journey (Craig et al., 2003). It is therefore likely that this observed correlation at low densities is an artifact of co-varying seasonal changes in SPL rather than a causal relationship between SPL and calf presence.

The local presence of whales in Hawai‘i fluctuates throughout the breeding season as a function of migratory trends, residency times, and movement patterns within the archipelago (Craig and Herman, 1997; Cerchio et al., 1998; Craig et al., 2001, 2003), a pattern that is reflected in daily chorusing levels (cf. Figure 4). Visual surveys traditionally conducted to estimate abundances are often limited to a relatively large temporal resolution, as they are constrained by feasibility (access to vessels, aircraft, land observation sites or weather), logistics, and funding (e.g., Mobley et al., 1999; Barlow et al., 2011). While these efforts are important for monitoring population trends, they provide “snapshots” and are often not able to capture fluctuations on smaller temporal scales. PAM allows the ability to survey continuously throughout the breeding season and therefore capture naturally occurring fluctuations to gain a better understanding of patterns of abundance. In the future, real-time acoustic monitoring could allow monitoring the local population without a time delay from off-season instrument recovery and data processing, relaying the current status of humpback whales in Hawai‘i to managers and the public.

The Hawai‘i humpback whale DPS is considered recovered from commercial whaling depletion following decades of population growth at a rate of ca. 6–7% (Mobley et al., 1999; Calambokidis et al., 2001), and, together with eight other DPSs, was delisted from the United States Endangered Species Act in 2016 (NOAA, 2016). However, population fluctuations that have occurred since then (Cartwright et al., 2019; NOAA, 2019; Kügler et al., 2020; Frankel et al., 2021) emphasize the importance of continued long-term monitoring and the necessity for tracking population trends in order to mitigate changes in the conservation status of populations in a timely manner. This becomes particularly important in light of accelerating global changes, which are expected to put many marine species under increasing stress through the combination of climate change as well as other anthropogenic impacts, including ship strikes, entanglement, tourism, and noise pollution (Moore, 2008; Lammers et al., 2013; Senigaglia et al., 2016; Tsujii et al., 2018; Cartwright et al., 2019; Guazzo et al., 2020; Currie et al., 2021). Although our results are specific to the west Maui area, they could be adapted for studying populations in other breeding areas that are subject to limited resources or access. For example, the Northwestern Hawaiian Islands, stretching 1,900 km northwest of the Main Hawaiian Island, are a significantly understudied, but presumed important humpback whale breeding habitat (Johnston et al., 2007; Calambokidis et al., 2008; Lammers et al., 2011). Access to the region is restricted due to its remoteness and protection as part of the Papahānaumokuākea Marine National Monument, and frequent visual surveys are not possible. Further refinement of the results presented here to account for different habitats and propagation conditions could help make them more broadly applicable, thus providing a valuable tool for monitoring humpback whale abundance in more remote parts of the archipelago. However, higher wind and/or vessel noise, lower whale densities, differences in behavior and/or migratory patterns could cause unique conditions that deviate significantly from the relationship between whale numbers and chorusing levels presented here.

Finally, this study demonstrates that PAM methods can be extended to estimate population sizes for species and environments where vocalizations are too abundant to identify individual calls. Our results warrant further studies on the applicability of this cost-efficient, non-invasive method to study humpback whale populations on other breeding grounds, as well as other species that provide similar challenges.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

AK performed the statistical analysis and wrote the first draft of the manuscript. All authors contributed to conception, design, and data collection of the study, manuscript revision, read, and approved the submitted version.

Funding

This work was supported by the SanctSound project, which is a collaboration between NOAA and the United States Navy to better understand underwater sound within the National Marine Sanctuary System. For more information, visit https://sanctuaries.noaa.gov/science/monitoring/sound/. Additional funding or in-kind support for this work was provided by Whale Tales/Whale Trust, the Linda and Jim Collister Fellowship, Ultimate Whale Watch, and Oceanwide Science Institute.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

The efforts described would not have been possible without the dedicated help of many individuals and we would like to thank them for donating their personal time to help with EAR recoveries and re-deployments, as computer operators during land-based surveys, and as observers and data recorders during vessel surveys. These volunteers include (in alphabetical order): Nanette “Chandra” Bowen, Lacey Brandt, Isabella Canepa, Sue Canfield, Christina Cantellops, William Carrier, Debra Caswell, Steve Chopek, Peter Colombo, Mark Deakos, Brad Eichhorst, Amy Eldredge, Tim Goeken, Ted Gruppenhoff, Russell Hansen, Lee James, Jeff Kallenbach, Jessie Kittel, Dani Kleinhenz, Matt Koller, Kirtley Mann, Laura Meades, Jason Moore, Marcus Richter, Jeannine Rossa, Julie Strong, Heather Spillane, Scott Splean, Liz Stahl, Larissa Treese, Marina Treese, Sara Wood, and Joan Zwar. We would also like to thank two reviewers for their insightful suggestions on the manuscript.

References

Ainslie, M. A., Miksis-Olds, J. L., Martin, B., Heaney, K., De Jong, C. A. F., Von Benda-Beckmann, A. M., et al. (2018). ADEON Underwater Soundscape and Modeling Metadata Standard. Version 1.0. Technical Report by JASCO Applied Sciences for ADEON Prime. Tokyo: JASCO.

Google Scholar

Allen, A. N., Harvey, M., Harrell, L., Jansen, A., Merkens, K. P., Wall, C. C., et al. (2021). A convolutional neural network for automated detection of humpback whale song in a diverse, long-term passive acoustic dataset. Front. Mar. Sci. 8:607321. doi: 10.3389/fmars.2021.607321