Some things never change: multi-decadal stability in humpback whale calling repertoire on Southeast Alaskan foraging grounds

Investigating long term trends in acoustic communication is essential for understanding the role of sound in social species. Humpback whales are an acoustically plastic species known for producing rapidly-evolving song and a suite of non-song vocalizations (“calls”) containing some call types that exhibit short-term stability. By comparing the earliest known acoustic recordings of humpback whales in Southeast Alaska (from the 1970’s) with recordings collected in the 1990’s, 2000’s, and 2010’s, we investigated the long-term repertoire stability of calls on Southeast Alaskan foraging grounds. Of the sixteen previously described humpback whale call types produced in Southeast Alaska, twelve were detected in both 1976 and 2012, indicating stability over a 36-year time period; eight call types were present in all four decades and every call type was present in at least three decades. We conclude that the conservation of call types at this temporal scale is indicative of multi-generational persistence and confirms that acoustic communication in humpback whales is comprised of some highly stable call elements in strong contrast to ever-changing song.


Results
Within 114.9 hours of recordings, we identified a total of 914 high quality calls that fit our inclusion criteria. Recordings spanned 140 unique days within five separate years across four decades (Table 1). A total of 175 individuals were photographically identified in conjunction with this study in 1997, 2007 and 2008. The minimum number of whales present during recording periods ranged from 6 to 90 (Table 1). Animals were not localized in this study; it is unknown whether all whales within the area were contributing to the recorded calls, or which specific individuals were calling.
Classification. Sixteen known call types nested within for vocal classes were described in Southeast Alaska in 2012 (Table 2). A comprehensive description of each class and vocal type is available elsewhere 36 . Aural-visual analysis (AV) assigned 367 calls to the Low Frequency Harmonic (LFH) vocal class, 303 to the Pulsed (P) vocal class, 79 to the Noisy-Complex (NC) vocal class, and 165 to the Tonal (T) vocal class ( Table 2).
Rotated principal component analysis (PCA) output indicated that the use of two principal components was adequate to encompass the variability of the data (χ 2 = 1409.98, p < 0.000001). The first rotated component (PC1) corresponded most closely to frequency parameters (lower, median, peak, and start frequency values; see Table 3 for variable descriptions), indicating that as PC1 increases calls generally increase in frequency. The second SCIENtIFIC REPORTS | (2018) 8:13186 | DOI: 10.1038/s41598-018-31527-x rotated component (PC2) corresponded most closely to temporal parameters (duration (negative relationship), amplitude modulation rate, frequency modulation rate, (Table 3)), indicating that as PC2 increases calls generally grow shorter and change more quickly.   A Classification and Regression Tree (CART) assigned 85% (root-node error, n = 754) of vocalizations to the same call type as AV classification (Table 4). Bout, duration, frequency ratio, PC1 and amplitude modulation were all important splitting variables. A random forest analysis correctly classified most of the calls (out-of-bag error rate (OOB) = 27%). Consistent with the CART analysis, the variables most important for splitting decisions were amplitude modulation, duration, bout, median frequency, PC1 and frequency ratio. Misclassifications were common among call types with low sample sizes 36 . These calls were re-reviewed manually by a second observer; upon observer agreement calls were categorized according to AV classification (Fig. 1).
Temporal Stability. Twelve of the sixteen call types were found in both 1976 and 2012 (Table 2, Fig. 1,   Fig. 2), and all call types were detected in at least three decades. Eight call types were represented in all four decades. Initial data analysis indicated that there were no significant differences in acoustic parameters between 2007 and 2008 (χ 2 > 4, p > 0.1) and those years were pooled to ensure a large enough sample size for statistical inference across call types and between decades. The acoustic parameters of calls in this study were variable, though within the margin of variability reported for each call type (Table 5) 36 , and in some cases fine-scale acoustic parameters varied significantly between decades (   Table 4. Confusion matrix of classification and regression tree (CART) output (top) compared to Aural-Visual (AV) call type assignment (left). Overall classification agreement was equal to 85%. For all call types except for droplets, mean PC1 values were significantly higher in calls produced in 1997 ( Fig. 3) than in either the 2007-08 or 2012 data sets. Mean PC1 values for droplets were highest in 1997, but mean PC1 values for droplets in 1997 were not significantly different than calls produced in the 2000's. In general, calls produced in 1997 were higher in frequency than in the other two decades. PC1 values tended to be lowest for calls produced in 2007-08 and intermediate in 2012 (Fig. 3). There was no obvious linear change in PC1 values over time. Where significant differences occurred, PC2 values showed a general downward trend between 1997 and 2012 (Table 5); PC2 estimates indicate that across call types calls produced in 1997 were shorter and had higher modulation rates, and calls from 2012 were longer and had slower modulation rates (Fig. 3).

Discussion
The results of this study indicate that some humpback whale non-song call types in Southeast Alaska are conserved in the acoustic repertoire at the decadal timescale. Twelve of the 16 call types were present in both the 1976 and 2012, demonstrating the persistence of these calls over 36 years. All other call types persisted over multiple decades between the 1997 and 2012. The absence of identical call types in recordings across the four decades does not necessarily indicate that certain call types were absent from the repertoire, since it is unlikely that the non-song calls aggregated for this study comprise the entire vocal repertoire of the entire population at any given time. Limited sampling effort (Table 1) might favor the recording of some call types over others; nonetheless, within the scope of this dataset the overall pattern of long-term conservation of call types within the acoustic repertoire is clear.
Demographic data are consistent with the assertion that stability at this temporal scale indicates call types are conserved across several generations. Given that females mature by the age of 13 and typically reproduce every 1-3 years 55-57 , new whales in our survey regions were born, reached sexual maturity, and gave birth to offspring that subsequently grew to sexual maturity and also gave birth over the 36-year duration of this study 58 . This is reflected in the dramatic increase in population size from the 1970's to the 2010's, which was estimated in 2008 to be more than five times as large as in 1979 52,53 . Demographic studies show that population growth in Southeast Alaska is primarily due to long-term maternally directed site fidelity and birth 54,55,59,60 with little evidence of immigration into the feeding areas of the population 59 . In further support of this point, in Glacier Bay National Park, where humpback whale monitoring efforts have been underway since the 1970's, almost half of the humpback whales first identified as calves returned to their maternal foraging grounds 55 . Thus, it is quite likely that throughout the study period multiple generations overlapped spatially and temporally within the survey region and were recorded by our hydrophones.
In this study, some acoustic parameters varied significantly, though non-linearly, across decades. Most vertebrate sounds demonstrate within-call variation, related to individual anatomy, behavioral or environmental context (see Fig. 2 for a visual example from this study), which does not contradict placement into call classes or  19,61 . In the eastern Australian humpback whale population call types classified with high degrees of confidence also exhibit fine-scale acoustic variability over time 37 . In the same population, acoustic parameters of calls also varied with social context, though call type assignment was robust 48 . Similarly, Deecke et al. 61 demonstrated changes in acoustic parameters of killer whale calls, although the overarching call structure remained stable 61 . Thus, although there were observed differences in acoustic parameters between decades, the persistence of the call types themselves in the acoustic repertoire should be considered stable.
There were some prominent trends in fine-scale acoustic behavior between years. Calls from 1997 were consistently higher and shorter than other time periods and calls of almost all types grew longer from 1997 to 2012. We offer a few possible explanations. Individuals may have adjusted vocal parameters in response to ambient noise conditions; to avoid acoustic masking individuals must either increase their calling amplitude or spectrally or temporally shift their vocalizations 62 . Alternatively, and equally plausibly, changes in acoustic parameters may reflect differences in social context rather than change over time. According to motivational-structural rules, mammalian vocal sounds encode information about a sender's motivational state, enabling a receiver to assess the likelihood of certain behaviors occurring 63 . It has been suggested that motivational state is encoded in humpback whale calls on migratory corridors 46 . In eastern Australia humpback whales used call types that were higher in frequency at times associated with increased arousal levels (i.e. affiliating groups for mid-level arousal, groups of competing males for high arousal) 46,48 . In the present study, recordings made in 1997 were typically made in association with large groups of whales engaged in coordinated foraging events. In contrast, despite dedicated observer effort coordinated foraging was not observed in 2007 or 2008 64,65 , and whales in 2012 were typically foraging alone, in low densities, and in some cases were vocalizing in isolation 41,48 . Coordinated foraging is likely to increase arousal and necessitates fine-scale cooperative interactions between individuals 66 . This indicates a difference in audience and plausibly a difference in motivational state between datasets in the different decades.  Table 5. Means (bold) and standard deviations of selected acoustic parameters for the subset of call types found across all four decades. PC1 corresponds to frequency components. PC2 corresponds to temporal components". Chi-squared statistic and P-values refer to results from Kruskal-Wallis test for equality of medians. See Table 6 for detailed results of Dunn's test for multiple comparisons. Significant differences are italicized and underlined. High-frequency, short-duration vocalizations that are frequency-modulated are thought to have an 'appeasing effect' on receivers 63 and have been associated with groups of affiliating humpback whales 46 . Consistent with this theory, calls from 1997, where groups of animals were engaged in complex coordinated activities, were generally higher, shorter, and more frequency modulated than in decades where social interaction was more limited during recording periods. Our result, that calls grew generally longer over the duration of this study, corroborates findings in other mysticetes, where call duration increased over time, possibly in response to elevated ambient noise 33,58,67 . However, statistically significant temporal differences in call features may or may not have practical significance. For example, PC2 values, which related to temporal features, were lower for droplet and teepee calls produced in the 2012 than in the 2000's, indicating that calls in 2012 were longer in duration; the differences in mean duration, however, were only a fraction of a second. Because these calls are pulsed and often occur in short bouts 36,68 , a difference in a fraction of a second may reflect a longer duration between call units, or could be an artifact of temporal smearing with distance from the hydrophone. While every effort was made to account for differences in environmental conditions and recording units, the temporal parameters of a call can still be affected by attenuation, reverberation, and other propagation effects. Should these very fine-scale measurements be indicative of true variation in temporal characteristics, they still may not represent enough biological variability to be detectable, or meaningful to a receiver. Additional research into auditory discrimination in humpback whales would be extremely valuable.
In the absence of calibrated ambient sound recordings, demographic information, and fine-scale behavioral sampling, it is not possible to quantify whether motivation, ambient noise conditions, or individual variability contributed most to changes in calling behavior. Future investigation into acoustic variability as a function of social context and ambient noise on feeding grounds would be useful for testing hypotheses about acoustic communication in foraging humpback whales.
The selective pressures that maintain and shape characteristics of communication signals are dictated by social structure 13,69 . Species like killer whales that live in discrete family units with little or no natal dispersal exhibit temporally stable and stereotyped pod-specific vocalizations [70][71][72] . Bottlenose dolphin societies characterized as 'fission-fusion' , with social affiliations ranging from ephemeral to long term, employ stable signature whistles to convey individual identity 73,74 . We have demonstrated that humpback whales also produce call types that persist across decadal time scales. Given the evidence of long-term affiliation between humpback whales at high latitudes 57,75,76 we suggest that some call types may be used communicate identity over time and space. Acoustic identity cues are common across taxa (e.g. northern fur seal (Callorhinus ursinus) mother-pup recognition calls 77 , king penguin (Aptenodytes patagonicus) contact calls 78 , Mexican free-tailed bats (Tadarida brasiliensis mexicana) isolation calls 79 ; see Tibbetts & Dale, 2007 for a review 80 ), and -as evidenced by savanna elephants -can be valuable among far ranging social animals 15,16 . Some of the most commonly documented call types identified in  Table 6. P-values from Dunn's test pairwise comparisons for differences in median rotated principal components between decades. PC1 corresponds to frequency components. PC2 corresponds to temporal components. Significant differences between decades are italicized and underlined.
this study, whups and growls, have been frequently identified elsewhere 38,39,49 , and may act as contact calls 34,45 . These call types are acoustically variable, and highly persistent over time on both migratory corridors and foraging grounds; this variability may reflect signature information used to confer identity over time. We propose that these call types are good starting points for investigating individual variation and recognition in this species. Calls may also persist in the vocal repertoire in Southeast Alaska because they are functionally specific to foraging activities in this region. Feeding calls are behaviorally linked to foraging on Pacific herring (Clupea pallasii) 41,51,60 and are closely matched to hearing capabilities of this prey, which are most sensitive in the 200-500 Hz range 81 . Under experimental conditions, playbacks of feeding calls have elicited a "flee and clump" response in Pacific herring, which presumably increases the whales' foraging efficiency 82 . Feeding calls are most commonly documented among groups of whales, where they may coordinate individuals; however, they are also produced by solitary animals, and may serve a prey manipulation function in this context 41 . This would explain both the call stereotypy -feeding calls were correctly classified 98% of the time in this study (Fig. 2, Table 4) -and also differences in call parameters between decades, which are likely related to social context (e.g. solitary foragers versus group foragers) or prey behavior (e.g. school size, school location). Although this call type does not exclusively serve a social function, it is likely to persist within the humpback whale repertoire because it is closely coupled with prey biology and provides a direct benefit to the individual producing it.
Lastly, humpback whale calling behavior in the feeding grounds stands in stark contrast to singing behavior. The composition and structure of song changes within the breeding season resulting in progressive seasonal and inter-annual change 21,22,28 , whereas this study and the work of Rekdahl et al. 37 using 11 years of audio data from the east Australia migratory corridor, make it clear that portions of the call repertoire persist with time 37 . Importantly, in the east Australian population stable call types were not used as song units in adjacent years 37 , implying -similar to this study where song was not documented concurrently with calls -that calls function independently of song. Finding this commonality across ocean basins and in contrasting portions of humpback whale migratory range is quite telling. We propose that conservation of call types may be as important to call function as novelty is to song function and encourage future investigations into humpback whale calling behavior to include this hypothesis.

Conclusion
This study provides the first evidence that humpback whale call types persist across multiple generations. The longevity of calls from humpback whale feeding areas stands in marked contrast to the ever-changing humpback whale breeding-season song. Further investigation is needed to better understand the role of temporally stable calls as song units, as well as individual variation and call use across age and sex classes, social context, noise conditions, and between humpback whale populations.

Methods
We compiled acoustic recordings and associated whale sighting data from four data sets spanning 1976-2012. Acoustic data were collected using passive acoustic recording devices during summer months (June-August) on humpback whale foraging grounds throughout Southeast Alaska (Fig. 4, see supplementary material for recording equipment specifications) and paired with various forms of sighting data collection in the same locations.
Acoustic recordings from 1976 were made in Frederick Sound, Southeast Alaska (Fig. 4), were continuous during each individual recording, and were of variable length -ranging from thirty-two minutes to ninety-four minutes -and included behavioral narration on a separate recording track. Narration was summarized to glean the number of individuals present during recording. Acoustic recordings from 1997 were made in Frederick Sound and Chatham Strait as part of an ongoing investigation into coordinated group feeding in humpback whales. When possible, individuals were photographed concurrently with acoustic recordings; photographs were compared to the collaborative Southeast Alaska Humpback Whale Catalog 83 for identification. A single experienced observer (FS) pre-processed these recordings prior to the inception of this study: Visually and aurally distinguishable calls were extracted from continuous recordings and were saved individually as sound "clips" with a 0.5 second buffer at the start and end of each call. Acoustic recordings from 2007 and 2008 were collected from a cabled hydrophone in Bartlett Cove, Glacier Bay National Park (Fig. 4) with a 30-seconds-per-hour recording cycle 45 . These data were reviewed by U.S. Navy acousticians to characterize the content of each sound sample. We analyzed only the samples annotated to contain humpback whale calls. Because these data are from a remote monitoring system and not a dedicated effort to record whales, the whales were often further from the hydrophone and as such, calls with high signal-noise-ratios (SNR) were rarer in 2007 and 2008 than in other years. For this reason acoustic samples from these years were pooled for analysis. Recordings from Glacier Bay were paired with photo-identification data collected in the directly adjacent regions as part of a long-term monitoring program 55 . Acoustic recordings from 2012 in Frederick Sound were collected in approximately 30 minute increments 36 ; recordings from 2012 were made in conjunction with shore-based counts conducted from the 18.3 meter-tall Five Finger Lighthouse, Frederick Sound 41 .
Acoustic samples from the 1976 were recorded on four consecutive days and represent a limited number of individuals. These samples were not assumed to be an adequate representation of humpback whale vocal behavior at that time. Thus, while all acoustic samples were classified to assess repertoire stability over time, acoustic samples from 1976 were omitted from comparative statistical analyses of acoustic parameters. Recordings from 1976Recordings from , 2007Recordings from , 2008Recordings from , and 2012 were sampled at 44.1 kHz. Recordings from 1997 were originally recorded with a sampling rate of 22.05 kHz and were resampled at a rate of 44.1 kHz for consistency with other years. A 10 kHz low-pass filter was applied to all recordings; in the case of data from the 1990's, the filter was applied both before and after resampling. We used a 10-pole Butterworth filter with corner frequency of 10 kHz. Because the frequency band of interest was below 10 kHz, up-sampling data from the 1990's to standardize data for equivalent analysis did not result in the interpolation of data beyond the original recording range.
Spectrograms of acoustic recordings were created with Raven Pro 1.5 (Cornell Lab of Ornithology, Ithaca, NY) using a 0.093 s window length (4096 samples; filter bandwidth 15.5 Hz), Hann window, frequency resolution uncertainty of +/−5.4 Hz, and 75% overlap, and were constrained to the 10 Hz to 3 kHz frequency range to facilitate analysis. Recordings were manually reviewed in their entirety by a single experienced observer (MF). Calls were visually and aurally identified within each recording and annotated in the time-frequency domain of the spectrogram. Recordings made in 1976Recordings made in , 2007Recordings made in , 2008Recordings made in , and 2012 were prepared to match the formatting of sound clips from 1997: acoustic samples containing calls were extracted from continuous files with a 0.5 s buffer adjacent to the calls' start and end times. The SNR of each extracted sample was calculated using the method described by  84 ; to be considered for analysis, calls had to have visually distinguishable start and end points, be non-overlapping, and have a SNR of at least 10 dB above ambient noise levels 36,37,39 .
Differences in recording equipment and sampling protocol no doubt manifest in the acoustic data. In particular, because high frequency components attenuate more quickly than low frequency components, recordings made further away from the animals are likely to have lost energy in the upper ranges, and temporal patterns may be less evident. To adjust for this, high frequency acoustic parameters that may be sensitive to differences in recording equipment were not incorporated into quantitative classification analyses (e.g. high frequency, bandwidth). To further account for variation in recording equipment, salient acoustic features were extracted using the Noise-Resistant Feature Set (NRFS) measurement suite included in the MATLAB-based program Osprey 36,84 (Table 3), with the change that the de-noising step was not performed before calculating the measurement values. The NRSF was designed for detection and classification of marine animal sounds across variable noise conditions. Rather than extracting measurements from an observer-drawn annotation box, the NRSF draws a smaller time-frequency region ("feature box") in which the energy is ranked and summed within the sound relative to background noise. By doing this the loudest parts of the spectrogram have the strongest influence over the measured values, which allows for more standard measurements across recording conditions. Finer scale time-frequency measurements were made in Raven Pro 1.5 (Table 3). Fine-scale measurements were made on the fundamental frequency for harmonic sounds; for amplitude-modulated sounds containing a broadband component, measurements were made on the lowest-frequency component of the call 37,39 . To account for the mammalian perception of pitch, which is approximately logarithmic rather than linear 85 , frequency parameters were log-transformed 36,39 (Table 3). The same time-frequency parameters were input into a Principal Component Analysis (PCA) in order to aggregate variables for classification and comparative analyses (psych package) 86 . A varimax rotation was applied to maximize loading and facilitate variable interpretation 39,60 .
Calls were classified aurally and visually (AV) into previously described vocal classes, and call types by a single experienced observer (MF) using the randomization method described by Fournet et al. 36 . Using the time-frequency parameters of Table 3 and rotated principal components, a non-parametric classification and regression tree (CART) with cross-validation was run in R on the aggregated dataset to assess the likelihood that calls were correctly classified by AV analysis (rpart package) 37,38,87 . A random forest method was then performed using the same acoustic parameters with (randomForest package) 88 . These two methods are emerging as the preferred method for classification of humpback whale calls as they are robust to non-normal datasets and outliers. Further, random forest analyses improve predictive accuracy by using a bootstrapping technique that determines the OOB error, or prediction uncertainty, associated with each classification tree, rather than just one. The number of predictors randomly selected at a node for splitting was set to three, and 1000 trees were grown 38,89 .
Non-song call types can be highly variable and in some cases appear to exist along a continuum 36,39,46,49 . In the event of low classifier agreement call types with very small sample sizes (>13) were manually re-reviewed by a second observer (DC); if observer agreement was consistent calls were grouped according to AV classification. Once calls were classified, the presence of call types was compared between years. To reduce subjectivity, only call types previously described for this population were included in analysis; as the goal of this paper was not to describe new call types samples that could not be classified to known types were omitted from analysis. For this reason it is possible that call types not analyzed in this study may exhibit long term stability, or lack therefor.
Within each vocal class a set of call types that exhibited stability across four decades and were found in large enough sample sizes were selected for fine-scale comparison of acoustic parameters over time. A Bartlett's test with a significance level (α) of 0.05 indicated that, in almost all cases, the assumption of equal variance between decades was not met. To account for this and for non-normally distributed data, the non-parametric Kruskal-Wallis test was used to assess significant differences in median call parameters between decades (α = 0.05). We used a post-hoc Dunn's test with a Bonferroni correction for all relevant pairwise comparisons, in the case of ties z-quantiles were used 90 . All analyses were conducted in R version 3.3.3 87 .
The dataset analyzed in the current study is available from the corresponding author on request. Acoustics samples of each call type can be found at mfournet.wordpress.com/sounds