Contribution to unravel variability in bowhead whale songs and better understand its ecological significance

Since the first studies on bowhead whale singing behaviour, song variations have been consistently reported. However, there has been little discussion regarding variability in bowhead whale singing display and its ecological significance. Unlike the better studied humpback whales, bowhead whales do not appear to share songs at population level, but several studies have reported song sharing within clusters of animals. Over the winter season 2013–2014, in an unstudied wintering ground off Northeast Greenland, 13 song groups sharing similar hierarchical structure and units were identified. Unit types were assessed through multidimensional maps, showing well separated clusters corresponding to manually labelled units, and revealing the presence of unit subtypes. Units presented contrasting levels of variability over their acoustic parameters, suggesting that bowhead whales keep consistency in some units while using a continuum in values of frequency, duration and modulation parameters for other unit types. Those findings emphasise the need to account for variability in song analysis to better understand the behavioural ecology of this endangered species. Additionally, shifting from song toward units or phrase-based analysis, as it has been suggested for humpback whales, offers the opportunity to identify and track similarities in songs over temporal and geographical scales relevant to population monitoring.

Bowhead whales (Balaena mysticetus) are highly vocal animals and during the coldest and darkest Arctic winter months they exhibit an elaborated acoustic display that rivals the complexity of the better known humpback whale songs 1 . Although the sex of the singers and the song function in bowhead whales are still unknown, it has been hypothesised that singing behaviour plays a major role in mediating sexual interactions [1][2][3][4][5] in a way similar to humpback whales and songbirds 6,7 .
Bowhead whale songs have been described following the same hierarchical framework established to study humpback whales, with sequences of repeated units composing phrases, phrases composing themes, and themes composing a song 1 . Contrarily to humpback whales that share song within the same population, bowhead whales appear to share songs within a limited cluster of animals 5,8 . In the first review of the earliest acoustic researches on bowhead whales Würsig and Clark 1 described a 'great deal of intra and inter individual song variation' and the following studies have consistently reported song variability. The study of the variation in behavioural traits such as the acoustic repertoire can provide insights on how sexual selection acts, at individual, group, population or species levels, and shape acoustic phenotypes (i.e. songs). For example, song variations in humpback whales have been used to investigate the function of song components in the context of inter-and intrasexual selection 9,10 , to examine population connectivity or group dynamics (e.g. [11][12][13], and to investigate mechanisms underpinning gradual song evolution and cultural revolution (e.g. [14][15][16]. In bowhead whales, variability is similarly present at different levels of the hierarchical song structure. At the unit level, the unit contour shape can change through splitting, missing components, frequency shifts, or modification in the number of inflection and modulation points 2,3,17 . Unfortunately, there are currently no clear criteria establishing the acceptable level of variability within a unit type and the classification procedure remains highly subjective. Manual or automated categorisation into unit types can be furthermore complicated by the graded nature of the bowhead whale repertoire. For example, Johnson et al. 8 described 'subnotes' (subtypes of units) along a graded unit sequence (similar to gradual changes in units described in humpback whale shifting themes 18,19 ), where both unit duration and frequency parameters evolved along the song. The authors eventually www.nature.com/scientificreports/ assigned the last units of the sequence to a different type, based on "trend of frequency", a criterion that has not been reported in other song studies. At the song level, one of the most variable components appears to be the number of repetition of units per phrase. All studies describing song acoustic and structural parameters report on a high level of variation, with repetition rates ranging from one to tens of times, and some units being sometimes entirely omitted or rarely present in the song 3,8 . Directly linked to unit repetition, variability across renditions of the same song have been previously reported including variability in song duration along a song bout 2 and the occurrence of 'primitive versions' during the initial development of a song as well as song variants based on unit or phrase repetition, order or absence 3,8 . Similar variants have been recently identified as "song fragments" in humpback whales songs 20 . Collectively, the bowhead whale acoustic studies indicate that variability at multiples levels is an important part of the bowhead whale singing behaviour. Although this variability has been consistently reported either as a qualitative descriptor of songs 1,21 or with a statistical measure of dispersion in previous publications 2, 3,8,17 , the study of variability itself has never been the central subject of bowhead whale acoustic research.
Research challenges related to the bowhead whale Arctic distribution have constrained the advances of acoustic studies and although the first reports of bowhead whale songs date back to the early 1980s there have only been a handful of publications since then describing their songs. Many of those studies have focused on song diversity using the number of different songs as a single metric 2,3,5,8,21,22 . At first glance, counting songs might seem straightforward, but defining at which point a song is different enough to be considered new represents a delicate problem when the signals to be categorised display high variability. Furthermore, integrating variability in the definition of what is new or different is essential to delve into the potential functions of a song and advances our understanding of singing behaviour in bowhead whales.
Here we used a one-year dataset that sampled bowhead whales from the critically endangered Spitsbergen population in a previously undescribed overwintering area in Northeast Greenland coastal waters to examine song variability. We selected songs shared by multiple individuals as a basis to focus on inter and intra individual aspects of variability. The shared songs, displaying comparable song structure and unit acoustic characteristics, were used to identify patterns of variations using semi-automated methods. We explored the potential ecological and behavioural significance of this song variability within the bowhead whale mating context.

Methods
Study site and data collection. Acoustic data were collected in coastal waters of Northeast Greenland (78°30′N and 10°0′W) on the continental shelf about 170 km offshore, close to the Northeast Water (NEW) polynya (Fig. 1). The NEW is an area of major interest as a summering ground for the Spitsbergen bowhead whale population 23 .
The recording equipment was deployed during the Oden Artic Technology Research Cruise 2013 on the 23th of August 2013 and retrieved on the 17th of September 2014. The hydrophone, an AGUAtech Low-Power Scientific Measurement Hydrophone (sensitivity -160 dB re 1 V/μPa), was directly connected to the acoustic recorder, an RTSys EA SDA14 (RTsys, Caudan, France). The system was suspended a few tens of metres above the sea floor using a subsurface float, in a water depth of 227 m. The recording was set on a duty cycle of 2 min on and 30 min off, at a sampling frequency of 39,062 Hz in 24 bits. The data was hosted in the LIDO (http://liste ntoth edeep .com) framework and software package 24 . Data annotation and measurements. The dataset was first subsampled to 25% (one out of four 2-min files) and then explored to detect the occurrence of shared songs. Files were aurally and visually analysed by an expert in Adobe Audition (Adobe Systems Inc., San José, CA, USA). Each 2-min file was displayed as a spectrogram (2048 point FFT, 75% overlap, Hann window for a frequency resolution of 19 Hz and a time resolution of 52 ms). The songs shared by multiple animals were identified based on two criteria: the songs overlapped in time and contained similar unit types temporally arranged in similar phrases. Figure 2 shows an example of three units with similar acoustic characteristics (frequency contour shape, medium frequency) overlapping in time with each other (i.e. a second unit starts while a first unit has not ended yet), and associated with lower frequency units also sharing between them similar acoustic characteristics. This indicates that three individuals are sharing similar songs.
Two unit types were consistently present in most of the detected shared songs: a long whistle-like unit displaying multiple inflections (thereafter referred to as M unit type) and a short whistle-like unit usually containing one inflection (S unit type). M and S units were temporally associated in a phrase, that was usually repeated multiple times. Subsequently, the analysis was focused only on the songs containing M and S units. All the song occurrences containing the M and S units were identified from the subsampled dataset. Songs were manually classified into song groups. Song groups were defined as songs presenting a high level of similarity on the following aspects of the song characteristics and structure: overall similarity of the frequency contours of the units M and S, presence of others unit types, order of units, and biphonation characteristics if present (overview given in Table 1). The units were named after the overall shape of the frequency contours, and units belonging to a particular song group were designed by adding the song group to their name (e.g. MSG11 refer to the units M from the song group SG11).
Manual analysis was then conducted on all files from the same days where the songs were initially identified to obtain more occurrences of these song groups. All occurrences, except the ones with a poor signal-to-noise ratio, were selected for further analysis. The minimum number of occurrences for a song group was three, the maximum was eight. www.nature.com/scientificreports/    25 ). Some of the measured units displayed biphonation, a nonlinear phenomenon consisting of the simultaneous production of two independent frequencies 26 . Due to degraded signal to noise ratio in the lower frequencies, biphonated components were not measured in Raven, but an estimation of their frequency bandwidth and duration were noted to allow comparison and define biphonation types. The unit measurements consisted of frequency, modulation and duration parameters that are commonly used in literature describing or classifying animal vocalisations. These parameters included minimum frequency (1), maximum frequency (2), frequency bandwidth (3), frequency bandwidth 90% (4), and duration (5). Unit contours were extracted using the Frequency Contour Measurements tool in Raven. The extracted contours were imported in Matlab v.R2017a (MathWorks, Inc., Natick, MA, USA) and a custom-made smoothing function (moving average) was applied to remove noise. Additional parameters (frequency start (6) and frequency end (7), number of inflections (8), median frequency (9)) were automatically extracted in Matlab from the smoothed contours. To describe the distribution of the parameters in each unit, the mean, standard deviation (sd) and coefficient of quartile variation (cqv) were used. The cqv was preferred over the coefficient of variation to take into account that some distributions were leptokurtic 27 .
Data analysis. After a normalisation process (centre by median, scale by range), the nine parameters described above were used as features (Orange data mining toolbox version 3.18.0, 28 ). Multidimensional scaling maps (MDS), based on pairwise similarities, were produced to examine the pattern of acoustic variation in units. MDS was used to provide an automatic way of investigating visually the grouping of units into clusters and to determine whether unit types fell into distinct subtypes without stipulating a priori categories. MDS analysis allows to view the clustering tendency of the data in a few dimensions, without forcing a certain number of clusters as is needed for many clustering algorithms. Unit clusters from the manual analysis were plotted against the MDS map to explore further the acoustic similarity/dissimilarity patterns. Kruskall stress and R-squared values were used to assess the quality of the MDS model 29,30 . Higher R-squared values indicate that the distances between points in the high dimensional feature space correlate well with the distances between the same points in the MDS projection. Similarly, a stress value inferior to 0.05 suggest that the goodness of fit of the MDS distance model can be considered as good 29 . Correlation between the selected features and the MDS dimensions was used to provide an interpretation of the MDS dimension. In order to compare the results to a fully automated clustering method, a hierarchical clustering (average linkage, based on Euclidean distances) was performed on the dataset. To take into account the higher variability regarding the M units, the dataset was split into two subsets, the M units and all the other units (named Except-M). The hierarchical clustering algorithm was run on each dataset. Cluster solutions from the semi-automated (MDS map) and fully automated method (hierarchical clustering) were assessed for performance using a cluster validity measurement, the Silhouette index. This index is an indication of clustering quality and takes values between -1 and 1. The higher the silhouette value is, the better is the quality of clustering, i.e. the intra-cluster dissimilarities are small compared to the inter-cluster dissimilarities 31,32 . www.nature.com/scientificreports/ www.nature.com/scientificreports/

Results
Bowhead whale songs were detected from the second week of October 2013 to the first week of April 2014, with almost daily detections during the peak of the season from December to February ( Supplementary Fig. S1 online). These detections indicate that a probable bowhead whale mating ground described from the Fram Strait data 22,33 could actually extend some 170 km further east towards the coast of Greenland. The manual analysis of the songs resulted in the identification of 13 song groups that comprised songs sharing similar acoustic characteristics (i.e. the presence of specific unit types M and S (see Methods and Fig. 3). The analysis revealed that up to 3 animals were singing similar songs simultaneously (Fig. 2). The song groups were detected between November and February (see Supplementary Fig. S2 online). An example of each of the song groups is shown in Fig. 3 and can be listened to on the LIDO website (link to song groups).

Song groups description.
A total of 804 units were measured and manually classified into 9 unit types using visual and audio properties of the units as well as the temporal arrangement of units in a song. An overview of the units and their main characteristics is provided in Table 2 and Supplementary Fig. S3 online, with their occurrence in song groups summarised in Table 1. Songs contained 2-6 unit types ( Fig. 3 and Table 1). The overall structure of the songs consisted of two main phrases. An M unit and a variable number of S units composed the first phrase, present in all the song groups. A second phrase consisting of two alternating short and narrowband low frequency units (sLFdown and sLFconst/sLFconst2) were detected in 6 song groups, the higher frequency sLFconst2 replacing sLFconst in half of them. Other units, such a s Vhigh and rumble-like units were present in different song groups, with the exception of sup occurring only in SG4.
Semi-automated clustering. The MDS map globally reflected the output of the manual analysis in terms of clustering the units into separate types. Interpretation of the distances and relative positions of the clusters from the two-dimensional MDS plot was regarded as valid based on the values of R squared (R 2 = 0.99) and goodness of fit (stress = 0.066). Close relationships (R 2 > 0.7) were found between the first MDS axis and the variables median frequency, high frequency, low frequency, contour start and frequency bandwidth, suggesting that this dimension was strongly related to the global positioning of the unit in the frequency domain. The second Table 2. Summary of key acoustic characteristics of 9 unit types of bowhead whale songs, including mean ± standard deviation, and coefficient of quartile variation (cqv in %). www.nature.com/scientificreports/ MDS axis was less correlated to the acoustic variables, although correlations with inflections and duration were noticeable (R 2 = 0.3-0.4) (Fig. 4). Units defined by the manual procedure appeared to be well separated on the MDS map. Short units (sLFdown, sLFconst, sup and sLFconst2) were closely positioned along the first dimension, i.e. they separated along a frequency gradient, respectively from lower to higher frequency (Fig. 4). Analysis at the song group level of the S cluster did not reveal any clear structuring although in some SGs the units showed a tendency to aggregate in the same spatial location in the cluster (e.g. SG4, Fig. 4 and Supplementary Fig. S4a online). Vhigh clusters presented a higher degree of scattering that could reflect the inaccuracy of measurements due to lower SNR of this unit type. Rumble were displayed along a continuum on the second dimension of the MDS map, with two subtypes based on duration parameters (short < 5 s and long > 5 s).
Unlike the densely grouped S units, the M units were more scattered. In the context of song sharing, this unit type was manually defined as a broad category, and high variability on some acoustic parameters (especially duration, inflection and frequency bandwidth) was consistently observed. The MDS projection suggested that M units formed five distinct subtypes ( Fig. 4 and Supplementary Fig. S4b online). The map showed one large and scattered cluster (Mo) containing 68% of the M units from 9 song groups overlapping each other, and 4 smaller separated clusters, each one containing M units from a specific song group. MSG11 was characterised by a very high number of inflection points and a reduced bandwidth (respectively 17 and 501 Hz, Table 2). The second cluster grouped MSG5 units and shared similar characteristics with MSG11 in terms of duration but displayed a three-time higher bandwidth, and only half of the inflection points. MSG1 had the longest duration, and intermediate values between MSG11 and MSG5 for inflection points. MSG2 was grouped together based on higher values for all parameters compared to the large cluster Mo. Inside the 5th cluster Mo containing all the other SGs, fine scale structuring was visible for some song groups, as noted during the manual analysis. Interestingly, SG7 was positioned at the Mo cluster edges (Fig. 4 and Supplementary Fig. S4b online), in a transition position between Mo and the MSG11 clusters and visual observation of the MSG7 unit contour shapes confirmed that MSG7 shared characteristics both with the large Mo cluster (large bandwidth, reduced number of inflections) and MSG11 cluster (reduced bandwidth, high number of inflections).
Standard deviation (SD) values for the M unit acoustic parameters were higher than for the other unit types, and still superior to S type values when considering only the Mo cluster (Table 2 and Fig. 5a). Short units had low SD values. Globally, the coefficient of quartile variation values per parameter showed that inflection, duration and frequency bandwidth were the most dispersed around the mean (Table 2 and Fig. 5b).
Clusters from the MDS map were used in combination with the manual analysis to create an improved classification of unit types and its validity was assessed using the silhouette coefficient. Values of 0.49 for both the  Fig. S5 online). A higher silhouette value represents a better clustering solution where each object (unit) is closer to its own cluster centroid. Silhouette values from the MDS map clusters were also higher over a range of cluster numbers, with one exception that clearly oversimplified the grouping of the units into only 3 or 4 distinct clusters.

Biphonation patterns on M and S units. All the analysed song groups presented a biphonation on one
or two unit types. Biphonation presented strong similarities between some song groups, and could be divided into 3 and 4 major types for M and S units respectively (Supplementary Table S6 online). Close ups of the units are shown in Supplementary Fig. S7 online. The main recurrent biphonation type on M and S units was a tonal frequency modulated downsweep, respectively present in 6 and 7 song groups. Mbi1 overlapped in time with the whole M unit, and sometimes even with the following S unit. Sbi1 tonal downsweep sometimes included frequency modulation. Remarkably different, consisting in short broadband "whoops" in sequence, Mbi2 presented a highly consistent and repetitive pattern on SG11 and SG13 but variable temporal characteristics in other song groups (SG2, SG8 and SG10). A similar biphonation pattern was present on S unit from SG1 (Sbi2). M from SG4 showed a distinct type of biphonation, a short constant tonal around 130 Hz with 3 harmonics, present at the end of the unit. Altogether, biphonation patterns resulting from the combination of biphonation types on M and S units were unique for almost half of the song groups.

Discussion
Bowhead whales singing behaviour has been studied for more than 3 decades and an increasing diversity in their song repertoire has been described from the simple patterned moan sequences reported by Ljungblad et al. 34 to the remarkably high song diversity described for the Spitzbergen population 22 . Although previous studies provided song or unit descriptions, there has been little discussion about variability and consistency in bowhead whale singing display and its impact on our understanding of bowhead whale acoustic ecology in a mating context.  21,22 , also described intense singing activity of bowhead whales in winter months, suggesting that this location was a mating ground for the Spitsbergen population 33 . Our data, recorded approximately 170 km from the Fram Strait, provides compelling acoustic evidence that the spatial range of the mating ground of this critically endangered population would likely extend from the western Fram Strait to the Northeast Greenland coastal area, hundreds of kilometres more west than previously reported.
Past studies describing the bowhead whale repertoire from different populations have reported the occurrence of the same song overlapping in time, demonstrating that bowhead whales do share songs 2,3,8 . Elaborating on this finding, Johnson et al. 8 estimated that during the migration of whales from the Bering-Chukchi-Beaufort population, a song type could be shared by up to nine animals. Temporal overlapping of songs in our data revealed that at least three animals were sharing the same song (Fig. 2). Structural characteristics of these shared songs were used to identify song groups based on similar unit types arranged in comparable temporal pattern, leading to the selection of 13 song groups. Analysis of the song groups revealed variability at multiple hierarchical levels of the songs, from units to phrases.
Variations in a given unit type have been described in previous acoustic studies of bowhead whale songs. Evolution of unit structure (loss of part of the unit, splitting into smaller units) or frequency (e.g. shifting) has been reported from all the populations that have been acoustically studied so far 2,3,8 . Our analysis revealed two major paths for intra-unit type variability, that were measured in terms of standard deviation and coefficient of quartile variation and were readily noticeable on the MDS maps. The first path of variation was related to duration and inflection points and was especially noticeable on two unit types, M and rumble. These two unit types both had a distinctive frequency contour displaying multiple inflections, thus variation in duration was coupled with variations in the number of inflection points. The second path was linked to variability in frequency parameters with variations mainly related to frequency bandwidth. Gradual changes on one or both of these paths of variation eventually led to the appearance of new unit types, and therefore larger vocal repertoire. Novel or complex vocal repertoires have been identified as important traits influencing female sexual selection in several studies on birds (e.g. 35,36 ) and in a similar way, sexual selection might be one of the mechanisms underlying innovation in humpback whale song 37,38 .
While ecological and behavioural functions of songs in bowhead whales are still unclear, the presence of intensive singing during the reproductive season indicate a function in sexual advertising. This would correspond to a polygynous mating strategy as seen in humpback whales (where one male mates with multiple females), although the exact nature of the bowhead whale mating system is still unclear to researchers 17 . In this context, a possible advantage of displaying calls with longer duration could be linked to an increase in stamina with age or as an honest advertisement of physical fitness as it has been described in bearded seals, Erignathus barbatus 39 . Variability in inflection points could similarly be associated to fitness since highly modulated sounds require higher motor control of the vocal apparatus. Interestingly, for most of the unit types, variability was contained inside a unique cluster corresponding to a single unit type (i.e. unit types were separable in the feature space as they had little to no overlap). However, combined variability on both frequency and modulation/duration parameters led to a gradual split of M units into subtypes, some of them matching the song groups labels. The M unit type, possessing characteristic of long duration, numerous inflections points and variable bandwidth could www.nature.com/scientificreports/ be identified as the unit that best represents abilities of the animal to produce physically challenging sounds. This would be similar to trilled songs in songbirds where production of trills, expressing vocal performance, defines a constraint on the song, represented by a trade-off between trill rate and frequency bandwidth [40][41][42] (Fig. 6). The upper boundary of this trade-off represents a vocal performance limit and songs containing trills with acoustic characteristics close to this limit are considered to express high vocal performance [40][41][42] . In the context of sexual selection, vocal performance has been related to male quality assessment by female birds and to male-male competition 40,43 . It appears that highly modulated sound produced by bowhead, such as the M units studied here, face a similar trade-off between number of inflection points (analogue to trill repetition rates) and the frequency range they can span (Fig. 6). Song groups containing M units at the limits of the performance constraints, such as SG1, SG5 and SG11, could then be associated with higher vocal performance (Fig. 6). If songs are used in a mate choice context, the presence of units with numerous inflections could function as an honest indicator of mate quality, by virtue of the "inflection rate/frequency band" trade-off. As such, these units could carry a higher level of inter-individual variability expressing different levels of fitness, and additionally higher levels of intra-individual variability in relation to the difficulties of accurately producing and reproducing such sounds. An additional source of unit-level variability described here comes from the presence of biphonated signals. Recent acoustic studies confirmed that bowhead whales use biphonation 4 and that the use of biphonated signals can be quite common in bowhead whale songs 8 . All the song groups studied here showed biphonation on one or two unit types, with the occurrence of different biphonation types including continuous whistle-like sounds, short impulsive broadband sounds and harmonics signals. Functions of biphonated signals remain unclear and they have been suggested to encode information on size 44 , individuality [45][46][47] or localisation of the signal emitter 48,49 . Mastering the production of biphonated signals, which require high vocal control of the laryngeal apparatus (Musculus diverticuli laryngei) 50 , could function as another honest signal of mate quality in a context of reproductive display 4 . If one structure is the source of two simultaneous sounds, then producing two very different sounds i.e. decoupling mechanisms of air movements in the larynx, would require strong vocal skills. Broadband discrete impulsive signals coupled with continuous M units could be more difficult to produce than whistle-like types. As such, the presence of type Mbi2 could indicate older or more experienced singers and, interestingly, this biphonation type occurred concurrently with M units than can be associated with higher vocal performance (SG11). Similarly, recent studies on bowhead whales' closest relatives, the North Atlantic right whales, have linked age with increasing structure in non-linear phenomena, from disorder (deterministic chaos) to increased control (biphonation and subharmonics) 51 .
On the upper level of the song hierarchy (phrases), temporal arrangement of units in phrases showed regular patterns of association between units in stable phrases. However, while phrases containing M and S were present in every song group, phrases containing the short alternating units (sLFdown, sLFconst, sLFconst2, sup) were only present in half the song groups. The duty cycle used in this study (2 min on, 38 min off) could have resulted in recording only partially the variability of each song group. The average song duration of a bowhead whale song is considered to be approximately one minute 1,3,8,21 but songs lasting almost two minutes have been reported. Therefore, in our study, most of the time, only one entire song was present in the file, and variability at the song bout level (differences in successive iterations of the same song by the same animal) could not be assessed. Nevertheless, the aforementioned phrases containing the short alternating units were present in half of the studied song groups, and can be considered as common and part of the song structure. In humpback whales, the only other cetacean species to produce such complex songs, omission or insertion of entire themes is considered part of the natural evolution of song types and themes present in 50% or more of the songs were included in song sequence study 14 . Similarly, temporal evolution of a song type over a period of a few months has been previously described in bowhead whales. Delarue et al. 3 reported the observations of primitive versions occurring during the initial development of the song until the final, stable version emerged. The overall song variability suggested a low level of constraint in the song structure (e.g. units missing, variability in the repetition rate of units, presence or absence of a phrase). In humpback whale songs, phrase structure is reported to be much more flexible than in bird songs and can even be progressively modified within individual songs 14,52 . While there are inherent recording limitations when conducting acoustic studies in the high Arctic, longer duty cycles are required to quantitatively describe song structure variability and evolution over a season.
To describe and quantify variability at unit level and to allow the description of song structure, we have used a semi-automated method combining multidimensional scaling maps and human expertise. MDS maps showed unit clusters corresponding well to human classification into nine distinct unit types. Based on this agreement, MDS was used to investigate subclustering in the M unit type. Units were positioned along a continuum on the MDS axes, following a gradation in the frequency, modulation and duration parameters of this unit type, which complicated the classification into subtypes. Graded vocal repertoire has been described for many mammalian taxa such as primates 53 , elephants 54 , sirenians 55 , odontocetes [56][57][58] , and mysticetes including humpback 59 and right whales 60,61 . Graded repertoires are challenging for classification tasks since there is no clear boundary where a grouping begins or ends. Manual (human) classification of animal vocalisations usually leads to accurate results 62 , but with a strong subjective component, whereas automatic methods can add objectivity and standardisation. It appears that a combined method such as the one we described here can greatly help the classification process at the unit level. The next step forward towards standardisation of bowhead whale song studies would be applying automated methods that have been validated on other species. For example, the Levenshtein distance, a similarity analysis that compares vocal sequences, has been successfully applied to humpback whale songs to objectively quantify differences among groups of songs [63][64][65] .
Earlier studies suggested that bowhead whales, unlike all other singing whale species, entirely change their song repertoire from year to year. This feature has been reported for the Bering-Chukchi-Beaufort (BCB) population on their migration route 1,3 and for the Eastern Canada-Western Greenland (EC-WG) population in Disko Bay in spring 5 . It is likely that these studies only sampled a small part of the song repertoire, neither locations www.nature.com/scientificreports/ being wintering grounds where most of the singing is thought to occur in relation to reproduction 66 . A recent study from the Spitsbergen population analysed data from a probable mating area in the Fram Strait during four successive winter seasons. As predicted, researchers recorded and described a much higher diversity of songs. Again, the results highlighted the lack of recurrence of song types from year to year, indicating that bowhead whales completely renew their singing repertoire each year 22 . Nevertheless, our results suggest that the high variability reported in this study could hinder the recognition of units and phrase types that are potentially reused from year to year, in a dynamic and progressively evolving singing display. Although songs may be considered different between years (depending of the level of variability researchers select to assess song diversity), some unit types and patterns of association between unit types may be maintained over time. Acknowledging this inherent variability could help to track patterns in units and song presence at different temporal and geographical scales. For example, close examination of songs from the Spitzbergen population from our dataset and from the Stafford et al. 22 study revealed that identical units forming highly similar songs were present at the two locations (170 km apart), during the same time period (December 2013-January 2014) and moreover, were present at the same location in the Fram Strait during two consecutive years (Fig. 7). Furthermore, we discovered the presence of highly similar units in our dataset from Northeast Greenland and in published song studies from the BCB population. The MSG11 unit described here is identical to the warble unit characterised in Delarue et al. 3 and to the high-frequency whistle reported in Würsig and Clark 1 . Figure 8 shows that units from theses 3 studies share similar median frequency (around 1500 kHz), shape (multiple inflections), and an identical highly specific biphonation pattern (broadband grunt or "whoop"). Additionally, the 1988 song described by Würsig and Clark 1 is composed by the aforementioned high-frequency whistle associated with another unit that closely resembles the S unit described here, both in terms of frequency parameters and contour shape (Fig. 8).
These findings are highly relevant in the context of population exchanges questions that arose from recent studies of the Spitzbergen population 22,23 . Whether the individuals recorded from the Greenland Sea are remnant of the nearly extirpated Spitzbergen population or immigrants from other populations (BCB or EC-WG) is not known. Drastic sea ice loss in the high Arctic has already made contact possible between once ice-separated bowhead populations and the continuous increase of this phenomena is likely to open new areas for bowhead whale migration 67 . Tracking acoustic similarities at the unit or phrase level will help to gather acoustic evidence revealing these movements.
Bowhead whales are undoubtedly amongst the most prolific singers, and the high diversity of the singing repertoire has just begun to unveil. In order to build up our understanding of bowhead whale acoustic ecology, we propose that future acoustic studies should consider (i) further developing automated methods using objective categorisation to facilitate inter study comparisons, (ii) shifting the framework of acoustic studies from a song-based to a unit or phrase-based approach as it has been suggested for humpback whales 19 to investigate song similarities in the context of song sharing, song evolution and cultural transmission, (iii) conducting collaborative studies to track acoustic patterns at larger temporal and geographic scales.  The blue box corresponds to the long whistle-like vocalisation, with multiple inflections (MSG11 described in this study). The red box shows the associated biphonation, a rapid sequence of shorts grunts or "whoops".