Humpback whale song occurrence reflects ecosystem variability in feeding and migratory habitat of the northeast Pacific

Summary

This study examines the occurrence of humpback whale (Megaptera novaeangliae) song in the northeast Pacific from three years of continuous recordings off central California (36.713°N, 122.186°W). Song is prevalent in this feeding and migratory habitat, spanning nine months of the year (September–May), peaking in winter (November–January), and reaching a maximum of 86% temporal coverage (during November 2017). From the rise of song in fall through the end of peak occurrence in winter, song length increases significantly from month to month. The seasonal peak in song coincides with the seasonal trough in day length and sighting-based evidence of whales leaving Monterey Bay, consistent with seasonal migration. During the seasonal song peak, diel variation shows maximum occurrence at night (69% of the time), decreasing during dawn and dusk (52%), and further decreasing with increasing solar elevation during the day, reaching a minimum near solar noon (30%). Song occurrence increased 44% and 55% between successive years. Sighting data within the acoustic detection range of the hydrophone indicate that variation in local population density was an unlikely cause of this large interannual variation. Hydrographic data and modeling of acoustic transmission indicate that changes in neither habitat occupancy nor acoustic transmission were probable causes. Conversely, the positive interannual trend in song paralleled major ecosystem variations, including similarly large positive trends in wind-driven upwelling, primary productivity, and krill abundance. Further, the lowest song occurrence during the first year coincided with anomalously warm ocean temperatures and an extremely toxic harmful algal bloom that affected whales and other marine mammals in the region. These major ecosystem variations may have influenced the health and behavior of humpback whales during the study period.

Keywords:

Biology and life sciences – Organisms – Eukaryota – Physical sciences – Psychology – Animals – Social sciences – Vertebrates – Amniotes – Mammals – Zoology – Physics – Behavior – Earth sciences – Ecology and environmental sciences – Ecology – Ecosystems – Animal behavior – Marine biology – Marine and aquatic sciences – Habitats – Seasons – Bioacoustics – Acoustics – Foraging – Animal sexual behavior – Marine mammals – Whales – Humpback whales – Animal migration – Seasonal variations

Introduction

Since first being characterized as song nearly a half century ago [1], the cyclic and hierarchically structured vocalization sequences of male humpback whales have been an important subject in behavioral ecology research [2]. The high degree of complexity and variation in humpback song has motivated investigation of its functional roles and temporal evolution. Studies have found evidence that humpback song may function in both intersexual and intrasexual communication [3–5], and it has been hypothesized that song maintains social contact during migration [6] and enables echoic perception [7]. Beyond potential immediate functions of song, its evolution across large spatial and long temporal scales reveals its role in cultural transmission within and between humpback whale populations [8–10].

Early interpretations that humpback song occurs only in breeding habitats have yielded to a growing body of observations of humpback song occurring in feeding habitats [11–19] and during migration between breeding and feeding habitats [6,20,21]. Use of acoustic recording tags that also characterize dive behavior have further shown that song can occur in close spatiotemporal proximity to feeding behavior [16]. This pervasive and complex vocal behavior is essential to humpback whale ecology, yet its occurrence in relation to ecosystem variations in feeding and migratory habitats remains largely unexplored. This study examines occurrence of humpback whale song in a feeding and migratory habitat of the northeast Pacific.

The study region is the Monterey Bay National Marine Sanctuary (MBNMS), in the central California Current System (CCS; Fig 1A and 1B). This region comprises essential feeding and migratory habitat for two of the fifteen distinct population segments (DPS) of humpback whales identified globally [22]. DPS are named according to primary breeding habitat, however their definition is based on a variety of information including estimates of exchange between breeding areas, geographic patterns of habitat occupancy, distinctions in the ecology of feeding and breeding habitats, and genetic differentiation. Studies examining maternally inherited mitochondrial DNA in these populations indicate natal fidelity to breeding grounds and maternally determined fidelity to feeding grounds [23–25]. Humpback whales that feed off the western United States migrate to breeding grounds off southern Baja, mainland Mexico and Central America (Fig 1A) [26–30]. The basis of this study is a three-year time series (August 2015 through July 2018) of continuous passive acoustic recordings, enabled by a cabled observatory (Fig 1B). The Monterey Accelerated Research System (MARS) cabled observatory sits on the continental slope outside Monterey Bay, California. The cable reaches shore at Moss Landing, the site of a former whaling station [31]. Humpback whale song is prominent in the local soundscape (Fig 1C).

Favorable foraging habitat for baleen whales in the CCS results ultimately from high levels of primary productivity caused by wind-driven upwelling of nutrient-rich water into the euphotic zone [32,33]. The physical signature of coastal upwelling is evident as a cooling gradient between the North Pacific subtropical gyre and the eastern boundary (Fig 1A). Eastern boundary upwelling ecosystems like the CCS are among the most productive of ocean margin habitats [32,34–36]. Beyond elevated primary productivity, favorable foraging habitat for baleen whales results from efficient trophic transfer. This efficiency is partly due to short food chains [32]. Between phytoplankton and humpback whales, the few trophic links are through euphausiids (hereafter ‘krill’) and small schooling fish that consume phytoplankton and microzooplankton grazers. Efficiency is also enhanced through development of foraging hotspots that host dense concentrations of plankton, micronekton, and fish [37–39].

Submarine canyons are particularly important foraging hotspots for krill in the CCS [39]. Monterey Bay coincides with Monterey Canyon, the largest submarine canyon along the western US (Fig 1B). The highest occurrences of krill hotspots in the CCS are associated with canyons between Monterey Bay and Bodega Bay, approximately 200 km further north [39]. Within this region, other forage species for humpback whales exhibit a latitudinal gradient, with anchovy and sardine being more abundant closer to Monterey Bay and maximum abundances of both within the bay [38]. Consistent with this richness in foraging habitat, summer-fall visual sighting data and associated modeling of population distributions off the entire western US show highest densities of humpback whales off central California [40–42].

The overall aim of this study is to understand the occurrence of humpback whale song in a biologically rich feeding habitat that supports relatively dense populations of this whale species. Beyond establishing a foundational description of song occurrence in this ecologically important region, this study aims to examine ecosystem variations that may affect humpback whale foraging ecology, and therefore also song occurrence.

Methods

Overview

The study approach begins with detection of humpback whale song from manual analysis of spectrograms, derived from three years of passive acoustic sensing data. Song occurrence is characterized at diel, seasonal, and interannual time scales. At the seasonal and interannual time scales, frequently acquired data from whale watching vessels serve as a proxy for local whale abundance, enabling comparison of song occurrence with statistics of humpback whale sightings in Monterey Bay. The sighting area comprises a portion of the geographic domain for which acoustic modeling defines detectability of humpback whale song at the hydrophone location. At the interannual time scale, ecosystem monitoring data enable characterization of biological factors of importance to humpback whale ecology, which may influence singing behavior. These factors include primary productivity, exposure to a neurotoxic algal compound, and forage species composition and abundance. Hydrographic data and modeling of acoustic transmission loss are applied to consider whether a specific change in geographic patterns of whale habitat occupancy, or changes in acoustic transmission, during an anomalously warm year would have altered song detection at the recording site.

Acoustic recordings

The acoustic recordings for this study are part of a recently established passive acoustic monitoring project [43]. Recordings were acquired through the Monterey Accelerated Research System (MARS) cabled observatory, in the center of MBNMS (Fig 1B). Since 28 July 2015 MARS has supported nearly continuous recording at a sample rate of 256 kHz using an Ocean Sonics icListen HF–an omnidirectional hydrophone with a bandwidth of 10 Hz to 200 kHz. Data stream directly to the Ocean Sonics Lucy software for shore-side recording. In this study we examine the first three full years of recordings, acquired between 1 August 2015 and 31 July 2018. Temporary network outages in the cabled observatory resulted in 94% temporal coverage during the three years.

Acoustic data analyses

The basis for analysis of humpback song was power spectral density (PSD) computation using the Long-Term Spectral Average (LTSA [44]), adapted for efficient execution in an HTCondor pool of computers [45]. LTSA computations were parameterized for routine processing across two frequency ranges (Table 1). Results from the two frequency ranges were merged to produce spectrograms spanning the frequency range of 0 to 6.4 kHz, retaining high resolution (1 Hz) below 1 kHz, the frequency range encompassing most humpback whale song vocalizations. Custom software written in MATLAB was used to aid manual examination of the entire time series of spectrograms, and to register the time periods of song presence for subsequent analyses, as detailed below.

Temporal patterns in song occurrence

The first phase of analysis involved labeling the times of song presence based on manual inspection of spectrograms. Multiple techniques were used to optimize visual detection of song. (1) Spectrogram representations from daily to hourly resolution were simultaneously examined to integrate visual information from longer recording periods covering song sessions with high-resolution detail of song structure. (2) Color scaling for spectrograms at all temporal scales was adaptively constrained to the 1^st and 99^th percentiles of PSD to maximize visual contrast. (3) The frequency axis of the spectrogram was log-scaled for greater detail in the frequency range where most humpback song energy occurs (Fig 1C). Song was distinguished from non-song vocalizations by requiring evidence of phrase and theme structure [2]; songs were almost exclusively detected within sessions (series of songs, e.g. Fig 1C). Song presence was identified only when it was unambiguous, whether from a single singer or multiple singers. Temporal overlap of multiple singers during many periods of song presence precluded consistent evaluation of pause duration between songs. However, song sessions of single singers consistently showed pauses much shorter in duration than the songs, and song was often nearly continuous within sessions. All time within a song session was labeled as song presence, including pauses. If song signal was weak relative to background for a portion of a song session, such that it was not possible to evaluate consistency relative to the rest of the song session, the uncertain portions were not labeled as song presence.

Analyses of song occurrence at all time scales–diel, monthly, and interannual–normalized song occurrence to the corresponding total recording time. Seasonal variation of song presence was examined for all years separately and as a monthly climatology based on all three years. Because song was absent in summer and peaked in winter, the reference year for the interannual analysis extended from 1 August of one year to 31 July of the next. Diel analysis was based on solar elevation: day > 0°, night < -12°, and dusk/dawn between 0 and -12° [17]. Solar elevation was computed for each minute of recording using the NASA SPICE toolkit for MATLAB [46] with the reference location being the ocean surface above the hydrophone. The proportion of minutes of song presence was computed in 12° solar elevation bins, such that one bin (between 0 and -12°) covered dusk and dawn. Because of the large seasonal variation in the histogram of solar elevation and the need to normalize song time to total recording time within each solar elevation category, this analysis would be biased by including months with very low levels of song occurrence, thereby obscuring diel patterns. Therefore, diel analysis is constrained to the months when most song (76%) occurs: November through January. Data from all three years across these months were aggregated for this analysis. An apparently linear relationship between song occurrence (y) and daytime solar elevation (x) was examined using Model 1 linear regression in MATLAB.

Seasonal variation in song length

Using the song presence data extracted in Stage 1, all periods of song were re-examined to quantify the length of individual songs and examine seasonal variation in song length. Consistent with published recommendations [2], the criteria for inclusion in analysis were (1) the ability to distinguish song across its full frequency range, and (2) the ability to define song start and end times consistently from a series of songs in a session, using a marker theme (e.g. Fig 1C). A marker theme for the start of a song was identified from the first song in a session and applied consistently to select start and end times from a subset of clearly defined songs in a session. Song periods containing multiple singers introduced uncertainty in defining individual songs and thus were not included. A total of 1,430 song length measurements were used to quantify variation in a monthly climatology. The entire sample set of song length exhibited a non-normal distribution (Shapiro-Wilk normality test, p < .001). Therefore, to evaluate whether song length changed significantly from month to month, a Wilcoxon rank-sum test was applied. Statistical tests used the Stats package in R, version 3.4.4.

Acoustic modeling informed by humpback song data

Acoustic modeling supports consideration of (1) the source domain for our recordings of humpback song, and (2) the potential influence of variation in habitat occupancy on interannual variation in song detection. To define the frequency at which to model acoustic transmission, a subset of 429 individual humpback songs was selected for analysis. Selection required (1) signal strength sufficient to accurately represent song frequency content, i.e. avoiding faint song signal from distant whales, and (2) absence of extraneous sound sources that would contaminate analysis of song frequency content. For each song a high-resolution (1 second, 2 Hz) LTSA was computed, and for each LTSA of a full song (Fig 2A) the frequency-dependent standard deviation of power spectral density defined a “song profile” (Fig 2B). Comparison of the temporal mean and temporal standard deviation of PSD illustrates why standard deviation (s_PSD) is a better measure of a song profile (Fig 2A and 2B). It provides greater distinction of the frequencies that song units occupy, it is not dominated by variation in the background, and it is less susceptible to tonal anthropogenic noise. The frequency of the maximum in mean s_PSD from all songs defined the frequency at which received levels at the hydrophone would be modeled (Fig 2C).

Acoustic transmission loss at 350 Hz was calculated from a wave-theory parabolic equation model that accounts for absorption in both the water column and the bottom, scattering in the water column and at the surface and bottom, geometric spreading (spherical and cylindrical), refraction, and diffraction [47]. Specification of regional ocean temperature and salinity was based on the September climatology from the US Navy Generalized Digital Environmental Model (GDEM). Bathymetry was specified at 250 m resolution. The source level of 175 dB re 1 μPa was based on published in situ measurements [48] and used to compute received levels at MARS, to characterize the spatial domain around the hydrophone over which song would be detectable under different noise conditions. Sound source depths were specified as 10, 20, or 30 m. Because the results were very similar for these depths, mapped results from the shallowest sound source (10 m) are presented to include results for sound sources in shallow water close to the coast (inner ~1–5 km), which can be important habitat for humpback whales in Monterey Bay.

Visual sighting data

To interpret temporal variations in humpback song detection in relation to visual sightings, this study used sighting data from commercial whale watch operations (Monterey Bay Whale Watch). The primary advantages of these data are frequent sampling, with visual surveys averaging 326 ± 21 days per year during the period of this study, and consistent coverage of an area near the hydrophone (southern Monterey Bay within ~15 to 35 km from the hydrophone, Fig 1B). These were not systematic surveys, but rather represent number of whales sighted on targeted trips, with destinations often being determined by where whales were known to have been sighted already. However, we expect that lower frequency (seasonal, interannual) patterns in regional presence will be reflected in statistics computed from these frequent observations. Although abundance of adult male humpback whales would be most relevant to song, the sighting data cannot support this degree of specificity. Since the whale-watching record variably included 1–3 vessels making 1–3 partial or full-day trips on each day, daily whale counts were normalized to a constant unit of effort, measured as humpback whales sighted per half-day trip. For comparison with humpback song occurrence at interannual time scales, temporal averaging was constrained to the annual period of song detection, September of one year through May of the next year.

Ecosystem data

Physical data

Variations in sea level reflect changes in the depth of nutrient source waters, while variations in wind-forced upwelling define the primary transport mechanism of nutrients from depth into the euphotic zone. These two physical variables are key predictors of the primary productivity supporting forage species assemblages that baleen whales ultimately rely upon [32,36,37]. Sea level anomaly (SLA) data were from the Integrated Multi-Mission Ocean Altimeter Data for Climate Research sea level data set, spanning 1993 through 2018 [49]. SLA data were subset for the CCS according to analytically defined boundaries [50] (23°N to 43°N within 800 km of the coast) and averaged spatially and seasonally, then the seasonal climatology was removed to emphasize interannual variation. Wind-forced upwelling was examined using the coastal upwelling index [51] for 36°N, 122°W (near Monterey Bay) obtained through the NOAA ERDDAP server [52].

While SLA data do not extend onto the continental shelf, ocean temperature data do and can thus describe physical variations within humpback whale foraging habitat near the coast, as well as pycnocline depth and strength that can affect acoustic transmission from shallow sound sources to a deep receiver. Regional ocean temperatures were examined using remotely-sensed and in situ data. Satellite data were Version 4.1 Multi-scale Ultra-high Resolution (MUR) sea surface temperature (SST) monthly anomaly, 2002 through 2018, produced by NASA Jet Propulsion Laboratory and accessed through the NOAA ERDDAP server [52]. In situ data were from three sources. The first was a glider section time series along California Cooperative Oceanic Fisheries Investigations (CalCOFI) Line 67, which extends offshore of Monterey Bay. This glider time series, 2007 through 2018, is part of a larger glider research program in the CCS [53]. Both satellite and glider time series were processed to seasonal anomalies. Pycnocline variation was examined using continuous observations from a CTD string on mooring M1, at the mouth of Monterey Bay (middle white diamond in Fig 1B) and periodic CTD cast data from station C1 (nearshore white diamond in Fig 1B). These locations are part of a time series observation program [54]. Stations in the bay are most relevant because our focal question is whether preferential habitat occupancy within Monterey Bay during a warm anomaly would affect acoustic transmission to the MARS hydrophone.

Biological data

Biological data included primary productivity, concentrations of the neurotoxic algal compound domoic acid (DA), and forage species abundances. Primary productivity measurements were from a time series program [54] regularly sampling three stations within and outside of Monterey Bay (white diamonds in Fig 1B). From the results of carbon uptake incubations, water column primary productivity was computed by integrating across the depth range spanning all light levels of incubation: 100, 50, 30, 15, 5, 1, and 0.1%. The CCS is prone to harmful algal blooms (HABs) that result from food web transfer of the biotoxin domoic acid (DA) [55–57]. DA is produced by several species of the diatom genus Pseudo‐nitzschia and can cause illness, heart and brain damage, reproductive failure, and mortality in marine mammals, including baleen whales [55,58–61]. For this study we use weekly measurements of particulate domoic acid (pDA) from Santa Cruz Wharf (location in Fig 1B). Quantification of pDA was by liquid chromatography‐mass spectrometry [62]; data were obtained through the Southern California Coastal Ocean Observing System data server [63].

To examine forage species consumed by humpback whales, standardized abundance data were derived from the NOAA-NMFS Rockfish Recruitment and Ecosystem Assessment Survey (RREAS). The RREAS conducts an annual mid-water trawl survey during late April through mid-June to assess ocean conditions and the abundance and distribution of epipelagic micronekton off California (1983-present). The RREAS samples a variety of epipelagic micronekton utilized by mid and upper trophic level predators, including krill, pelagic juvenile rockfishes (Sebastes spp.) and northern anchovy (Engraulis mordax) [64]. Micronekton samples were collected at fixed sampling locations (Fig 1B) during night using a modified Cobb midwater trawl with a 9.5-mm cod-end liner; 15-minute tows were made at each station with a headrope depth of 30 m [64]. After each haul, all taxa were enumerated, partitioned into either young-of-the-year (YOY) or adult classes as necessary, and relative species abundance was measured as catch-per-unit-effort (CPUE) per station. For a description of the spatial distribution and temporal variability of numerically dominant taxa, community structure, and an assessment the variability of biodiversity, see [38], [65], and [66], respectively. Due to their importance in the diet of humpback whales [67], we quantified the variability of abundance for total krill and northern anchovy during 2015–2017 using data from all net-tow stations [68,69] (locations in Fig 1B). Anomalies of abundance were calculated from the spatial mean CPUE across all stations sampled between Bodega Bay and Pt. Sur (Fig 1B) and were standardized according to the long-term mean and standard deviation over 1990–2017 [64,68]. Forage species data are available from NOAA [70].

In examining ecosystem variations in relation to interannual changes in song occurrence, leading temporal lags are appropriate because ecosystem conditions preceding the annual (fall) rise of singing behavior could influence whale population dynamics, foraging ecology, and behavior. Beginning with forcing of primary productivity by wind-driven upwelling, we quantify cumulative upwelling between the annual onset in spring through the end of the calendar year. For the three winters examined, strong downwelling caused by winter storms consistently ended by March. Therefore, a consistent measure of cumulative upwelling spanned March through December of each year. This same annual period was used for mean sea level in the CCS, spatially and temporally integrated primary productivity in the Monterey Bay region, and particulate domoic acid (pDA) concentrations. For pDA we quantify the annual temporal span of detection, maximum concentrations, and a potential exposure metric integrated annually by summing the product of each concentration measurement and the time period it represents (1 week). Forage data are from spring/summer of each year and provide an indicator of the overall condition of forage species populations within the central CCS [71]. To describe ecosystem physical patterns from satellite and in situ temperature data, the averaging period was constrained to the second half of the year because glider sampling was missing for the second quarter of one year, and the consistently sampled period of the third and fourth quarters spans late-summer foraging and the annual period when song begins and rises to its annual peak.

Results and discussion

Humpback song detection

The region around the hydrophone from which humpback song would be detectable under different levels of background noise is characterized from received level for a sound source having a frequency (350 Hz; Fig 2) and source level (175 dB re 1 μPa, peak-to-peak) representative of humpback song [48]. Model results for 10, 20, and 30 m deep sound sources showed very similar patterns. We show the results from 10 m (Fig 3A) because the 10 m isobath reaches closest to the coast, and it is important to evaluate reception of humpback song signal from throughout the bay, including nearshore habitat (Fig 3B). These results show that under relatively loud background noise conditions (110 dB re 1 μPa) humpback song would be detectable at MARS only if the singing whale was within ~ 20 km of the hydrophone (Fig 3A). Under relatively quiet background noise conditions (70 dB re 1 μPa) humpback song would be detectable at MARS if the whale was within a much larger area, including most of Monterey Bay (Fig 3A). Our humpback whale sighting data are from southern Monterey Bay, which is within the domain of MARS acoustic detection of humpback song under background noise conditions below ~ 97 dB re 1 μPa (Fig 3A). Thus, comparison of sighting and song data is appropriate.

Time series overview

The three-year record of humpback song detection shows large seasonal and interannual variations (Fig 4, black bars). Song began during fall of each year (September–October) and rose steeply to a peak in November or December. The temporal coverage of song reached a maximum of 86% during November 2017. Song occurrence dropped steeply after January of each year, exhibiting highly variable levels during February through May. Song was not detected during June through August in any year. Integrated over the annual cycle (August through July), total song occurrence increased by 44% between the first and second years, and by 55% between the second and third years. Considering the recording coverage throughout the time series (Fig 4, gray bars), song occurrence patterns described at monthly resolution for the full time series are robust. Since seasonal variation informs analyses of diel and interannual variation, we begin with describing seasonal aspects of song occurrence, then diel and interannual variability. At each temporal scale, we evaluate hypotheses for observed patterns of song occurrence.

Seasonal

Monthly mean song levels show that 76% of all song occurs within one quarter of the year, November through January (Fig 5A). Between the start of song in fall through the end of the annual peak in January, song length increases significantly each month (p < 0.01, Wilcoxon rank sum). This increase is evident in minimum, maximum, and quartiles of song length (Fig 5B). Median song length increases from 4.9 minutes in October to 10.1 minutes in January. The annual peak in song coincides with the winter trough in day length (Fig 5C), and with sighting-based evidence of population movement (Fig 5D). Following high levels of humpback sightings during summer and through the start of the song peak in November, sighting levels drop steeply during December and January (Fig 5A and 5D). This seasonal pattern in sightings during our study is consistent with that based on the full sighting record, 2003 through 2018 (Fig 5D). Our sighting data cannot define how much of the humpback population remains in the region, and within the detection range of the hydrophone, after leaving Monterey Bay. However, these indications of population movement during November–January (Fig 5D), coincident with the annual song peak (Fig 5A) and day length trough (Fig 5C), are consistent with seasonal patterns of migratory and vocal behavior driven by changing ecosystem conditions leading into the breeding period.

The South Pacific analogue to the CCS is the Humboldt Current System (HCS). Like the CCS, physical distinction of the HCS caused by coastal upwelling is evident in the cooling gradient between the subtropical gyre and the eastern margin (Fig 1A), and it is important habitat for humpback whales. In the HCS the first continuous acoustic monitoring off southern Chile also demonstrated humpback song occurrence in feeding habitat [17]. Song was detected throughout the 130-day study during February–June 2012. Although this study did not sample a full year, it revealed increasing song occurrence during austral fall and a maximum in winter (June), consistent with our results. Studies of humpback song in the western and eastern North Atlantic also exhibit seasonal patterns consistent with those observed in the northeast Pacific during our study. A decade of acoustic monitoring west and north of the British Isles, eastern North Atlantic, detected humpback song during October–April and equatorward progression during January–March, consistent with migration to low-latitude breeding habitat [13]. Research in Stellwagen Bank National Marine Sanctuary, western North Atlantic, showed increasing song occurrence during the fall–winter transition, when males prepare to migrate to low-latitude breeding grounds [15]. Seasonal variation in hormone levels and associated neurological changes have been hypothesized as a regulatory factor in the seasonal variations in humpback song occurrence [14,15,72]. Measurements of testosterone levels in blubber samples from North Pacific adult male humpback whales show minimum levels in summer, maximum levels in winter, and an increasing trend during fall [73], supporting the hypothesis of hormone levels as a regulatory factor in singing behavior.

Seasonal increase in the length of song sessions has been reported from analysis of recordings from the northwest Atlantic [15]. To our knowledge, seasonal variation in the length of individual songs has not been examined in previous studies. The steady seasonal increase in song length we observed motivates consideration of possible causality. Primary hypotheses include learning, which is the basis for song evolution [8–10], and behavioral influence of seasonally rising hormone levels [73]. Singing behavior during periods of heightened song activity within a population may lead to increasing song length through cognitive development. Conditions that covary with population density, such as prevalence of song activity and diversity of song unit types, may also play a role in learning and associated changes in song attributes (complexity, length). Supporting relatively high population densities, favorable feeding habitats such as Monterey Bay may enhance density-dependent influences on learning. Seasonally increasing hormone levels may affect singing behavior and song attributes, including length. Beyond the scope of this study, how song complexity and composition change with increasing song length is of interest in understanding humpback whale behavior.

Diel

As described in the methods, it is appropriate to constrain examination of diel variation to the time of year when song is prevalent. During the annual peak in song occurrence, November through January (Fig 5A), diel variation is pronounced (Fig 6). Maximum song occurrence is during the night, averaging 69% temporal coverage. Song occurrence decreases during dawn and dusk (52%), and further decreases with increasing solar elevation during the day to reach a minimum near solar noon (30%). Decreasing song occurrence with increasing solar elevation (Fig 6) is strongly linear (Model 1 linear regression, r² = 0.99). Our results on diel variation are consistent with those from the study in the HCS [17]. In that study, when song detection was elevated (April—June), the highest levels were during dark and lowest levels were during daylight. Observations in breeding habitat around Hawaii have shown diel variation in song occurrence, maximum at night and minimum during the day, as well as increasing male competitive interactions during the day [74,75]. Alternative hypotheses for diel variation off Hawaii included whales singing louder or moving closer to the nearshore recorder during the night, or more whales singing at night [74]. Assuming more whales singing at night, the authors note that vision in daylight plays a key role in competitive group formation, and that males may switch to vocal advertisement as the primary mating strategy in darkness. These hypotheses proposed to explain diel variation in song occurrence within breeding habitat are relevant to its occurrence in the feeding and migratory habitat of our study. However, the importance of foraging activity to humpback whales in the Monterey Bay region introduces an additional consideration. If foraging strategies rely on vision in daylight, shifting behavior to foraging during daylight may reduce singing. It is not known whether humpbacks in this region sing in close spatiotemporal proximity to feeding behavior, as observed off Antarctica [16].

Interannual

Interannual variation in humpback song occurrence was marked by 44% and 55% increases between consecutive years (Fig 4). Between the first and second years, the increase was associated with a longer temporal window of song detection between fall and spring, increasing from 197 to 233 days, and with greater occurrence during the peak months of November through January (Fig 4, Table 2). Between the second and third years, the increase was associated with yet greater occurrence during the peak months as well as a broader peak that started ~1 month earlier and ended ~1 month later (Fig 4, Table 2). Quantified on an annual basis, song temporal coverage increased from 11% in the first year to 25% in the third year, more than doubling during this period (Table 2). Alternative hypotheses for the interannual variation in song detection include interannual differences in the (a) local abundance of singing whales, (b) geographic distribution of singing whales in relation to the hydrophone location–in a way that would affect detection of song, (c) sound transmission to the hydrophone–as influenced by changes in hydrography, and (d) levels of singing activity, as may be influenced by animal health and time allocation to singing behavior. We consider these hypotheses using observational data and acoustic model results.

Environment

Consideration of all hypotheses for the strong interannual trend in song occurrence requires the context of interannual variations in ecosystem conditions. The first year of recording coincided with the highest sea level anomalies (SLA) in the CCS in the last 25 years (Y1 in Fig 7A). While SLA remained positive during the second and third years of recording, the magnitude of the anomaly was greatly reduced after the first year (Y2 and Y3 in Fig 7A). The exceptionally high SLA was due to a strong and persistent warm anomaly in the northeast Pacific [76]. Strongly positive SST anomalies prevailed throughout the CCS during 2015, and they ranged between weakly positive to weakly negative during 2016–2017 (Fig 7B). Subsurface temperature data extended between the coastal upwelling zone and the inner CCS off Monterey Bay (vertical sections in Fig 7B). These data show that the warm anomaly affected much of the upper water column across the oceanic margin, and that it increased toward the coast and toward the surface (Fig 7B, 2015). During 2016–2017 the near-coastal warm anomaly was absent (Fig 7B, east of ~123°W) although subsurface warm anomalies persisted further offshore. The major physical changes within ~ 80 km of the coast extended across primary foraging habitat of humpback whales.

Changes in local abundance of whales?

The first hypothesis for causality of the positive interannual trend in song is a positive trend in the local abundance of singing whales. Our sighting data from southern Monterey Bay are within the domain of MARS acoustic detection of humpback song (Fig 3A). Although our sighting data cannot constrain consideration to possible singing (adult male) whales, they represent extensive sampling. For the three years of recordings the sighting data span 82%, 82%, and 95%, respectively, of all days during months when song was detected—September through May (Fig 5A). Unlike the strong positive trend in song occurrence, sighting rates do not show a positive trend (Fig 8A). The increase in sightings between the second and third years may partly explain the increase in song detection between those years, however variation in the local abundance of whales is not indicated as causal in the song increase between the first and second years and cannot explain the steady positive trend over three years (Fig 8A). Further, coincidence of the lowest levels of song occurrence with the highest levels of sightings during the first year (Fig 8A) emphasizes the potential for behavioral regulation on the interannual time scale, as indicated at the seasonal time scale (Fig 5).

Changes in the geographic distribution of whales?

The second hypothesis is that interannual variation in song detection was influenced by changes in the geographic distribution of singing whales in relation to anomalous habitat conditions (Fig 7). A previous study suggests that rorqual whales, dominated by humpbacks, preferentially occupied habitat within Monterey Bay during a warm anomaly caused by El Niño and associated with reduced krill abundance [77]. Effects of the warm anomaly during our first year of recording share some similarity with El Niño effects, and a key question is whether preferential occupancy within Monterey Bay would reduce detection of song at the deep MARS site. Considering a 60° arc extending 30 km onshore from the hydrophone into Monterey Bay, mean received levels decrease 21 dB re 1 μPa between the hydrophone and inner shelf waters (Fig 8B). For a similar domain mirrored in the offshore direction, mean received levels decrease 20 dB re 1 μPa (Fig 8B). Thus, attenuation of song signal between singing whales and the hydrophone would not have been significantly greater for whales in Monterey Bay compared to offshore. These results do not support the hypothesis that preferential occupancy within Monterey Bay, if it occurred during the severe warm anomaly in first year of recording, would have reduced detection of song at the MARS site and thereby contributed to the positive interannual trend.

Changes in the acoustic transmission of whale song?

The third hypothesis is that changes in hydrography during this highly anomalous period (Fig 7) may have affected acoustic transmission to the hydrophone. If pycnocline depth was anomalously deep during the warm anomaly, relatively more song may have originated above the pycnocline and encountered greater attenuation due to the pycnocline. This is particularly relevant to song originating inside the bay, which would have a shallow angle of incidence relative to the hydrophone. Also, if pycnocline density gradients were anomalously strong during the warm anomaly, transmission loss may have increased, in turn reducing song detection at MARS. The potential role of variation in the pycnocline was evaluated using density profiles from stations M1 and C1 (white diamonds in Monterey Bay in Fig 1B). Although water density was anomalously low due to presence of the warm anomaly during the first year, the magnitude and vertical distribution of pycnocline density gradients was not anomalous. At both locations the maximum average vertical density gradient was shallower than 30 m during the song periods of all three years. Ship CTD profiles at the inner bay station were infrequent and inconsistent, with 8, 4, and 3 profiles during the song periods of each year, respectively. In contrast, CTD data from M1, at the mouth of Monterey Bay (middle white diamond in Fig 1B), was continuously recorded throughout all three song periods and is a better resource for comparison of magnitude. Across the three annual song periods, the magnitudes of the maximum average vertical density gradient were similar: 0.016, 0.016, and 0.018 kg m^-4, respectively. The similarity in the depth and magnitude of pycnocline gradients across years indicates that hydrographic variation was unlikely to have influenced interannual variation in song detection at MARS.

Changes in singing behavior in relation to ecosystem variations?

Extensive observations show major ecosystem variations relevant to humpback whale ecology. Further, these variations exhibit trends consistent with the observed trend in song occurrence. The trend of integrated primary productivity (IPP) in the Monterey Bay region closely parallels that of song occurrence (Fig 9A and 9B). IPP increased 60% between the first and third years (Fig 9B). Through trophic linkages in the short food chains between phytoplankton and baleen whales, such large changes in primary productivity could influence food resources for humpback whales. The trend in primary productivity followed that of cumulative upwelling (Fig 9C), i.e. wind-forced nutrient supply to the euphotic zone. It is also consistent with sea level variation (Fig 7A). Elevated sea level corresponds to a deeper thermocline (nutricline) and would cause relatively low levels of nutrient recruitment from wind-driven upwelling [78]. The decreasing trend in sea level thus represents shoaling of the nutricline and a positive trend in potential upwelling nutrient supply.

Consistent with the increasing trend in primary productivity, krill transitioned from a strong negative anomaly in the first year, to near normal in the second year, to a strong positive anomaly in the third year (Fig 9D). Adult anchovies exhibited depressed abundances during all three years (Fig 9E). The warm anomaly of the first year coincided with a large positive anomaly in young-of-year (YOY) anchovy (Fig 9E). One way these biological anomalies may have influenced humpback whale ecology is through prey switching. Large scale climatic variations in the CCS cause changes in the abundance and diversity of forage species [66,68], which in turn leads to prey switching in many species, including humpback whales [69]. If the forage species anomalies increased temporal and energetic demands of foraging, time and caloric energy available for singing behavior may have decreased.

Beyond the magnitude of primary productivity, a potentially influential factor is the extent to which primary productivity generates algal toxins that can affect marine mammal health and behavior. Our first year of recording coincided with a severe toxic algal bloom. Detection of the toxin, domoic acid (DA), in marine mammals—including whales, dolphins, porpoises, seals, and sea lions—spanned the largest geographic extent ever recorded [61]. Concentrations of DA in the particulate fraction (pDA), essential to trophic transfer of toxin, reached the highest levels ever recorded in Monterey Bay [79]. At a fixed location in northern Monterey Bay from which weekly monitoring is conducted (Fig 1B), integrated pDA was highest leading into and during this first year of recording, and it exhibited a decreasing trend thereafter (Fig 9F). Between the first and third years, integrated pDA decreased by 82%, and similar trends existed in the temporal span of pDA detection and maximum pDA concentrations (Fig 9F, Table 3). This neurotoxin can cause disease and mortality in marine mammals, including baleen whales [55,58–61]. Coincidence of highest potential DA exposure (Fig 9F) with lowest song activity amid the highest levels of humpback sightings (Fig 8A) suggests the possibility that high levels of DA exposure may suppress singing behavior.

Conclusions

This study extends the growing body of observations showing occurrence of humpback song in feeding and migratory habitats. In the northeast Pacific, the Monterey Bay region is important feeding and migratory habitat for humpback whales. The abundance of humpback whales off the U.S. West Coast has increased from about 900 individuals in the early–mid 1990s to about 2,900 whales in 2014, as this species has recovered from historic whaling impacts [80]. Sighting data indicate increased humpback whale abundances in Monterey Bay since 2013 (Fig 10), further supporting the critical role of this biologically important area [40] to humpback whale populations of the northeast Pacific.

Levels of song occurrence in this feeding and migratory habitat are high, occupying 9 months of the year with temporal coverage reaching as high as 86% of the time during the annual winter peak. Absence of song during summer, when visual sighting levels are high, emphasizes behavioral regulation of song activity. Coincidence of the annual trough in day length with both the annual song peak and evidence of population movement out of Monterey Bay is consistent with the underlying seasonal ecosystem variations and associated migratory patterns required for successful foraging and reproduction. These seasonal patterns observed in the northeastern Pacific are similar to those observed for humpback whales in the southeastern Pacific, northeastern Atlantic, and northwestern Atlantic. The rise of song activity leading into the breeding season is consistent with a central role in reproductive ecology, yet the prevalence of song in feeding and migratory habitats motivates further consideration of the greater spectrum of possible roles in social interaction, cultural transmission, and perception. Favorable foraging habitats that support relatively dense populations of whales may enhance density-dependent influences on song development and learning.

Interannual variation in song was great, exhibiting a positive trend over three years during which the temporal coverage of song doubled. Espousing the method of multiple working hypotheses [81], and within the limits of model and observational data, we considered alternative hypotheses. The hypothesis having some support is that exceptionally great physical and biological variations in the ecosystem influenced singing behavior, perhaps through changes in foraging ecology. While examination of causality between complex ecosystem variation and singing behavior can only be speculative, possible influences of the observed biological anomalies on humpback whale singing behavior are consistent with the observed changes in song detection. During the first year of recording, song levels were the lowest despite having the highest rate of humpback sightings in an area within the acoustic detection range of the hydrophone. This coincided with the greatest physical anomalies in the CCS in the last quarter century, severe depletion of a primary food resource–krill, and extremely high levels of a neurotoxin known to harm whale health–domoic acid. Relaxation of these extreme anomalies over the following two years, likewise, may have resulted in increased song activity by supporting greater health and reducing the time and energy requirements of foraging.

NOAA's National Marine Sanctuaries and National Marine Fisheries Service play a central role in understanding and protecting species and their habitats, with increasing interest in natural and anthropogenic sounds [82]. Passive acoustic monitoring is an essential and effective method of studying variations in marine mammal presence and behavior. The new knowledge of humpback whale song activity from this study has management implications, for example informing decisions on when to allow anthropogenic sounds through permitting processes. Interdisciplinary observations of complex marine ecosystems are essential to understanding ecosystem health and its long-term changes [83]. The integrated sound and ecosystem analysis presented in this study provides an example approach to comparative studies across regions having different mammal species, conservation threats, and opportunities to augment stewardship of protected species and habitats.