Body density of humpback whales (Megaptera novaengliae) in feeding aggregations estimated from hydrodynamic gliding performance

Many baleen whales undertake annual fasting and feeding cycles, resulting in substantial changes in their body condition, an important factor affecting fitness. As a measure of lipid-store body condition, tissue density of a few deep diving marine mammals has been estimated using a hydrodynamic glide model of drag and buoyancy forces. Here, we applied the method to shallow-diving humpback whales (Megaptera novaeangliae) in North Atlantic and Antarctic feeding aggregations. High-resolution 3-axis acceleration, depth and speed data were collected from 24 whales. Measured values of acceleration during 5 s glides were fitted to a hydrodynamic glide model to estimate unknown parameters (tissue density, drag term and diving gas volume) in a Bayesian framework. Estimated species-average tissue density (1031.6 ± 2.1 kg m-3, ±95% credible interval) indicates that humpback whale tissue is typically negatively buoyant although there was a large inter-individual variation ranging from 1025.2 to 1043.1 kg m-3. The precision of the individual estimates was substantially finer than the variation across different individual whales, demonstrating a progressive decrease in tissue density throughout the feeding season and comparably high lipid-store in pregnant females. The drag term (CDAm-1) was estimated to be relatively high, indicating a large effect of lift-related induced drag for humpback whales. Our results show that tissue density of shallow diving baleen whales can be estimated using the hydrodynamic gliding model, although cross-validation with other techniques is an essential next step. This method for estimating body condition is likely to be broadly applicable across a range of aquatic animals and environments.


Introduction
The body condition of animals influences survival rate and reproductive success and thereby impacts the dynamics of entire populations. Body condition also affects an animal's behavioural decisions related to foraging, predator avoidance, migration, and reproductive strategies (e.g. [1][2][3]). Many marine mammals undergo substantial changes in lipid-store body condition as a result of annual fasting and feeding cycles [4,5]. For migratory species, the cost of reproduction at breeding grounds is supported by energy gained on feeding grounds. Thus, the amount of energy stored during a feeding season strongly influences reproduction via pregnancy rate [6], foetal development [7], body condition and survival of offspring [8][9][10][11] and the competitive capabilities of males. It is also likely that body condition influences the foraging decisions made by baleen whales relative to where prey items are located in the water column [12,13]. Because body condition is an important factor affecting fitness, measuring body condition of free-ranging cetaceans is essential for understanding their ecology as well as for designing effective conservation plans [14][15][16].
Baleen whales (parvorder Mysticeti; order Cetartiodactyla) are a group of marine mammals that cycle fat stores on an annual basis, substantially changing their appearance, behaviour, and fitness [17]. Given these dramatic changes, developing methods to quantify their body condition in the field has great value. Traditional approaches to examine variations in body condition and energy store of baleen whales involved anatomical measurements that were often made in conjunction with whaling operations [4,18,19]. Blubber thickness of whale carcasses has been used as a proxy of body condition [4,18,20], since most of the energy is stored in the form of blubber [21] although a considerable amount of energy is also stored in muscle and intra-abdominal fat [4,21]. Blubber lipid content of whale carcasses was also important in the assessment of condition of cetaceans [19]. As the thickness as well as lipid content and fatty acid composition of blubber has been shown to vary across the body of cetaceans, multiple-site measurements of blubber thickness are particularly useful to examine total body condition of cetaceans [22][23][24]. Many studies have investigated seasonal trends in energy storage of several species of baleen whales by means of blubber thickness and morphometric data, reporting that seasonal fattening varies with different sex and age classes, reproductive stages, as well as prey availability [4,6,25].
Although carcasses have provided many insights into the physiology and body condition of baleen whales, a key limitation is that temporal changes of the same individual cannot be measured. Also, studies using carcasses may not be widely applicable to cetaceans because they require lethal sampling or collection of samples from stranded animals or fisheries bycatch. To collect blubber and other tissue samples from free-ranging cetaceans, biopsy darting is commonly used where modified dart tips are delivered using a crossbow or a pneumatic rifle [26]. The percentage lipid content of blubber from carcasses is considered to be an informative measure of fattening [19]. However, the biopsy blubber samples may not be useful to measure body condition of free-ranging cetaceans because (1) the force of darting can damage adipocytes causing lipids to be squeezed out of samples, or to seep out of blubber biopsies while in seawater [27], and (2) the sample only penetrates a short distance into the blubber layer. In addition, it is difficult to obtain multiple biopsy samples from a whale whose blubber thickness and composition vary across the body [22][23][24]. Visual assessment of external shape and appearance based on boat-based photographs has been used for evaluating body condition of right and grey whales [28,29]. Photogrammetric measurements of body width, reflecting blubber thickness, using vertical aerial photographs taken from aircraft or unmanned aerial vehicles has also been used to assess nutritive body condition of some whale species [11,30,31] although measurements of such body shape patterns are limited to the visible 2-dimensionall shape of surfacing of whales, and may not be suitable for other more cryptic species.
An alternative approach is to use body density of diving animals as a proxy of lipid-store body condition [32]. Lipids are less dense than seawater while other non-gas body components are denser than seawater. Body composition, particularly the ratio of lipid to lean tissue, therefore strongly influences body density and hence the buoyancy of diving animals [5]. It has been shown that buoyancy influences swimming behaviour and energetics of diving animals [33]. For instance, buoyancy forces affect stroking efforts [34,35] and swimming patterns, with more gliding occurring in the direction aided by buoyancy [15,34,[36][37][38][39].
Buoyancy also influences gliding performance by altering vertical speeds during inactive drift periods [5], prolonged glides [39] or short-duration glides [38]. This effect of buoyancy on gliding performance has led to the development of tag-based methods to quantify the body density of diving animals via hydrodynamic analysis. This approach was first developed for free-ranging elephant seals (Mirounga spp.): body density was quantified by analysing the vertical speed during inactive drifting periods (i.e. drift rate) at which the buoyancy force is assumed to be equal to the drag force [5]. The drift dive method has proven useful for longterm monitoring of body lipid-stores in elephant seals providing new insights into when and where they gain or lose lipid stores [5,40,41]. However, use of the drift dive method is limited to a few pinniped species that routinely perform drift dives (Mirounga angustirostris [34,42]; M. leonina [5,40]; Arctocephalus forsteri [43]; Cystophora cristata [44]). Gliding during the descent or ascent phase of a dive, on the other hand, is commonly observed across a range of diving taxa [33,37]. A more widely applicable approach, the glide model, was introduced by Miller et al. [38] to estimate body density of sperm whales using a hydrodynamic glide model that predicts how drag and buoyancy forces influence acceleration (or deceleration) during short-duration glides. Aoki et al. [34] conducted a validation analysis using isotope dilution and confirmed a strong correspondence in body density estimates of elephant seals obtained from the drift dive and the glide models [45].
In the glide model, acceleration during a glide is determined by the difference between drag and net buoyancy forces along the swimming path of the animal [38]. The force of non-neutral buoyancy or 'apparent weight' (difference in mass of the diving animal and the displaced water) acts vertically on diving animals, and depends on the density of body tissues as well as the volume of air carried within the body (the diving gas volume). While body tissues are relatively incompressible at depth, the volume of air in the body progressively decreases with increasing depth, thought to closely follow Boyle's Law for marine mammals [46]. Thus, tissue-derived buoyancy can be separated from air-derived buoyancy when gliding data is available over a wide depth range. To date, the glide model has been demonstrated to be useful to estimate the body density of several species of marine mammals, including elephant seals [34] and some deep diving toothed whales (sperm whale, Physeter macrocephalus [38]; Northern bottlenose whale, Hyperoodon ampullatus [15]; long-finned pilot whale, Globicephala melas [47]) that routinely perform dives deeper than 200 m where the effect of air-derived buoyancy is considered to be negligible [38].
In this study, we apply the hydrodynamic glide model to estimate body tissue density of humpback whales (Megaptera novaeangliae) in two geographically distinct feeding populations (the Gulf of St Lawrence, Canada and the Western Antarctic Peninsula, WAP). In comparison with deeper diving toothed whales, humpback whales may not seem ideal candidates for the glide model because they routinely dive only to relatively shallow depths at which gas volumes are likely to more strongly influence net buoyancy. For example, the mean dive depth per tag record in this study ranged from 22.8 to 180.8 m, with the deepest dive recorded being 388.3 m. Apart from a shallower diving depth range, humpback whales tend to dive and glide at relatively shallower pitch angles, requiring the generation of lift. The large flippers of humpback whales are well-suited for this purpose [48], but the need to generate substantial lift forces may raise concerns about the applicability of the glide model because the current model does not include the potential effect of lift-induced drag which was shown to be negligible in deep divers that maintain steep pitch during glides [34].
The objective of this study was to examine whether the hydrodynamic glide model can be applied to shallower diving baleen whales by examining the precision of body density estimates obtained from a narrow depth-range dataset. Our results show that we were able to obtain estimates of humpback whale body density using this method. Though the precision of the estimates was not as fine as was previously reported for a deep-diving toothed whale [15], the precision of individual body density estimates was substantially finer than the variation across different individual whales, including some differences between the geographic locations where tags were attached. We conclude that the glide method has potential to be used to track the body condition of shallow diving baleen whales, enabling future applications as a tool to study their health and how body condition relates to reproductive status, animal behaviour and the influences of environmental change and variability.

Data collection
Field studies were carried out at two geographically distinct summer feeding grounds of humpback whales (Megaptera novaeangliae): the Gulf of St Lawrence in Canada (49.7-50.0N, 63.4-65.0W) and the western side of Antarctic Peninsula (64.4-64.9S, 61.8-63.1W). Animal-borne archival tags used in the study were either 3MPD3GT loggers (Little Leonardo Co., Tokyo, Japan) or sound and movement recording DTAGs ( [49]; Table 1). The 3MPD3GT loggers were programmed to record depth, temperature, flywheel swim speed and 3-axis magnetism at 1Hz, and 3-axis ± 3 g acceleration at 32Hz. The DTAG sampled pressure and a 3-axis ± 2 g acceleration at 50Hz, which was later downsampled to 5Hz. The 3MPD3GT loggers have the ability to measure flow speed using a front mounted impeller (flywheel). To ensure that speed is measured in the direction of travel, 3MPD3GT tags are mounted in hydrodynamic (tear shaped) floats with a single suction cup mounted at the anterior end, and vertically mounted tail fin at the posterior end. The location of fin and suction cup ensure that the force acting on the tag cause the tag housing to swivel on the animal and orient into the direction of flow. DTAGs are attached to the animal with four suction cups. Tagging was conducted from rigid-hull inflatable boats and either a 5 m or an 8 m handheld carbon fibre pole was used to attach the tag.

Analysis of tag data
Pressure data recorded by archival tags were converted to absolute values of hydrostatic pressure using calibration values and converted to meters. A dive was defined as any submergence to a depth of > 10 m. Dives were broken into descent, bottom and ascent phases based on changes in pitch following Miller et al. [38]. As tags were attached to whales at random orientations, the 3-axis acceleration data recorded by the tags was converted to a whale-centred, whale fixed reference frame (whale-frame) using established methods [38,49]. The accelerometers recorded both specific (e.g. stroking) and gravity-based accelerations (i.e. changes in response to posture change). Under the assumption that changes in the posture of the tagged whale occurred at lower frequency than changes in accelerations resulting from body motions such as thrust, a frequency-based filter (low-pass finite impulse response filters with tag-specific thresholds set at 0.12-0.15 Hz) was applied to the entire acceleration time-series to separate these two components. Then, pitch and roll angles of the whales were calculated from the low-frequency component of accelerations [37,39,50], while the high-frequency component was used to identify stroking versus gliding periods. For 3MPD3GT dataset, stroking was identified when oscillation on the high-frequency component of surge accelerations indicating fluke beats exceeded a threshold that was set for each deployment (0.1-0.2 m s -2 ). Speed sensor data was visually inspected to confirm the presence of stroke-derived acceleration. For DTAG dataset, stroking was detected using high-frequency accelerations at both surge and dorso-ventral axis with thresholds set for each deployment and each axis (0.1-0.2 m s -2 ). Gliding periods were automatically detected as the period when the tagged animals did not stroke. The speed sensor of the 3MPD3GT logger recorded swim speed as the rotation of an external impeller mounted on the anterior end of the logger, which correlates linearly with the speed of water flow passing through the impeller. The rotation rate (number of rotations per second) was converted to speed (m s -1 ) using a calibration line obtained in-situ for each deployment [37]. The calibration line was obtained from a linear regression of rotation rate against swim speed that was calculated from vertical depth change divided by sine of the pitch at 5 s intervals when absolute mean sine of pitch was greater than 0.7-0.9. For the DTAG data, speed during glides was estimated using the rate of change of depth divided by the sine of pitch [38].
Data during glides were extracted in 5 s duration segments [15]. Glides shorter than 5 s were excluded from the analysis and glides longer than 5 s were broken into 5 s sub-glides. For each 5 s sub-glide, mean depth (d), speed (v) and pitch angle (p) were calculated. Acceleration (a) was measured by regressing speed versus time over each 5-second interval (S1 Fig). The variance of the acceleration measurement during each 5 s sub-glide was quantified as the root mean square of residuals from the fitted regression line. Seawater density (ρ sw ) for each sub-glide was calculated from a CTD cast that was made close in time and location to each tagged whale. In this analysis, we only used stable glides (circular variance of roll < 0.1) that were at steep pitch angle (absolute pitch > 30˚) to enable robust estimates of speed for DTAG records. In addition, any glides associated with lunge feeding were excluded from the analysis because body form and kinematics of whales drastically change during this feeding behaviour [51]. Lunge feeding events were detected as peaks in jerk (i.e. differential of acceleration) for DTAG records [52]. For 3MPD3GT records with speed data, a lunge was detected as peak in speed when the speed exceeded the threshold of mean speed plus two standard deviations followed by a rapid deceleration. According to a fine-scale kinematic study of lunge-feeding humpback whales, whales stroke throughout lunges but glide at the end of feeding once the mouth has been closed [52]. To exclude any feedingrelated glides, we excluded any glides recorded within 46 s after the lunge from the analysis because it has been reported that humpback whales spend an average of 46 s for filtration and prey handling [52].

Hydrodynamic performance model
We used the equation presented by Miller et al. [15] where acceleration (m s -2 ) along the swimming path is determined by drag force (the first term) and buoyancy forces derived from body tissue (the second term) and gasses carried by each whale (the third term): Here, C D is the drag coefficient, A is the relevant surface area (m 2 ), m is the mass of the whale (kg), ρ sw is the density of the surrounding seawater (kg m -3 ), v is swim speed (m s -1 ), ρ tissue is the density of the non-gas component of the whale body (kg m -3 ), g is acceleration due to gravity (9.8 m s -2 ), p is animal pitch (radians), V air is the volume of air at the surface (m 3 ), ρ air is the density of air (kg m -3 ), d is glide depth (m), and r is compressibility for animal tissue (i.e., the fractional change in volume per unit increase in pressure). The value 101325 converts pressure in atmospheres to pressure in Pascals, so that the units of body tissue compressibility are proportion per Pascal x 10 −9 .
The first additive term of the equation represents the effect of drag on the forward motion of the whale during a glide, which is primarily a function of speed itself. C D Am -1 is the unknown term that is treated as a single quantity in this approach with units of m 2 kg -1 . The second term quantifies the effect of net buoyancy derived from unknown tissue density (ρ sw ) on speed during a glide. The third term quantifies the influence of net buoyancy derived from the unknown volume of gas per unit mass carried in the dive (V air m -1 ) on speed during a glide. As gas compartments of whales are compressed during dives, the volume and density of gas carried by the animal are modelled to change with hydrostatic pressure following Boyle's Law. The model also includes the effect of tissue compressibility (r) that was fixed as 0.38 x 10 −9 Pa -1 based on the value estimated for northern bottlenose whales [15].

Bayesian estimation
The unknown parameters in the hydrodynamic glide model (mainly ρ tissue , V air m -1 and C D Am -1 ) were estimated by Bayesian Gibbs sampling with the freely available software JAGS within R (coda, R package v0.17-1 2015, http://cran.r.project.org/web/packages/coda/index.html) and R2jags (R package v0.5-7 2012, https://cran.r-project.org/web/packages/R2jags/index.html) using data extracted for each 5-s sub-glide. Acceleration during glides was measured using a linear regression line of speed versus time. Observation error measured from variance of acceleration for each 5 s was incorporated in the model by treating acceleration as a normal variable with a precision parameter (1/variance) [15]. A small increment (0.001) was added to the standard errors to ensure finite values for the precision parameter. For the Bayesian estimation, a specific prior distribution must be set for each unknown parameter. A non-informative uniform prior from 800 to 1200 kg m -3 was set for body tissue density (ρ tissue ). An informative prior was set for the combined drag coefficient term (C D Am -1 ) based on several sources of information: drag coefficient (C D ) was estimated to be 0.0026 based on the value estimated for a fin whale (Balaenoptera physalus) swimming at 4 m s -1 [53]. Based on body lengths (L) ranges from 6 to 15 m, body mass (m) was estimated as 20005 kg on average (range 3253-48556 kg) using an equation derived for humpback whales: m = 0.016473L 2.95 x 1000 [54]. Surface area (A) was estimated as 47.4 m 2 (range 15.3-89.0 m 2 ) using a prediction equation obtained from bottlenose dolphins (Tursiops truncatus): A = 0.08m 0.65 [55]. Thus, an expected value for the combined drag term (C D Am -1 ) would be 7 x 10 −6 m 2 kg -1 , with a range from 5 x 10 −6 m 2 kg -1 for large whales to 12 x 10 −6 m 2 kg -1 for small whales. In order to capture uncertainty around this expected value, we specify the prior to be a normal distribution with a mean of 7 x 10 −6 m 2 kg -1 and standard deviation of 2 x 10 −6 m 2 kg -1 that was truncated at 1 x 10 −6 m 2 kg -1 and 20 x 10 −6 m 2 kg -1 . For diving gas volume (V air m -1 ), a uniform prior from 5 to 80 ml kg -1 was set based on the total lung capacity (65-72 ml kg -1 ) estimated for 6 to 15 m long whales using an equation derived from various marine mammals: total lung capacity = 0.10m 0.96 x 1000 [56].
Following Miller et al. [15], we explored variability of unknown tissue density, combined drag term and diving gas volume by evaluating a total of 12 model structures. We fitted a model in which the quantity of the unknown parameters ρ tissue , V air m -1 and C D Am -1 remained constant across the tags and dives (global estimates). We also fitted hierarchical models in which the individual-specific estimates of tissue density and/or drag term, and the dive-specific estimates for diving gas volume are sampled from each global (i.e. individual-average or diveaverage) distribution that was estimated for each parameter. See the JAGS script in the appendix of Miller et al. [15] for the detailed structure of the hierarchical model. All models were sampled in three independent chains, with 24,000 iterations each. The first 12,000 samples were discarded for burn-in, and the remaining posterior samples were downsampled by a factor of 36 to remove any serial correlation in the samples. We report the mean and 95% percentile, hereafter termed posterior mean and credible interval (CI), of the posterior samples as the best estimates of the parameter value and its uncertainty. The 95% credible interval is the Bayesian analogue for the more traditional (frequentist) confidence interval, and defines the range of values within which the true parameter value lies with 95% probability, given the observed data. Convergence was assessed for each parameter, using trace history and Brooks-Gelman-Rubin diagnostic plots [57]. The best model was selected based on the deviance information criterion (DIC), with a lower value indicating a better model fit relative to model complexity.

Results
A total of 33 tag datasets were analysed (Table 1). In the Gulf of St Lawrence, archival tags were deployed on 12 whales in the Jacques-Cartier Passage and adjacent waters between July and September 2011. All tagged whales were part of a long-term photo-identification study that has been carried out at the study site since 1984 [58]. Photographic and field observations of behaviour and known associates suggest that at least two adult females (H002 and H584_2) were pregnant when the tag data were collected. Pregnancy of H002 was also confirmed by hormonal analysis of blow samples and blubber samples. One adult male (H607) was tagged twice at the beginning of the feeding season (July 22, 2011) and later the same season (September 1, 2011). At the Antarctic field-site, 19 whales were tagged over the course of two field seasons that ran between May and June in both 2009 and 2010. Antarctic animals were tagged in Wilhelmina and Andvord Bays along the WAP and inshore waters of the Gerlache Strait. Two pairs of tagged whales were found to be mother-and-calf pairs based on visual observation from the tag boat and biopsy samples ( Table 1).The whales conducted dives to a maximum depth of 388.3 m. Mean swim speed throughout dives was 1.5 ± 0.4 m s -1 (± SD, Table 2). Gliding was observed both during descent and ascent phases although the percentage of time spent gliding varied among whales ranging over 1.5-45.2% and 2.8-60.0% during descent and ascent phases, respectively. Pitch angles during descent and ascent phases were -39.8 ± 20.6˚and 30.6 ± 22.4˚on average, respectively ( Table 2). From the whole dataset, we extracted a total of 18546 5-s sub-glides that were not associated with a lunge. However, 73.7% of these glides were filtered out due to shallow pitch angle (< 30˚) and 0.1% due to high variability in roll (circular variance of roll > 0.9). In addition, 1.4% of glides were removed due to lack of speed and/or acceleration data throughout the 5-s glides. As a result, 24.7% of the total 5-s sub-glides met the criteria for the use of hydrodynamic glide model. The number of 5-s sub-glides that could be used for the hydrodynamic glide model was positively correlated with the duration of tag dataset (Spearman's rho = 0.633, p < 0.001; Fig 1A) although the number of useable glides also varied depending on the behaviour of the tagged whales (foraging, resting, etc.). Eight tag datasets were excluded from the Bayesian estimation of tissue density because of insufficient sample size (<20 sub-glides in each dataset; Table 1). Data of Mn11_H698 was also excluded because an in-situ calibration of the speed sensor was not applicable for this deployment.
Twenty-four of the 33 tag datasets (10 from the Gulf of St Lawrence and 14 from Antarctica) were used to estimate tissue density and the other unknown parameters. Of the 12 Bayesian models, the model with the lowest DIC indicated global plus individual variation in tissue density and drag terms, and global plus dive-by-dive variability in diving lung volume ( Table 3). The difference in DIC from the next-best model was 1657.5 units.
The global body tissue density was estimated with a posterior 95% credible interval (CI) of 1029.5-1033.6 kg m -3 (mean = 1031.6 kg m -3 ). Individual posterior mean values ranged from 1025.2 to 1043.1 with ±95% CI of 0.04-2.3 kg m -3 . The 95% CI range for individual tissue density estimates decreased with increasing number of 5-s sub-glides in the dataset (Fig 1B). There was no significant relationship between the 95% CI range and the average depth at which the sub-glides occurred (Spearman's rank test, p = 0.22); depth of glides ranged from 5.1 to 343.2 m with individual mean ranging 25.2 ± 10.8 to 97.3 ± 55.4 m. There was a tendency for the 95% CI ranges to be smaller for the whales tagged in Antarctica using DTAGs (0.6 ± 0.5 kg m -3 ) than in the Gulf of St Lawrence using 3MPD3GTs (2.5 ± 1.3 kg m -3 ). It is possible that different sampling frequencies and resolution of sensors as well as speed determination methods (measured/ estimated) of 3MPD3GTs and DTAGs might influence the precision of tissue density estimates. Yet, the effect of the two different archival tag models could not be fully addressed due to the differences in location itself and longer data duration for the Antarctic DTAG dataset (22.9 ± 1.6 h) compared to the Gulf of St Lawrence 3MPD3GT dataset (3.5 ± 1.5 h; Table 1 (Table 1). Tissue densities of two pregnant females were estimated as the lowest (1026.5 ± 0.5 kg m -3 for Mn11_H0 02) and the second lowest (1028.6 ± 0.7 kg m -3 for Mn11_H584_2) among the whales from the Gulf of St Lawrence (Table 1). There was a significant negative correlation between relative tissue density to seawater and percent time spent gliding during ascent vs decent phases of non-feeding dives (Spearman's rho = -0.72, p <0.001; Fig 3).
The posterior mean of the global drag term was 11.8 x 10 −6 ± 1.6 x 10 −6 m 2 kg -1 (± 95% CI). The posterior mean was higher and the distribution had little overlap with the prior distribution that had a mean of 7.0 x 10 −6 m 2 kg -1 (Fig 4). The posterior means of individual drag term values ranged from 6.0 x 10 −6 to 25.5 x 10 −6 m 2 kg -1 , but most of them were near 12.5 x 10 −6 m 2 kg -1 ( Table 1). The posterior mean of global diving gas volume was 27.7 ± 1.1 ml kg -1 (±95% CI).

Discussion
To date, the hydrodynamic glide model has been used to estimate tissue density of deep diving marine mammals such as elephant seals, sperm whales, northern bottlenose whales and longfinned pilot whales [15,34,38,47]. In this study, we successfully applied this method to substantially shallower-diving humpback whales to estimate tissue density from two geographically distinct feeding populations. To examine the variability of the unknown parameters (tissue density, drag and diving gas volume), we fitted 12 models with different model structures. The best model included individual variation in tissue density and drag, supporting our expectation that each whale had different tissue density. The best-fitting model also included dive-by-dive variation in diving gas volume. Although there was no apparent overall relationship between diving gas volume and dive duration, it is possible that whales change the amount of inhaled air before dives depending on their activity [15]. The gliding patterns of whales correlated with their estimated tissue density, with denser whales spending relatively more time gliding during descent and less-dense whales spending more timing gliding during ascent phases (Fig 3). The significant correlation of tissue density and gliding patterns provides a degree of validation that the tissue density estimates, or at least their relative values, were accurate.

Drag term estimates
The drag coefficient is one of the key parameters to estimate tissue density using the hydrodynamic glide model. Following Miller et al. [15], the combined drag term (C D Am -1 ) was estimated using a relatively narrow Gaussian prior that was determined based on auxiliary published data in order to improve the precision of tissue density estimates. However, the global (individual-average) estimate of the drag term in the best-fitting model (11.8 x 10 −6 m 2 kg -1 ) did not concentrate within the distribution of the prior (7.0 x 10 −6 m 2 kg -1 ).
As in previous studies [15,34,38], we neglected any specific effect of lift, although liftrelated induced drag may not be negligible in the case of humpback whales due to their large pectoral flippers [48] and propensity to glide at shallow angles. It is possible that the influence of induced drag due to lift generation may explain the mismatch of the prior expectation of the combined drag term and its posterior estimate from the data. Adding the induced drag to the hydrodynamic glide model, the drag part of the equation can be expressed as Model structure refers to the allowed variation in the model for the unknown terms, with G referring to global (i.e. individual-average) parameter only, I referring to individual specific estimates included, and D referring to dive-by-dive variation included. The column head refer to ρ tissue as tissue density (kg m -3 ); C D Am -1 as combined drag term (m 2 kg -1 ); V air as volume of air (ml kg -1 ). Data are presented with ±95% CI in parentheses. For global parameter estimates, .g refers to the global parameter and .var refers to individual or dive-by-dive variance. https://doi.org/10.1371/journal.pone.0200287.t003 where A Flipper is flipper surface area (m 2 ), AR is flipper aspect ratio and C L is the lift coefficient [59]. Because both the parasite drag and the induced drag are a function of speed-squared, the equation can be rewritten as Thus, the structure of the equation is unchanged just with the addition of induced drag to that of the parasite drag term C D Am -1 .
We suggest that the model estimated higher global C D Am -1 values due to the effect of induced drag by assuming that the model estimated the combined term in parenthesis, instead of the parasite drag term (C D Am -1 ) alone. The lift coefficient of a humpback whale flipper is estimated as 0-0.9 through wind tunnel measurements [60]. Based upon literature values for the surface area (A Flipper , 12.20 m 2 ) and the aspect ratio (AR, 5.67) of a humpback whale flipper [48], A Flipper C L 2 /(πARm) is estimated as 0-22 x 10 −5 for a 12-m long whale. Adding this value to 7 x 10 −6 (i.e. mean of the C D Am -1 prior), the combined drag term in the parenthesis is expected to range between 7 x 10 −6 and 29 x 10 −6 m 2 kg -1 which overlaps with the global drag term estimates in this study (11.8 x 10 −6 ± 1.6 x 10 −6 m 2 kg -1 , ± 95% CI). This suggests that the mismatch between the prior and C D Am -1 estimates derived from the addition of the induced drag and that lift-related drag forces should not be ignored for this species. Relationship between gliding patterns and relative tissue density. The y-axis indicates differences in the percentage of time spent gliding during ascent and descent phases of non-feeding dives by each whale. Vertical and horizontal error bars show standard deviation and 95% credible interval range, respectively. A relative tissue density of >1 indicates that tissue density was denser than surrounding seawater. https://doi.org/10.1371/journal.pone.0200287.g003 In comparison with deeper diving marine mammals in previous studies [15,34], the effect of lift seems particularly important for humpback whales that glide at shallower pitch angles ( Table 2) where lift generation increases with correspondingly greater induced drag. Humpback whales have large flipper with a high aspect ratio that can produce lift forces to support their acrobatic movements such as high-speed turning and banking that are associated with feeding [48,61]. In addition, the scalloped leading edge of their large flippers serves to delay stall angles and increase lift [60,62]. Recent studies using animal-borne video camera reported that humpback whales also perform lift-generating flipper strokes for propulsion during lunge feeding [63]. In our study, we only used data during stable glides (circular variance of roll < 0.1) to minimize the influence of lift during maneuvering. However, the influence of lift-induced drag is detectable in our dataset possibly because humpback whales likely use their wing-like flippers to produce lift during stable glides at non-vertical pitch angles. Yet, it is noteworthy that our general results about tissue density seem to be robust because the model quantified the combined effect of parasite and induced drag. As a sensitivity analysis, we refitted the model using a non-informative wide range prior for the drag term instead of a narrow Gaussian prior. The resulting global average drag term was 13.1 x 10 −6 ± 2.4 x 10 −6 m 2 kg -1 and global average tissue density was 1031.6 ± 2.1 kg m -3 , which differed very little from the estimated values with a narrow prior. Similarly, individual tissue density estimates were nearly identical to the result of our original model, supporting the robustness of the tissue density estimates to the prior specification. Thus, the general results about tissue density seem to be robust because the model appears to have estimated a reasonable value for the combined effect of parasite and induced drag.

Body tissue density
Estimated individual-average (global) body tissue density of humpback whales (1031.6 ± 2.1 kg m -3 ; Table 3) was similar to that of other cetaceans reported to date (1030.0 ± 0.8 kg m -3 for Physeter macrocephalus [38]; 1031.5 ± 1.0 kg m -3 for Hyperoodon ampullatus [15]), indicating that non-gas body tissues are typically denser than seawater. However, long-finned pilot whales were estimated to have even denser tissues of 1038.8 ± 1.60 kg m -3 [47]. For humpback whales in this study, a large variation was detected in individual-specific body tissue density ranging from 1025.2 to 1043.1 kg m -3 , as we expected, because individual tissue density at feeding grounds would change depending on factors such as age, sex, reproductive status, prey availability and the number of days since arrival at the feeding ground [4,6,18,24]. In a study of fin whales conducted using Icelandic whaling data, pregnant females had the highest rate of fattening during the feeding season as they increased their total body energy content by nearly 80% [4]. A similar trend was reported for minke whales (Balaenoptera acutorostrata) in Iceland: the blubber volume of pregnant females almost doubled over the feeding season [24]. Using the hydrodynamic glide model, high lipid-stores of two pregnant female humpback whales (Mn11_H002 and Mn11_H584_2) were indicated by low tissue density estimates of 1026.5 kg m -3 and 1028.6 kg m -3 that were the lowest and the second lowest, respectively, among all of the tagged whales in the Gulf of St Lawrence. A decrease in tissue density over the feeding season due to accumulation of lipid stores was also detected in this study: tissue density of a repeated sampled adult male (Mn11_H607) decreased from 1037.0 to 1031.2 kg m -3 in 40 days. Based on extrapolation from elephant seals, the proportion of lipid content (P lipid ) corresponding to these tissue densities of Mn11_H607 would be 36.3% and 39.0%, as determined from ρ tissue = ρ lipid P lipid +ρ lipid-free (1-P lipid ), where ρ lipid and ρ lipid-free are 900.7 and 1114.6 kg m -3 , respectively [34]. Our results also showed that one of the two calves had low body tissue density of 1025.2 kg m -3 (Mn10_155b) in agreement with general expectation that calves are more buoyant because they deposit fat during the lactation period [64]. The other calf (Mn10_139a), however, had relatively high tissue density of 1040.8 kg m -3 that was supported by its gliding pattern suggestive of negative buoyancy: the whale spent more time gliding during descent (61.3%) than ascent (45.7%) phases of non-feeding dives (Fig 5D). It is possible that Mn10_139a had a poor body condition, reflecting its mother's poor condition indicated by its relatively high body density (Mn10_139b, 1029.4 kg m -3 ) compared to the other mother in the study (Mn10_155a, 1027.6 kg m -3 ).
During the feeding season, it is essential for humpback whales to accumulate a sufficient amount of energy for survival, growth and/or reproduction. Previous studies estimated the amount of energy gained by baleen whales over the course of a feeding season via anatomical measurements and chemical analysis of multiple whale carcasses [4,21]. More recent work has described dynamic foraging patterns of whales throughout the course of the foraging season suggesting that whales alter their feeding behaviour (rates and dive depth) commensurate with changes in the availability of prey [65]. This could lead to non-linear changes in the accumulation of energy, and, combined with body density estimates collected over similar time periods, can offer insights as to the most critical times and locations for whales to regain energy stores and how different life history classes vary. This information is critical to understanding how environmental changes and potential human disturbance can significantly impact individual and population-level health of marine mammals and other animals.
Changes in tissue density lead to changes in buoyancy that influence swimming patterns of diving animals given strong selection for them to travel efficiently to and from depth [33,39]. For example, it is expected that animals with higher density should glide more during descent aided by negative buoyancy whereas less dense positively buoyant animals should employ more glides during ascent. In agreement with the expectation, a negative correlation between tissue density estimates and percent time spent gliding during ascent vs. descent phases of non-feeding dives was observed (Fig 3), suggesting that the model successfully detected relative differences in individual tissue density. Because there is greater variability in tissue density of humpback whales than deep diving toothed whales [15,38,54], the relatively low precision of tissue density estimates obtained here (±95% CI of individual tissue density ranged up to 2.3 in this study whereas to 0.4 in [15]) seem to be sufficient to detect individual and/or temporal variation. Lower precision tissue density estimates may be expected if there is high variability in acceleration that is not accounted for in the model. Therefore, sample size is particularly important to consider for the body density estimation of humpback whales for which induced drag may cause variability in gliding acceleration. The result showed that the number of 5-s sub-glides in each dataset is one of the key factors affecting the range of 95% CI for the posterior estimates of individual tissue density. Specifically for humpback whales, >200 sub-glides in each dataset seem to be needed to obtain highly precise estimates with 95% CI range of 1 kg m -3 (Fig 1B).
In this study, we estimated tissue density of humpback whales at two geographically distinct feeding grounds. Tissue density of whales from Antarctica and Gulf of St Lawrence largely overlapped, but there was a tendency for Antarctic whales to have lower tissue density (Fig 2), indicating that animals in that location at that time had larger lipid reserves than did the animals tagged in Canada. It is possible that the geographic differences reflected different temperature and prey conditions on two feeding grounds. However, as numerous factors can affect individual tissue density, more data, including basic information of individuals such as sex, age class and reproductive status that can be obtained from photo-ID and biopsy studies would be essential to identify factors that cause these geographic differences.
The methods used in this study closely followed methods published for other deep-diving odontocete species. Some studies used an estimated value of 0.06 for the entrained mass of water m e which is moved forward along with the body of the animal [38]. Because we did not have specific measurements of animal mass in this study, we also did not include estimates for entrained mass. We also do not expect addition of a constant mass proportion of mass to all estimates would affect their relative values. However, the absolute values obtained in this study could be made more accurate with finer estimates of the length, mass and surface area of each whale as was done using photogrammetry by Miller et al. [38]. We do recommend incorporation of such data when available, and to include estimates of entrained mass to obtain accurate absolute values of body tissue density.

Diving gas volume
The posterior mean of global (dive-average) diving gas volume was 27.7 ± 1.1 ml kg -1 (±95% CI). This value is substantially lower than the estimated lung volumes of mysticete fin (29-61 ml kg -1 ) and sei (Balaenoptera borealis) whales (116-151 ml kg -1 ), whose lung volumes were measured via inflation of excised lungs [66,67]. This could indicate that: 1) lung volumes of mysticete whales are smaller than estimated from excised lungs in which some amount of air is likely to be trapped [46]. Piscitelli et al. [67] noted that the mass specific volume of sei whales from that study were outliers on a comparative basis relative to smaller cetaceans; 2) humpback whales in this study dove with less than their full lung capacity; or 3) our estimate was incorrect and too low.
As cetaceans appear to inhale immediately prior to diving, the diving volume of cetaceans is thought to be close to the total lung volume [68]. In fact, the calculated diving lung volume of deep-diving northern bottlenose whales (27 ml kg -1 , [15]) was similar to the measured total lung volume of 28 ml kg -1 ([69] reviewed in [68]). However, shallower diving species may not always dive with full lungs; for example, the diving lung volume and the total lung volume of bottlenose dolphins are 40-50 ml kg -1 [70] and 50-91 ml kg -1 , respectively (reviewed by [67]). Differences in lung sizes and thoracic morphology of shallow and deep diving cetaceans have been reported [46]. As the effect of air-derived buoyancy is stronger at shallower depth, it is possible that shallower-diving whales do not always dive with full lung capacity. A large variation in dive-by-dive estimates of diving gas volume found in this study and, albeit weak, the positive relationship between diving gas volume and dive depth would support this hypothesis. While no systematic variation of diving gas volume in relation to dive duration was detected, further detailed analysis of dive-by-dive variation in diving gas volume could provide new insights into their diving physiology.
Another possible explanation for the low estimate of diving gas volume is that the amount of gas stored in the body might decrease during dives. It has been reported that humpback whales actively exhale underwater in some situations. For example, humpback whales have a diverse repertoire of feeding behaviours, including "bubble feeding" that involves underwater exhalation to form bubble clouds, nets or curtains to corral prey [61,71]. Bubbling is also observed in non-feeding situations such as play [71] and social interactions [72]. Although apparent bubbling was not detected from acoustic audits of the DTAG datasets in this study, it is possible that some air might passively escape from the body during dives. If such underwater exhalation and/or passive loss of air occurred, our estimate of diving gas volume would be too low because the majority of the glides used in the analysis were recorded during ascent phases of dives, and thereby the estimate reflects the amount of gas in the body at latter part of dives.

Conclusion and future directions
We demonstrated that the hydrodynamic glide model can be used to detect individual and temporal variation in body tissue density of humpback whales, suggesting that it is likely to be broadly applicable across a range of aquatic animals including shallow diving baleen whales. The important next step is validation with other techniques such as visual assessment [28,29], biopsy sample measurements, and photogrammetric measurements of body width versus length using overhead images [11,30].
This study represents a cross-sectional design, in which the tissue densities of multiple animals were measured. Longitudinal tracking of changes in individuals' tissue density as has been done with elephant seals drift dives [5], or repeated measurements as we made for whale Mn11_H607, may be a more powerful approach to determine specific factors that affect the lipid-store body condition of humpback whales. Considering that humpback whales are less difficult to tag such that multiple tagging of the same individual is possible, this tag-based minimally invasive approach may provide an effective tool to monitor body tissue density as a measure of body condition. By integrating life-history data of individuals (e.g. age, sex, size, reproductive status) as well as prey availability at feeding grounds, this approach can be helpful to understand bioenergetics and health of individual whales within increasingly human-altered ecosystems. Ultimately, tracking the tissue density of individual whales using longer duration tags could be a powerful technique to relate their body condition to how they interact with features of their natural environment.