Diel patterns in swimming behavior of a vertically migrating deepwater shark, the bluntnose sixgill (Hexanchus griseus)

Summary

Diel vertical migration is a widespread behavioral phenomenon where organisms migrate through the water column and may modify behavior relative to changing environmental conditions based on physiological tolerances. Here, we combined a novel suite of biologging technologies to examine the thermal physiology (intramuscular temperature), fine-scale swimming behavior and activity (overall dynamic body acceleration as a proxy for energy expenditure) of bluntnose sixgill sharks (Hexanchus griseus) in response to environmental changes (depth, water temperature, dissolved oxygen) experienced during diel vertical migrations. In the subtropical waters off Hawai‘i, sixgill sharks undertook pronounced diel vertical migrations and spent considerable amounts of time in cold (5–7°C), low oxygen conditions (10–25% saturation) during their deeper daytime distribution. Further, sixgill sharks spent the majority of their deeper daytime distribution with intramuscular temperatures warmer than ambient water temperatures, thereby providing them with a significant thermal advantage over non-vertically migrating and smaller-sized prey. Sixgill sharks exhibited relatively high rates of activity during both shallow (night) and deep (day) phases and contrary to our predictions, did not reduce activity levels during their deeper daytime distribution while experiencing low temperature and dissolved oxygen levels. This demonstrates an ability to tolerate the low oxygen conditions occurring within the local oxygen minimum zone. The novel combination of biologging technologies used here enabled innovative in situ deep-sea natural experiments and provided significant insight into the behavioral and physiological ecology of an ecologically important deepwater species.

Keywords:

Hidden Markov models – Sharks – Oxygen – Body temperature – Animal migration – Swimming – Oxygen metabolism – Dissolved oxygen

Introduction

Diel vertical migration, where organisms migrate through the water column for feeding, predator/competitor avoidance or to maximize energetic efficiency, is a widespread phenomenon in aquatic ecosystems [1–5]. Based on physiological tolerances, organisms may alter their behavior in response to changing environmental conditions encountered during migration [6–9]. In the subtropical waters off Hawai‘i, bluntnose sixgill sharks (Hexanchus griseus; hereafter referred to as sixgill sharks) undergo a distinct diel vertical migration descending to depths below 500 m during the day and ascending to 200–350 m during the night and encounter a wide range of temperatures (5–20°C [10–11]) and dissolved oxygen conditions (9–85% saturation [12]; <60–200 μmol kg^-1 [10]). Sixgill sharks appear to be highly tolerant of low dissolved oxygen, remaining in waters that are generally considered hypoxic (defined as <61 μmol kg^-1 [13–15]) within the local oxygen minimum zone (OMZ) throughout their deeper, daytime distribution [10,12,16].

Diel vertical migrators may tolerate hypoxic conditions during daytime forays into OMZs via: (1) enhanced oxygen extraction and transport in the gills, blood and tissues; (2) upregulation of anaerobic metabolism; and (3) suppression of total energy consumption by reducing energetically-expensive processes such as locomotor activity [17,18]. Further, these are not mutually exclusive, and many organisms employ a combination of these strategies to accommodate local oxygen conditions and expand vertical habitat [17,18]. Therefore, sixgill sharks may restrict locomotor activity as a strategy to reduce overall metabolic expenditure in order to facilitate brief (hours) forays into cold, hypoxic water [14,17,18].

Ectothermic species, such as sixgill sharks, must also contend with changes in environmental temperature, which has a strong influence on physiology and fitness from individual to community scales [19–22]. Sixgill sharks off Hawai’i spend their shallower nighttime distribution (200–350 m) at variable temperatures (11–15°C) within the thermocline and the majority of their daytime distribution occupying stable, cold temperatures (5–7°C) at depths below 500 m [10,11]. As such, sixgill sharks may modify their behavior to enable them to take advantage of diverse thermal habitats [23,24]. For example, the large body mass of sixgill sharks may provide sufficient thermal inertia to slow heat loss during their descent into cold, deep waters [25,26]. Thermal inertia may enable sixgill sharks to maintain relatively high activity levels for longer periods at colder temperatures, and thus exploit deep prey resources more effectively compared to resident (non-vertically migrating) species. Sixgill sharks could also migrate between foraging and non-foraging habitats to select optimal temperatures for bioenergetic and/or digestive efficiency [23,27,28].

Recent advances in biologging have enabled the recording of in situ dissolved oxygen contemporaneously with the behavior of the tagged animal [12]. Here, we combined this novel technology with an accelerometer, temperature-depth recorder and intramuscular thermistor to examine for the first time the fine-scale swimming behavior, activity, and thermal physiology of sixgill sharks in response to environment changes experienced during diel vertical migrations. Specifically, we hypothesize that sixgill sharks reduce activity levels during their deeper daytime distribution in response to low temperature and dissolved oxygen conditions.

Materials and methods

Data collection

Five sixgill sharks (2 male, 3 female; Table 1) were captured using demersal longlines set off Kāne‘ohe Bay, O‘ahu, Hawai‘i (21.46°N, 157.80°W) (for additional details of longline fishing methods see [11]). Both males (326 and 352 cm total length [TL]) were mature and the three females (309, 346, and 412 cm TL) were large subadults based on size at maturity estimates [29]. To quantify swimming behavior, we used a TDR10-X accelerometer (57 × 38 × 24 mm, 69 g; Wildlife Computers, Redmond, WA, USA), which recorded tri-axial acceleration at 1/32-s intervals. Environmental data were obtained using a prototype dissolved oxygen pop-up satellite archival tag (DO-PAT: 170 × ⌀ 60 mm, 85 g; Wildlife Computers, Redmond, WA, USA; for additional details see [12]) or a prototype TDR10-DO inline finmount tag (109 × 53 × 21 mm, 118g; Wildlife Computers, Redmond, WA, USA) equipped with a micro dissolved oxygen probe (27 × ⌀ 12 mm; OxyGuard International A/S, Farum, Denmark). In total, this combination of sensors recorded depth, water temperature, and dissolved oxygen saturation at 1-s intervals. Intramuscular temperature was recorded at 10-s intervals using a TDR-Mk9 archival tag (71 × 17 × 17 mm, 120 mm stalk length, 29 g; Wildlife Computers, Redmond, WA, USA) with a temperature sensor incorporated on a sensor stalk that was inserted approximately 8 cm below the skin into the dorsal musculature. These devices were embedded in small, syntactic-foam floats equipped with an Argos satellite transmitter and UHF radio tracking beacon (MiniPAT: 115 × ⌀ 40 mm, 53 g; or SPOT-216E: 80 × 20 × 11 mm, 30 g; Wildlife Computers, Redmond, WA, USA), and a timed-release mechanism (19 × ⌀ 16 mm, 10 g; Little Leonardo Co., Tokyo, Japan). The floats provided just enough overall buoyancy to return the instrument package to the surface on release. Total masses of the devices (400 and 528 g in air) were less than 0.3% of the shark’s estimated mass (Table 1). To attach the package a hole was pierced on the central dorsal surface through the skin with a sharp probe and the device was secured parallel to the longitudinal axis of the body by a stainless-steel cable tie (360 mm, 8 g; Little Leonardo Co., Tokyo, Japan) passed through the hole. After ~96 h, a pre-programmed release timer severed the cable tie and detached the package; however, the TDR-Mk9 sensor stalk took an additional 0.5–55 h to pull out of the dorsal musculature before allowing the package to float to the surface where the satellite and UHF transmitters enabled us to locate and retrieve the devices. Procedures were approved by the ethics committee at the University of Hawai‘i (Institutional Animal Care and Use Committee Protocol #05–053).

Data processing

The first 15–40 h of each record were omitted from analysis because the behavior of the shark was inconsistent with the remainder of the record and indicated a recovery period after capture and handling. Analyses of acceleration and intramuscular temperature data concluded ~96 h post-deployment when the pre-programmed release timer severed the cable tie and detached the instrument package from a fixed position, which was evident in the acceleration data. The final 31 h of environmental data recorded by the DO-PAT from HG3 were corrupted, so only the first 120 h were analyzed. The accelerometer battery died after 10 h for HG4; therefore, only environmental and intramuscular temperature data were analyzed for this shark. After removing these data from the analysis, there was a total of 265 h (range 55.6–80.4 h) of accelerometer data from 4 sharks and 321.5 h (range 56.0–80.7 h) of intramuscular temperature and 402 h (range 56.1–105.0 h) of environmental data (depth, water temperature, and dissolved oxygen saturation) from 5 sharks.

Dissolved oxygen saturation data were smoothed relative to corresponding depth data using a robust local regression (LOESS) algorithm with a 5-s window in MATLAB (The MathWorks Inc., Natick, MA, USA) (sensu [12]). The smoothing algorithm effectively reduced the noise from the last deployment (HG5), but some outliers persisted for the remaining deployments. These outliers were identified and removed using studentized residuals and new values were based on nearest neighbor interpolation.

Data analysis

Diel patterns

To describe overall diel patterns of habitat use, the proportion of time spent at depth (5 m bins), water temperature (0.5°C bins), and dissolved oxygen saturation (1% bins) by time of day was calculated at 15 min intervals for every complete 24 h period for each shark and then averaged across all sharks. Vertical displacement rates were calculated by dividing the change in depth between each record by the sampling rate. The majority (98%) of vertical displacement rates calculated using the high-resolution depth time series (recorded at 1-s intervals) ranged from 0–1 m s^-1 (maximum 2.5 m s^-1) with a depth sensor resolution of 0.5 m. Therefore, data were sub-sampled at 1-min intervals to better describe diel patterns in vertical displacement rates. The depth time series for each individual was visually examined and divided into four categories (sensu [31]): shallow night phase, deep day phase, and two crepuscular vertical migrations (dawn descent and dusk ascent). To describe patterns of habitat use for shallow night phases and deep day phases, the proportion of time spent at depth, water temperature, and dissolved oxygen saturation was calculated using the high-resolution time series for each shark and then averaged across all sharks.

Intramuscular temperature

To describe changes in intramuscular temperature across diel vertical migrations, we estimated changes in the whole-body heat-transfer coefficient k that are necessary to account for the observed rates of body warming and cooling using a function of heat exchange with the environment and internal (metabolic) heat production [7,32]. Heat loss (or gain) is proportional to the difference between the intramuscular temperature of the shark and ambient water temperatures:

• k = a constant value;

• k={k1;Ta(t)<Tb(t)k2;Ta(t)≥Tb(t),

where k₁ (cooling) and k₂ (heating) are two values for the heat-transfer coefficient. The optimized parameters for each model were estimated based on minimization of the sum of squared errors (SSE).

Diel patterns in the proportion of time spent at intramuscular temperature (0.5°C bins) by time of day and for shallow night and deep day phases were calculated as described in the previous section. We measured the difference in intramuscular temperature from ambient water temperature as ΔT and calculated the proportion of time spent at ΔT (0.5°C bins) during shallow night and deep day phases for each shark and then averaged across all sharks.

Acceleration data

Acceleration data were analyzed using Igor Pro 6.3 (WaveMetrics Inc., Portland, OR, USA) with the ‘Ethographer’ package [33]. Raw acceleration data include both low-frequency gravity components (caused by the shark’s changing body inclination), and high-frequency specific components caused by dynamic movements such as tail beating. We separated specific and gravity components of acceleration using a 0.01 Hz low-pass filter, which removes high frequency waves to reveal the gravity component. We used the gravity component of acceleration to calculate the shark’s body inclination (pitch and roll; sensu [34]). Differences in the pitch angle between the longitudinal axis of the shark’s body and accelerometer were corrected using the relationship with vertical displacement rate as described in [35]. Acceleration data were categorized into three swimming phases (ascents, descents, level) based on pitch angles following [36], assuming pitch angles above 5° to be ascents, below -5° to be descents, and intermediate angles to be level swimming. We used continuous wavelet transformation to generate a spectrogram of swaying acceleration, classified the dominant peak as tail beat cycles, and calculated tail beat frequency and amplitude of acceleration at every second. Overall dynamic body acceleration (ODBA) was calculated by summing the absolute specific components of acceleration across all three orthogonal axes. Average (per minute) metrics of pitch and acceleration (tail beat frequency, amplitude of acceleration, and ODBA) during diel depth phases were compared for each individual.

We used generalized additive mixed models (GAMMs) to determine which temporal, physiological, and environmental factors influence ODBA of sixgill sharks. ODBA was used as a measure of activity (as opposed to tail beat frequency) since it measures change in behavior in all axes and has been used as a proxy for energy expenditure across taxa, including sharks [37–40]. We incorporated time of day, depth, water temperature, intramuscular temperature, dissolved oxygen saturation, and swimming phase (ascent, descent, level) as candidate predictor variables. For GAMM analyses, we sub-sampled data over 1-min means. Prior to model fitting, data exploration was carried out per [41]. Collinearity between candidate predictor variables was assessed with Pearson correlation coefficients and variance inflation factors (VIF) using the ‘corvif ‘ function [42] in R [43]. High absolute Pearson correlation coefficients (>0.73) and VIF (>3) indicated high collinearity between all predictor variables except for time of day and swimming phase. Therefore, a principal component analysis (PCA) was conducted on the correlated predictor variables in R. Principal components with greater than 10% of variance explained were included with time of day and swimming phase as candidate predictor variables in GAMM analyses. Time of day and swimming phase were also included with each of the correlated predictor variables in separate models to determine individual environmental effects on ODBA. Interactions between smooth terms were examined using tensor product smooths (t2; [44]). As observations were repeated measures collected from the same individuals, we modeled individual sharks as a random effect. We included a correlation structure with an autoregressive process of order 1 (AR1) to account for serial correlation in time series data [42]. ODBA was modeled using a Gaussian distribution with a log link function. Thin plate regression splines were estimated for each physiological and environmental variable whereas time of day was modeled using a cyclic cubic regression spline, which constrains the start and end points of the smooth term to be the same [45]. GAMMs were constructed in R using the ‘mgcv’ package [45] and model fits were compared using Akaike Information Criterion corrected for small sample size (AIC_c) and Akaike weights (w) [46], calculation of percent deviance explained, and residual diagnostic plots. Models with substantial support were selected based on a ΔAIC_c < 2 from the model with the lowest AIC_c.

We used Hidden Markov models (HMMs) to distinguish between different activity states based on ODBA. HMMs connect observed behavior data, such as ODBA, to an underlying latent process, generally interpreted as the animal’s unobserved behavior [47,48]. HMMs directly account for the serial dependence in our dataset [48,49] and allow modeling with a higher sampling resolution (i.e., seconds instead of minutes). In particular, HMMs provide a measure of the probability of sharks switching or persisting within behavioral states [48] and how these are driven by temporal, physiological, and/or environmental conditions. We applied a two-state HMM with one state post hoc interpreted as representing less active behavior and the other relatively more active behavior following [48] and [50]. For HMM analyses, we sub-sampled data over 10-s means corresponding to the longest sampling rate (10 s) of intramuscular temperature data. Mean ODBA values calculated over a period of 10 s capture 1–2 dominant tail beat cycles for sixgill sharks [11], thereby reducing the variability in activity estimates due to peaks and troughs in acceleration data during the tail beat cycle [51,52].

The state-dependent probability distributions were estimated using a gamma distribution. We considered several sets of initial values for the parameters of the state-dependent distributions and as a result are confident that we found the global optimum of the likelihood function. To examine potential temporal, physiological, and environmental effects on activity states, we allowed the entries of the transition probability matrix to be functions of up to three covariates including swimming phase, time of day, and either a correlated predictor variable or principal component (as described above) and their interaction. Time of day (t) is represented by two trigonometric functions with period 24 h, cos(2πt/8640) and sin(2πt/8640) (24-h periodicity is represented by 8640 time points since data were collected every 10 s). HMMs were constructed in R using the ‘momentuHMM’ package [53]. The ‘momentuHMM’ package currently does not have the capacity to include individual-level random effects; however, we performed a comparable analysis by including individual-level fixed effects on the state transition probabilities [53]. We calculated the HMM likelihood using the forward algorithm, which allows for parameter estimation via a numerical maximum likelihood approach [54]. Model fits were compared using AIC, Akaike weights, and pseudo-residuals [54]. Models with substantial support were selected based on a ΔAIC < 2 from the model with the lowest AIC. We used the Viterbi algorithm to decode the optimal state sequence underlying each time series, thus connecting each observation to one of the two states [54]. To determine how activity varied by depth, water temperature and dissolved oxygen saturation, we related the decoded state sequences to a grid of depth (5 m bins) and water temperature (0.5°C bins) cells and a grid of depth and dissolved oxygen saturation (5% bins) cells. For each grid cell, we calculated the percentage of decoded states corresponding to state 1 or state 2 following [48] and [50]. We estimated the stationary (equilibrium) state probabilities (i.e., the marginal probability of a state at a fixed value of the covariate; see [49]) as a function of each covariate included in the best-fit model. As sixgill sharks exhibit a distinct diel shift in habitat use, we set each covariate to their overall median value for deep day and shallow night phases.

Results

Diel behavior

All sharks undertook diel vertical migrations and spent the majority of their deep daytime distribution (median 11.1 h, IQR 10.5–11.4 h) between 500–650 m (80%) at relatively stable ambient temperatures between 5–7°C (85%) and dissolved oxygen saturations between 10–25% (85%), with a peak between 15–20% (46%) (Fig 1 and S1 Fig). HG5 was the only shark to descend below 700 m and encountered the lowest temperatures (4.7°C) and dissolved oxygen saturations (10%; Table 1); however, only 7% of its overall vertical distribution was below 700 m. During the night, all sharks ascended to shallower waters and spent the majority of their time between 200–350 m (78%), with a peak depth use between 250–300 m (43%) and regularly encountered more variable ambient temperatures between 10–16°C (79%) and dissolved oxygen saturations between 65–80% (55%) (Fig 1 and S1 Fig).

Sixgill sharks underwent extensive changes in depth (descent: median 297 m, IQR 236–304 m; ascent: median 262 m, IQR 204–284 m) alongside considerable changes in ambient temperature (descent: median 7.1°C, IQR 4.4–7.8°C; ascent: median 5.8°C, IQR 5.1–8.5°C) and dissolved oxygen saturation (descent: median 45%, IQR 32–56%; ascent: median 41%, IQR 31–65%). There was no significant difference (Wilcoxon rank sum test; z = 1.439, P = 0.1522) in the duration of dawn descents (median 59 min, IQR 47–80 min) and dusk ascents (median 81 min, IQR 62–100 min). Vertical displacement rates were significantly higher during crepuscular vertical migrations (dawn descents and dusk ascents) compared to shallow night and deep day phases (S1 Table). Overall vertical displacement rates were centered around zero during deep day phases indicating that the sharks were less vertically dynamic compared to shallow night phases where they exhibited a larger range of vertical displacement rates (Fig 2).

Intramuscular temperature analyses

Sixgill sharks demonstrated a thermal inertia that slowed their rates of heating and cooling across the large changes in ambient water temperature experienced during diel vertical migrations. This thermal inertia is probably attributable to their large body size. Overall, sixgill sharks underwent considerable changes in intramuscular temperature (5.8–14.9°C), but these were slightly narrower than the range of change of ambient water temperatures (4.7–17.7°C). There was no significant difference (Wilcoxon rank sum test; z = 0.3723, P = 0.7097) in absolute ΔT between the beginning of dawn descents (median 0.4°C, IQR 0.2–1.1°C) and dusk ascents (median 0.4°C, IQR 0.2–0.7°C). Similarly, there was no significant difference (Wilcoxon rank sum test; z = -0.5584, P = 0.5765) in absolute ΔT between the end of dawn descents (median 4.5°C, IQR 3.5–6.1°C) and dusk ascents (median 5.2°C, 3.3–7.0°C). There was also no significant difference (Wilcoxon rank sum test; z = 0.4964, P = 0.6916) in the overall absolute change in intramuscular temperature during dawn descents (median 0.7°C, IQR 0.3–2.1°C) and dusk ascents (median 0.8°C, IQR 0.5–1.2°C).

Variable k models (2) had lower SSE values compared to constant k models (1) for all sharks except HG5 (Table 2). However, intramuscular temperatures estimated by models (1) and (2) were very similar to observed body temperatures (Fig 3) and the confidence intervals for optimized values of k₁ produced by the variable k model (2) were relatively large (Table 2), suggesting that these models overfit the data. The output from models (1) and (2) were not significantly different (t-test; P > 0.05 for all individuals), and a linear regression of observed intramuscular temperature versus estimated values from constant k models (1) explained ≥99% of the variance in intramuscular temperature (slopes 0.99–1.03, y-intercepts -0.21 to 0.10) for all individuals. Therefore, the constant k model (1) was selected as the most parsimonious explanation. A comparison of optimized constant k values derived from model (1) to estimated body mass indicated that larger individuals had lower whole-body heat-transfer coefficients (Fig 4).

Sixgill shark intramuscular temperature was >0.5°C above ambient water temperature for the majority (69%) of their deeper daytime distribution (Fig 5). In contrast, intramuscular temperature had a difference of -1.5 to 1.5°C from ambient water temperature for the majority (63%) of their shallower nighttime distribution. HG3 and HG4 experienced the warmest ambient water temperatures of 17.3 and 17.7°C, respectively. However, they spent very little time (2% and 1%) at water temperatures above 16°C. Excursions into warm waters were very brief for HG3 (median 5 min, IQR 3–9 min, maximum 41 min) and HG4 (median 2 min, IQR 1–3 min, maximum 11 min) where median intramuscular temperatures were 3.4°C (IQR 2.2–7.1°C, minimum 1.2°C) and 7.8°C (IQR 5.2–8.6°C, minimum 3.7°C) cooler than ambient water temperatures, respectively.

Accelerometer data analyses

Acceleration data indicated sixgill sharks (excluding HG4) exhibited relatively shallow pitch angles during dawn descent (-1.5° to -11.8°) and dusk ascent (5.0° to 8.6°) segments of diel vertical migrations (S2 Table). HG1 and HG5 exhibited significantly greater tail beat frequencies during dusk ascents, whereas HG2 and HG3 exhibited higher tail beat frequencies during dawn descents. All individuals had significantly higher amplitude of acceleration (tail beat “effort”) and ODBA during dawn descent than dusk ascents (S2 Table).

Temporal, physiological, and environmental factors were examined for their influence on ODBA. The first principal component (PC1) comprised of depth, ambient water temperature, intramuscular temperature, and dissolved oxygen saturation explained 89.1% of variance (S3 Table) and was included with time of day and swimming phase in the GAMM that best fit the data (S4 Table and S5 Table). The model explained 20.7% of deviance in ODBA by sixgill sharks indicating there are additional factors beyond our model explaining the majority of variation in their activity (S4 Table). Residual analyses indicated a satisfactory fit for small to moderately large ODBA values; however, a departure from normality was observed at very large ODBA values (i.e., bursts) and therefore must be interpreted with caution (S2 Fig).

GAMMs revealed sixgill sharks exhibit higher rates of activity during descent compared to ascent and level swimming phases. GAMMs further identified distinct diel patterns in activity by sixgill sharks with an increase in activity from 02:00–10:00 (peak ~07:00) and 12:30–18:00 (peak ~15:00) and reduced activity levels from 18:00–01:30 (low ~21:30) (Fig 6). A peak in activity at PC1 scores ranging from -2.3 to -1.3 corresponded to the shallower nighttime distribution at a median depth of 294 m (IQR 249–277 m), median water temperature of 12.6°C (IQR 12.0–13.3°C), median intramuscular temperature of 11.8°C (IQR 10.7–12.3°C), and median dissolved oxygen saturation of 69% (IQR 67–72%) (Fig 6). An equivalent peak in activity at PC1 scores ranging from 1.8 to 2.8 corresponded to the deeper daytime distribution at a median depth of 604 m (IQR 579–622 m), median water temperature of 5.9°C (IQR 5.6–6.2°C), median intramuscular temperature of 6.8°C (IQR 6.4–7.3°C), and median dissolved oxygen saturation of 18% (IQR 17–20%). ODBA had the greatest decline at PC1 scores >2.8 corresponding to a median depth of 759 m (IQR 747–774 m); however, this result is biased by one individual (HG5) that exceeded a depth of 700 m.

For GAMMs that included swimming phase and time of day with each of the correlated predictor variables, the models including swimming phase, time of day, and depth or water temperature best fit the data (ΔAIC_c = 0.8) and explained 15.4% and 17.9% of deviance in ODBA, respectively (S4 Table and S5 Table). Swimming phase and time of day had similar effects on ODBA as the model including PC1 (see above). There was a peak in activity at depths between 550–700 m with a secondary peak between 275–375 m, and activity was lowest at depths shallower than 275 m and below 700 m (HG5 only) (S3 Fig). The highest peak in activity was observed between temperatures of 11–13°C with a secondary peak between 5.0–6.3°C, and a rapid decline at water temperatures exceeding 16°C (S4 Fig).

HMMs were constructed by means of two behavioral states: relatively low levels of ODBA (state 1) and overall higher activity levels including ‘bursts’ in ODBA (state 2) (S5 Fig). The HMM that included individual shark, swimming phase, time of day, dissolved oxygen saturation and the interaction term (time of day × dissolved oxygen saturation) best fit the data (S6 Table). Similar to GAMM residual patterns, HMM pseudo-residual analyses indicated a satisfactory fit for small to moderately large ODBA values and a departure from normality at very large ODBA values (bursts) (S6 Fig). Despite the described lack of fit, the HMM captures most of the relevant structure and these bursts will still be assigned to the high-activity state (state 2) allowing for biologically meaningful inference [48].

Overall sixgill sharks spent more time in the less active state (range 66–91%) across their vertical range except for HG2 (47%). Sixgill sharks also remained in the less active state for longer durations (means 8.8–26.3 min, range 1.2–199.3 min) compared to the more active state (means 2.5–9.8 min, range 0.3–180.0 min). Similar to the GAMMs, HMMs identified clear diel patterns in activity with the majority of time spent in the more active state (state 2) during shallow night phases (range 27–60%) and deep day phases (range 29–66%) (Fig 7); however, differences among individuals were observed. Moreover, sixgill sharks spent higher proportions of their deeper daytime distribution (range 36–73%) in the more active state compared to shallower nighttime distributions (range 11–37%) except for HG5 (deep: 6%, shallow: 11%), which spent the least amount of time overall (9%) in the more active state (Fig 8 and S7 Fig).

We estimated the stationary (equilibrium) state probabilities as a function of individual shark, swimming phase, time of day, and dissolved oxygen saturation and set each covariate to their overall median value for deep day (time of day, 12:44:30; dissolved oxygen saturation, 19%) and shallow night phases (time of day, 01:03:30; dissolved oxygen saturation, 71%). Individual shark and swimming phase were set to HG1 and level, respectively. Overall sixgill sharks had a lower probability of being in the high-activity state during ascents when compared to descent and level swimming phases, and during shallow night phases compared to deep day phases. Sixgill sharks showed the highest probabilities of being in the high-activity state between 14:00–18:00 (up to 60% probability) and at low dissolved oxygen saturations below 19% (up to 65% probability) during deep day phases (S8 Fig). The probability of being in the high-activity state is approximately 30% at higher dissolved oxygen saturations (>45%) experienced during shallow night phases and increases in the hours following dusk ascent before declining again to a probability of ~20% prior to dawn descent (S8 Fig).

Discussion

We combined a novel suite of biologging technologies to examine the fine-scale swimming behavior, activity, and thermal physiology of a relatively understudied deepwater species in relation to its environment. Sixgill sharks undertook pronounced diel vertical migrations and spent a considerable amount of time in cold, low oxygen conditions during their deeper daytime distribution. Analyses of intramuscular temperature revealed sixgill sharks maintain constant heat-transfer coefficients during diel vertical migrations, but their large body mass provides sufficient thermal inertia to slow heat loss during their descent into cold, deep waters. Sixgill sharks exhibited relatively high rates of activity during both shallow night and deep day phases and, contrary to our predictions, did not reduce activity levels during their deeper daytime distribution while they were experiencing low temperature and dissolved oxygen. Further, the highest overall rates of activity were observed during dawn descents and deep day phases.

Biologging technologies overcome the challenges associated with direct observation of behavior, especially for large deep-sea species that are not amenable to captivity [55] and/or suitable for laboratory studies. In particular, high-resolution accelerometers provide a useful tool for examining behavioral responses of free-ranging animals to large changes in environmental conditions. Previous studies examined changes in activity for sixgill sharks using either vertical or horizontal displacement rates via acoustic [56–58] or satellite telemetry [10,59,60]. However, direct comparisons of vertical displacement rates (derived from one-dimensional depth data) and activity (derived from multidimensional acceleration data) have revealed distinct differences in diel patterns [61–63]. For example, [10] demonstrated reduced vertical displacement rates during the deeper daytime distribution of sixgill sharks in Hawai‘i and suggested more active behavior in the shallow nighttime habitat. Here, we demonstrate that although sixgill sharks do exhibit peaks in activity and are more vertically dynamic during their shallower nighttime distribution, they are most active during their deeper daytime distribution. We suggest the decrease in vertical displacement rates during the day may be due to selection of a preferred isobath or to forage on a prey field that is less mobile and dispersed through the water column than prey hunted during the night.

ODBA and amplitude of acceleration were greater during the dawn descent than the dusk ascent segments of diel vertical migrations for all sharks. By contrast, tail beat frequency showed an opposite trend. Although two sharks (HG1 and HG5) displayed higher tail beat frequencies during dusk ascents, consistently higher amplitude of acceleration suggests sharks were beating their tail harder [64] during descents and the higher ODBA values indicate more energy was utilized during descent segments, despite lower tail beat frequencies for two individuals. Although we did not measure and record swimming speed, our findings are consistent with the interpretation that sixgill sharks exert more energy on descent than ascent segments of vertical migrations and periodically glided ‘uphill’ indicative of positive buoyancy [11] (S9 Fig). In addition, sixgill sharks exhibited shallow descent (-1.5° to -11.8°) and ascent angles (5.0° to 8.6°) during vertical migrations typical of movement designed to minimize the cost of horizontal transport [51,65,66]. Intramuscular temperature was close to equilibrium with ambient water temperature at the time of dusk ascent and this cold temperature may necessitate reduced energy expenditure during ascents aided by intermittent passive gliding.

GAMM and HMM analyses revealed diel changes in activity with sixgills most active during their shallow nighttime and deeper daytime distributions. Most notably, GAMMs revealed that sixgill sharks were least active in the warmest waters (>16°C) encountered at shallower depths or in deep waters (>700 m) exposed to low temperature (<5°C) and oxygen conditions (10% saturation). These results are biased by only two individuals (HG3 and HG4) encountering temperatures above 16°C and one individual (HG5) descending below 700 m; however, the lack of occurrence of other individuals at these environmental ranges suggests they may be outside the preferred vertical habitat range for sixgill sharks in Hawai‘i. Moreover, the diel behavior of sixgill sharks in Hawai‘i is fairly predictable and high-resolution (3–15 s) archival records obtained from recovered pop-up satellite tags revealed similar habitat preferences over longer time periods (53–200 days; [10], D. M. Coffey unpublished data). Therefore, the behaviors we observed over shorter time periods are likely representative of long-term behavior. Although sixgill sharks in Hawai‘i have been observed at water temperatures up to 20.6°C, [10] suggested a thermal ceiling at temperatures above 16°C. In other locations, the warmest temperature at which sixgill sharks were observed was 16°C [56] and ~17°C [67]; however, the length of exposure to these temperatures was not reported. Our empirical results support the hypothesis of a preferred thermal habitat at water temperatures below 16°C. This is supported by the fact that the largest decline in activity occurs at temperatures above 16°C and observing that exposures into warmer waters were too short to considerably warm intramuscular temperatures. However, longer duration deployments are needed.

Interestingly, the highest probability of being active was observed during the deeper daytime distribution (500–650 m) at cold temperatures and low oxygen conditions. Sixgill sharks spent the majority of their deeper daytime distribution with intramuscular temperatures warmer than ambient water temperatures providing them with a significant thermal advantage over non-vertically migrating and smaller-sized prey. Therefore, thermal inertia may allow sixgill sharks to maintain higher activity levels for long periods at cold temperatures, and thus exploit deep food resources more effectively compared to more resident species that are at thermal equilibrium. Moreover, oxygen conditions below the thermocline in Hawai‘i are not very variable [10,68] and reach a mean (± SD) saturation of 12 ± 2% (minimum ~9%) between 600–800 m (Station KAHE, Hawai‘i Ocean Time-series [69]) corresponding to the lowest saturation values recorded in this study and [12]. The increase in activity at these depths demonstrates that sixgill sharks are highly tolerant of low oxygen levels occurring within the local OMZ. Though, this activity may incur an oxygen deficit requiring a return to shallower waters with higher oxygen content as suggested for diel vertically-migrating swordfish (Xiphius gladius) [6].

Although sixgill sharks exhibited relatively high rates of activity under cold, hypoxic conditions, their large body size and comparatively slow tail beat frequency (mean 0.12–0.21 Hz) and swimming speeds (mean 0.3–0.4 m s^-1 [11]) suggest they have low metabolic rates relative to more active, shallower-dwelling species [70,71]. Low metabolic rates could allow sixgill sharks to maintain aerobic metabolism under low oxygen conditions [17,18]. However, regional comparisons have revealed that some deep-sea species have low metabolic rates in regions with high oxygen concentrations and exhibit a reduced tolerance for low oxygen conditions (i.e., high critical oxygen thresholds, P_crit) [17,18]. Therefore, low metabolic rates are not necessarily an explicit adaptation to low oxygen conditions and may better reflect ‘sluggish’ behaviors that result from relaxed selection pressure for rapid locomotory capacities for visual predator-prey interactions in environments with limited visibility [72–75].

Given the greater energy yield of aerobic metabolism, we expect anaerobic metabolism (e.g., anaerobic glycolysis) to merely support activity levels above routine metabolic rates for short periods of time, as opposed to completely supporting metabolic expenditure throughout the daytime distribution under hypoxic conditions [17]. However, the large body size of sixgill sharks and consequent lower mass-specific metabolic rate may allow for greater glycogen stores and a slower accumulation of deleterious anaerobic end-products (e.g., lactate and hydrogen ions [H⁺]) compared to smaller species [76–78]. In addition, sixgill sharks may possess mechanisms for enhanced oxygen extraction and transport, which could enable them to maintain aerobic metabolic rates at the observed low oxygen conditions (i.e., low critical oxygen threshold, P_crit) [17,18]. For example, some vertical migrators possess high-affinity respiratory proteins (e.g., hemoglobin) that are both temperature and pH sensitive (i.e., Bohr effect), which facilitate oxygen extraction in cold, low oxygenated waters and enhance oxygen release in warm, high oxygen conditions at shallower depths [17,18]. In addition, some residents of OMZs possess adaptations to increase gill diffusion capacities in order to facilitate oxygen uptake in low oxygen conditions including large gill surface areas and thin blood-to-water diffusion distances [17,18,72,79]. Identifying the oxygen partial pressures (pO₂) at which blood is 50% saturated (P₅₀) (e.g., [80]) and the pO₂ at which sixgill sharks can no longer maintain aerobic metabolism at a level independent of ambient pO₂ (P_crit) [18] will better define the oxygen tolerance of sixgill sharks, though the latter is technically challenging for large-bodied species, particularly in the deep-sea.

The enhanced thermal inertia of sixgill sharks may confer ecological and physiological advantages at cold temperatures as suggested for fishes with endothermic capacities. Endothermic fishes have evolved specialized anatomical and physiological mechanisms for retaining metabolic heat enabling higher levels of activity for longer durations when diving below the mixed layer to exploit deep forage resources more effectively than ectothermic species [81–83]. However, some large ectothermic sharks also undertake diel vertical migrations into cold temperatures below the thermocline to exploit deep forage resources (e.g., [4,31,84–86]). Though, ectothermic blue sharks (Prionace glauca) and ocean sunfish (Mola mola) make frequent ascents to surface waters well before intramuscular temperatures reach equilibrium with ambient water temperatures experienced at depth and employ variable heat-transfer coefficients where intramuscular temperatures warm more quickly than it cools [9,84], possibly due to changes in blood flow rates as demonstrated by diel vertically migrating swordfish [87] and bigeye tuna (Thunnus obesus) [7,32]. In contrast, sixgill sharks demonstrate prolonged exposures to cold temperatures and low dissolved oxygen concentrations at depth and maintain constant heat-transfer coefficients with comparable rates of heating and cooling, which may be attributed to low metabolic heat production as a result of relatively low activity rates [88,89]. We suggest the combination of relatively slow swimming speeds [11], low activity (as a proxy for low metabolic rates), and large body size with enhanced thermal inertia may enable sixgill sharks to actively search for sparse or patchily distributed forage resources across distinct habitats during diel vertical migrations.

ODBA had many extreme values (i.e., bursts) that needed to be accommodated in both modeling approaches. This issue could be resolved in HMM analyses by estimating the state-dependent densities nonparametrically, particularly for the high-activity state [48,50,90]. Furthermore, the HMMs provide a data-driven, unsupervised learning approach to analyzing ODBA; therefore, we cannot directly link the two activity states to specific behaviors. For example, changes in activity may be the result of foraging, predator avoidance, and/or social interactions, and these behaviors may also not be exclusive to a single state. Moreover, GAMM analyses revealed the best fit model including swimming phase, time of day, and PC1 (comprised of depth, ambient water temperature, intramuscular temperature, and dissolved oxygen saturation) only accounted for 20.7% of deviance in ODBA. Therefore, a significant portion of sixgill shark activity could be explained by an unobserved ecological or physiological covariate such as sporadic feeding events that may be independent of the environmental variables measured here. The relative high rates of activity during both shallow night and deep day diel phases suggests sixgill sharks do not conform to a ‘hunt warm, rest cold’ or ‘hunt cold, rest warm’ behavioral routine [23,24,28]. Disentangling the movement and behavior of sixgill sharks, and other marine species, from their prey remains a significant challenge, especially with limited local dietary information (e.g., see [91]). Recent advances in biologging technologies capable of measuring ingestion and digestion [92] and integrating animal-borne cameras [11,93] would provide meaningful insight into other potential drivers of activity (e.g., foraging) and the adaptive significance of vertical migrations, as well as significantly improve interpretation of these modeling approaches [24,50].

Conclusions

The novel combination of biologging technologies used in this study provides significant insight into the behavioral and physiological ecology of a deepwater species. The different modeling approaches provide a robust framework for empirically testing the influence of environmental drivers on the activity of sixgill sharks across diel vertical migrations and empirically demonstrated their tolerance of cold, low oxygen conditions in deepwater habitats. Obtaining similar data in regions with more variable and/or lower dissolved oxygen content may help define the oxygen tolerance for this species and resolve the relative impacts of cold temperatures and low oxygen on behavior [17,72]. Significant advances in biologging technologies enable innovative in situ deep-sea natural experiments, thus improving our understanding of the biological response of deepwater species to contemporary and future ocean conditions.

Supporting information

S1 Fig [pdf]
Diel patterns in habitat use.

S2 Fig [pdf]
Diagnostic plots for the best-fit generalized additive mixed model on overall dynamic body acceleration.

S3 Fig [pdf]
Response curves and partial coefficients from the generalized additive mixed model including swimming phase, time of day, and depth.

S4 Fig [pdf]
Response curves and partial coefficients from the generalized additive mixed model including swimming phase, time of day, and water temperature.

S5 Fig [pdf]
State-dependent densities of overall dynamic body acceleration (ODBA) for sixgill sharks.

S6 Fig [pdf]
Diagnostic plots for the two-state hidden Markov model with the best fit on overall dynamic body acceleration.

S7 Fig [pdf]
Activity of sixgill sharks along a vertical profile of depth and ambient water temperature.

S8 Fig [pdf]
Stationary state probabilities from the best-fit hidden Markov model.

S9 Fig [pdf]
Swimming performance during a vertical movement.

S1 Table [pdf]
Results of Kruskal-Wallis test examining the influence of diel depth phase on vertical displacement rate.

S2 Table [pdf]
Results of Kruskal-Wallis test examining the influence of diel depth phase on pitch and acceleration metrics.

S3 Table [pdf]
Principal component analysis (PCA) loadings for components (1‐4).

S4 Table [pdf]
Ranked generalized additive mixed models of effects on overall dynamic body acceleration.

S5 Table [pdf]
Results from the best generalized additive mixed models fitted to overall dynamic body acceleration (ODBA).

S6 Table [pdf]
Ranked hidden Markov models of effects on activity states based on overall dynamic body acceleration.