Recovery, body mass and buoyancy: a detailed analysis of foraging dive cycles in the European shag

Please cite this article in press as: Carlsen, European shag, Animal Behaviour (2021), h Foraging dives in birds and mammals involve complex physiological and behavioural adaptations to cope with the breaks in normal respiration. Optimal dive strategies should maximize the proportion of time spent under water actively foraging versus the time spent on the surface. Oxygen loading and carbon dioxide dumping carried out on the surface could involve recovery from the consequences of the last dive and/or preparation in anticipation of the next dive depth and duration. However, few studies have properly explored the causal pattern of effects within such dive cycles, which is crucial prior to any assessment of optimal dive strategies. Using time depth recorders and global positioning system loggers, we recorded over 42 000 dives by 39 pairs of male and female European shags, Phalacrocorax aristotelis. Dives either involved a straight descent and ascent, presumably reflecting an unsuccessful search for prey, or a descent followed by horizontal movement followed by an ascent, presumably reflecting active hunting pursuit of pelagic prey. Males were larger than females, but we were unable to distinguish between sex effects and the nonlinear effects of body mass on dive behaviour. Path analysis showed that within-individual dive-to-dive variation in surface times can best be explained as recovery from the previous dive. As expected in a pelagic hunter with unpredictable dive durations, there was no evidence of anticipatory preparation of oxygen stores in predive surface durations. Among-individual variation in dives showed that body mass directly affected descent durations, but individual variation in all other dive and surface durations was driven by variation in descent duration, suggesting a critical role for dive depth in overcoming body mass-dependent effects of hydrodynamic/wave drag and buoyancy. Our analyses test for the first time certain critical assumptions for studies assessing optimal dive strategies in birds and mammals, thereby revealing new details and avenues for research concerning adaptive diving behaviour. © 2021 The Author(s). Published by Elsevier Ltd on behalf of The Association for the Study of Animal Behaviour. This is an open access article under the CC BY license (http://creativecommons.org/licenses/ by/4.0/).

Foraging in marine systems can be unpredictable for animals that obtain their food by diving, due to large variation in prey abundance and location (Schreiber & Schreiber, 1989). Many animals foraging in marine systems therefore have slow reproductive rate and late maturation, which are life history characters associated with stochastic and demanding environments (Ailsa, Bernie, & Richard, 2001;Huang, Chou, Shih, & Ni, 2011;Schreiber & Schreiber, 1989). We still lack detailed knowledge of the various behavioural and physiological adaptations involved in their foraging dive strategies (Green, Halsey, & Butler, 2005). For example, it is unclear when the respiratory surfaceetime costs of foraging dives are being paid for, through either preparation or recovery. Surface durations could be at least partly 'preparatory', as has been suggested following analyses of patterns of predive surface durations in shags and cormorants (Lea, Daley, Boddington, & Morison, 1996), penguins (Sato et al., 2002;Wilson, 2003), ducks (Butler & Woakes, 1979;Stephenson, Butler, & Wokes, 1986), guillemots (Elliott et al., 2008) and sea lions (McDonald & Ponganis, 2012), or entirely 'recovery' based as is common in diving mammals (Leeuw, 1996). A purely recovery-based diving strategy entails extra time spent on the surface after long dives, representing a 'lost opportunity cost' in terms of time that could have been spent hunting while prey was present (see Stephens, Krebs, Brown, Vincent, & Ydenberg, 2007). It has therefore been suggested that some diving animals instead use a strategy of extensive 'preparation' before dives (see species listed with references above, Butler & Woakes, 1979;Stephenson et al., 1986;Sato et al., 2002;Wilson, 2003;Wilson & Quintana, 2004;Ponganis, Meir, & Williams, 2010).
Determining whether dive durations are mainly driven by 'preparation' or 'recovery' is crucial for testing theoretical models concerning optimal foraging strategies, such as the marginal value theorem (MVT; Charnov, 1976) as one can estimate the cost of each dive (i.e. surface duration, plus travel durations to and from the foraging depth) versus the gain (i.e. dive duration or actual duration of active foraging). Predicting the optimal foraging dive versus surface durations (see Houston & Carbone, 1992;Kramer, 1988) has therefore been the topic of many studies (e.g. Foo et al., 2016;Stephens et al., 2007;Walton, Ruxton, & Monaghan, 1998). Beyond simple pulmonary exchange, these MVT predictions are crucially based on the expectation of an exponential increase in surface duration costs. This is because the rate of gas exchange that depends upon the partial pressure differences between the tissues and the atmosphere in blood as oxygen (O 2 ) loading (Carbone & Houston, 1996;Walton et al., 1998) and carbon dioxide (CO 2 ) expelling follows diminishing returns (Halsey, Reed, Woakes, & Butler, 2003). The time spent 'recovering' or 'preparing' for a dive should in this way be directly proportional to the preceding or subsequent dive duration (or at least its anticipated length), respectively (Kramer, 1988).
Strategies for adjusting surface durations to the corresponding dive duration may vary from species to species depending on respiratory and cardiovascular systems (e.g. avian versus mammal), diving physiology, the foraging strategy and the predictability of prey capture on a dive-to-dive basis. In dive preparation, more time and/or effort is invested in preparing for longer dives that are more likely to provide maximal rewards in terms of prey capture. Diving animals have the ability to prepare for longer dives by expelling more CO 2 and loading up on more O 2 (e.g. arterializing venous blood and increasing respiratory volume) beforehand, compared to resting levels (Ponganis et al., 2010;Sato et al., 2002), thus postponing the critical point of surfacing. Key differences in physiology between diving mammals and birds are important here (Butler & Jones, 1997). For example, unlike mammals, birds have a more efficient gas exchange using a cross-current system, allowing them to more efficiently load O 2 and reduce body CO 2 using less time (Scheid & Piiper, 1972), and this may affect predictions regarding strictly preparatory-versus recovery-based surface times. Through dive preparation, diving animals can potentially increase dive durations (Butler & Woakes, 1979) while staying well within their aerobic capacity for repeated dives within a bout (Wilson & Quintana, 2004). This strategy should be most efficient if the individual can accurately anticipate the necessary length of the next dive, because if not (i.e. the individual uses more O 2 than it has prepared for) then the cost of preparation will be added to the cost of recovery. The ability to prepare and store extra amounts of O 2 before a dive will always be limited based upon physiological traits such as body size, whereas recovery is more flexible but may more often involve greater accelerating energy and/or time costs (Scholander, 1940). However, preparation does have an additional cost due to increased buoyancy from increased O 2 stores during the following dive (Lovvorn & Jones, 1991;Watanuki et al., 2005). Preparation is not only about adaptively increasing dive duration, but also about reducing the cost of dives that are likely to be short, such as information-sampling dives used to locate prey without active foraging (Sato et al., 2002). To optimize levels of individual dive preparation, a diver must have good information about the depth, position and density of its prey, such as in the case of benthic-feeding ducks foraging repeatedly in the same location (Butler & Woakes, 1979;Stephenson et al., 1986). Thus, our hypothesis is that species should instead utilize more recovery-based dive strategies when hunting freely moving prey in more stochastic foraging environments with unpredictable prey depths and changing patch quality, locations and densities.
In recovery-based diving, each dive duration will have a corresponding postdive surface duration, which will be heavily dependent upon the duration and energetic demands of the previous dive (Carbone & Houston, 1996). During restitution periods, CO 2 needs to be exhaled to regain homeostasis with respect to blood pH and partial pressures of gases in blood and muscles (Krebs & Johnson, 1937). Diving animals have been observed utilizing a greater proportion of their stored O 2 than at resting state (Ponganis et al., 2010). As noted above, refilling O 2 stores in blood haemoglobin and myoglobin takes considerably longer with accelerating time costs, as compared with restoring lung and air sac deposits (Walton et al., 1998), and likewise getting rid of lactate is a far more timeconsuming and costly process than simply expelling CO 2 (Butler & Jones, 1997).
Many animals increase their dive durations by decreasing total metabolism and thus rate of O 2 usage, utilizing the 'dive response' (Irving, Solandt, & Solandt, 1935) leading to accumulation of blood lactate (Carbone & Houston, 1996;Scholander, 1940), which determines their aerobic dive limit (ADL, Kooyman, Wahrenbrock, Castellini, Davis, & Sinnett, 1980). Crossing the ADL increases possible dive durations, but with exponentially increasing costs in terms of recovery times this is largely not profitable in terms of foraging efficiency, especially for species dependent on a sequence of multiple back-to-back foraging dives (Wilson & Quintana, 2004). In lengthy dives, however, individuals may approach the limit of their ADL, which is largely determined by the size of body O 2 stores (Schreer & Kovacs, 1997), and therefore larger individuals are expected to dive for longer durations before reaching their ADL. Body mass also affects how animals experience certain physiological and environmental effects, such as rates of O 2 uptake, hydrodynamic properties (Liu, Kolomenskiy, Nakata, & Li, 2017;Lovvorn & Jones, 1991;Webb et al., 1998), and availability of fast mobile fish prey (Lovvorn, Liggins, Borstad, Calisal, & Mikkelsen, 2001). The differences in body mass and composition within and between species, and even between the sexes of sexually dimorphic animals (Cook, Lescro€ el, Cherel, Kato, & Bost, 2013), are expected to have an impact on individual physiological capabilities, including maximum O 2 storage capacity (Stephenson, Turner, & Butler, 1989), buoyancy (Lovvorn & Jones, 1991), and heat loss and hypothermia (Enstipp, Gr emillet, & Lorentsen, 2005). Buoyancy can affect optimal dive strategies with its varying costs depending upon dive depth and the ratios of body mass to pelage/plumage or pulmonary/air sac volume. These effects will interact with patterns of prey predictability and availability at different depths, and thus influence whether individuals need to utilize 'preparation' versus 'recovery' dive strategies.
Here we used the European shag, Phalacrocorax aristotelis (hereafter shag) as an example of a pelagic foraging diving animal (Cramp & Simmons, 1977) to examine the extent to which their foraging dives are based upon 'preparation' versus 'recovery', using detailed behavioural analyses of natural variation in the durations of different parts of the dive cycle (predive surface, descent, bottom, ascent and postdive surface). Shags are opportunistic foragers with substantial geographical and temporal variation in their main prey (Gr emillet, Argentin, Schulte, & Culik, 1998;Hillersøy & Lorentsen, 2012). However, they mostly hunt elusive pelagic fish prey, and so even if the depth of dives can be anticipated for some more benthic prey types in certain locations, the duration of most dives is expected to be very variable. In the shag, males are also 15e19% larger in body mass than females (Cramp & Simmons, 1977; this study), and dive up to 50% deeper (Christensen-Dalsgaard et al., 2017), allowing the sexes to potentially exploit different foraging patches, even at the same location. The buoyancy of diving shags decreases drastically with depth due to compression of respiratory air along with the compression of plumage air (Lovvorn & Jones, 1991). The plumage of shags and cormorants is semipermeable to water (Gr emillet, Chauvin, Wilson, Le Maho, & Wanless, 2005;Wilson, Hustler, Ryan, Burger, & Noldeke, 1992), which allows it to collapse during dives enabling individuals to become streamlined, efficient swimmers (Enstipp et al., 2005). Although a thin insulating layer of air is retained in the plumage during dives Lovvorn, Croll, & Liggins, 1999), an increase in diving depth drastically increases metabolic rates in shags (Enstipp et al., 2005).
Most previous studies have looked directly at the physiological changes before a dive, such as heartrate and respiratory rate without testing its effect on the following dive durations. We used path analysis (see Henshaw, Morrissey, & Jones, 2020) to specify the causal relationships between behavioural and physiological parameters within the dive cycle. Importantly, we applied path analysis to both within-and among-individual variation (see van De Pol & Wright, 2009) in a variety of dive cycle parameters (see Table 1). If dives are preparation based, we expected withinindividual variation in surface duration before a dive (predive duration) to explain variation in the length of all following dive parameters for the individual. If, however, the dives are recovery based, we expected the preceding within-individual dive durations to explain the following surface duration (postdive duration). Path analysis on the among-individual variation in dive parameters then allowed us to properly disentangle the effects of sex versus body mass differences on individual shag diving and foraging strategies.

Study Site
The Sklinna archipelago, situated about 20 km off the coast of Vikna in Trøndelag, Central Norway (65 12 0 N 10 59 0 E), holds one of the largest shag colonies in Norway with ca. 2000 breeding pairs in 2017.

Data Collection
The fieldwork was conducted during JuneeJuly 2013e2018, including 78 birds (39 pairs) over six different breeding seasons. Chick-rearing shags were chosen based on their nest accessibility and how 'protective' the pairs were, as those that aggressively stayed around the nest were easier to capture/recapture. Parental birds were fitted with loggers when nestlings were approximately 5e35 days old. Nestling age was determined using morphological criteria determined from control nests (from nesting areas in similar habitat within the Sklinna colony) checked every fifth day. The shags were captured and then recaptured at their nest by hand or using snares. Each individual was fitted with a GPS logger (i-gotU GT-120, Mobile Action Technology, New Taipei City, Taiwan; refitted in heat-shrink tubes) and time depth recorders (TDR, G5, CEFAS Technology Ltd, Lowestoft, U.K.). TDR loggers were attached to the GPS logger prior to instrumentation, and the loggers were attached to three to four middle tail feathers using TESA tape. The maximum logger deployment weight was 30.6 g, corresponding to 1.6% and 1.8% of mean body mass of males and females, respectively. The GPS loggers recorded location (±10 m) every 30 s, and the TDR recorded water depth below (±0.1 m) every 1 s. The loggers were removed during recapture after approximately 2e5 days. Deployment of loggers normally required less than 3 min of handling and retrieval less than 10 min, and no disturbance effects were noted in either adults or their chicks. In cases where there were signs of parental disturbance in the form of decreased nestling provisioning, then the second parent was not captured, and so these pairs were not included in the study.
The sex of adults was determined initially by body size features and ultimately via their vocalizations (Koffijberg and Van Eerden, 1995;Cramp & Simmons, 1977), because males and females made very distinct types of calls while defending the nest at our approach (Snow, 1960). At capture, body mass was obtained using a Pesola spring balance (accuracy ± 10 g). Both adults in the pair were fitted with recording instruments during the same breeding season, although not overlapping in time, usually within only a few days of each other. At recapture, biometric measures were obtained (wing length (ruler ± 1 mm), head and bill length (digital calliper ± 1 mm) and body mass (see above)). Adult female average mass was 1610 g (range 1370e1860 g), while average adult male mass was 1920 g (range 1660e2280 g). Growth data (i.e. captureerecapture difference in chick weight) were collected for all nests during the time of recording and these measurements were compared to the control area within the same colony (see above) containing 50 nests where adults were not fitted with loggers. There was rarely indication of parents reducing their provisioning rates or changing any patterns of nestling feeding while fitted with loggers. There were no obvious differences in the number of surviving chicks in experimental versus neighbouring control nests, aside from impeded survival due to gull predation.

Data Handling
Data handling and simulations were programmed in R 3.5.1 (R Core Team, 2018) and the TDR raw data were analysed with the package DiveMove (Luque, 2007). The total number of dives in this study was 46 103. The surface for dives was calibrated at ±1 m, so that no dive movement less than 1 m depth was counted as a real dive, which helped to remove possible nonforaging 'cleaning' dives (Christensen-Dalsgaard et al., 2017). The time submerged during foraging dives was divided into vertical descents and ascents involving <1 m horizontal movement versus >1 m horizontal movement. This horizontal movement was calibrated with the package DiveMove's zero-offset corrected (ZOC) method (Luque, 2007), smoothed using ±4 m depth filters, and registered as dive bottom duration. Dives were classified into two types according to the presence/absence of this horizontal dive bottom duration: U-shaped (with a horizontal dive bottom) versus V-shaped (with no horizontal dive bottom) dives (see Appendix 1, Fig. A1). Pre-and postdive durations at the surface longer than 360 s were used to separate dive bouts (i.e. distinct sequences of successive dives at one location) whenever surface durations were too long to be explained by simple replenishment of O 2 storage or momentary resting within a dive bout. GPS coordinates for each dive were assigned as the closest coordinates recorded within 30 s before and/or after the dive (i.e. GPS locations could not be recorded during dives). GPS data were processed using R library ggmap (Kahle & Wickham, 2013). Merging, combining and sorting of the data set were preformed using the package dplyr (Wickham, Romain François, Henry, & Müller, 2018), and plots were generated by ggplot2 (Wickham, 2016). A total 24 'locations' were identified as distinct places where most dives occurred (i.e. clusters of dives surrounded by areas with no dives), and distinguished as areas of uniform average depth and foraging conditions as determined from a topographical base map by Kystverket (https://kart. kystverket.no/). GPS coordinates were thus abbreviated to two decimal labels based on these dive locations, while the geographical size and number of observations varied between locations.

Statistical Analyses
Data analyses were performed mainly with mixed-effect models, conducted with R library lmer-function in lme4 (Bates, Maechler, Bolker, & Walker, 2015). Models were fitted to dive parameters: maximum depth, total dive duration, ascent duration, descent duration, bottom duration, predive duration, postdive duration and bout length (see Table 1 for a glossary). We fitted the fixed effects of sex, dive type (V-shaped or U-shaped) and body mass, while controlling for the random effects of individual ID, location, day and year. Body mass squared was included alongside body mass to detect any curvilinear effects but removed if nonsignificant. The random effect individual ID allowed us to control for all types of individual variation, including body mass when it was not a fixed effect. Similarly, any seasonal effects were also controlled for via the inclusion of the random effect of day. Estimated effect sizes are given as ±95% confidence intervals (CI) and random effects are presented as proportions of total variation explained. Akaike information criterion (AIC) values, P values and the principle of parsimony (i.e. the simplest model were chosen when two or more models were separated by an AIC score of <2) were used for model simplification (see Forstmeier & Schielzeth, 2011). Residual distributions for all models were checked and normality of data inspected with qqplot (Wickham et al., 2016).
Path analyses were based upon mixed-effect models with multiple fixed effects that thus controlled for adjacent path effects, similar to simpler partial correlations and regressions (Crespi & Bookstein, 1989). They allowed us to explicitly infer the direction of cause-and-effect relationships in extensive multivariate data sets. All variables in the path analyses were log-transformed to allow easier comparison of effect sizes within and between the different models. Log-log statistical comparisons also anticipated the types of exponential nonlinear effects predicted within dive cycles, by reducing everything to simple linear functions. The posttransformation linearity of the data was confirmed by assessing and rejecting nonlinear squared terms and inspected using ggplot2 (Wickham et al., 2016). The path analyses were based on a slightly reduced data set of 42 014 observations, where the first and last dives in each bout were removed to allow for only whole dive cycles involving both pre-and postdive durations. Among-and withinindividual effects were explored by using mean and meancentring dive parameter values, respectively, for each individual (see van De Pol & Wright, 2009). The effect sizes and CIs presented in the path diagrams are therefore summaries of the most important aspect of various mixed-effect models.

Effects of Sex and Dive Type
The effects of sex and dive type, as well as their possible interaction, were estimated in separate mixed-effect models for each of the main dive parameters separately (Table 2, Fig. 1; and see Appendix 2, Table A1). Dive depths ranged from 1 m to 63 m, dive durations 1 se154 s, descent durations 1 se89 s, ascent durations 1 se77 s, bottom durations 1 se106 s and surface durations 1 se365 s. There was a strong effect of dive type, with V-shaped dives being shorter than U-shaped dives, reflected in all the dive parameters, although interestingly not in surface durations. This validates our categorization of two types of dives equating to shorter prey-searching 'sampling' V-shaped dives where no prey was caught versus longer 'hunting' U-shaped dives where prey was actively pursued and probably caught. There was also a strong effect of sex on all dive parameters, including surface durations (Table 2, Fig. 1; and see Appendix 2, Table A1). This might be reflective of sex differences in body mass, although we could not statistically separate between these two effects (see Appendix 3, Table A2). There were also significant interactions between dive type and sex (Table 2), which complicates any interpretation of these different effects. For these reasons, all subsequent analyses were carried out separately for each sex and dive type.

Effect of Body Mass
Body mass had a significant effect on the depth and duration in both V-shaped and U-shaped dives in females, with effect sizes of over 4 m in mean depth and 5 s in mean durations for V-shaped dives between the lightest and heaviest females (Table 3; and see Appendix 2, Table A1). The effects in males were less clear, although male body mass had an effect on bottom duration, which then affected the total dive duration in U-shaped dives. Note, however, that any effect of body mass squared in females was negative (i.e. decelerating), meaning that overall the effect of body mass appeared to decrease at larger body mass values, perhaps   explaining the lack of body mass effects in the larger males. The heaviest female was estimated to have had dive durations equal to males of the same weight, emphasizing that the statistical effects of sex and body mass cannot be distinguished here (see Appendix 3, Table A2). The effects of body mass on descent versus ascent durations were similar in scale. Interestingly, there were again no direct effects of body mass on pre-or postdive surface durations in either sex (Tables 3, 4; and see Appendix 2, Table A1). However, there may have been indirect effects of surface durations (Table 2), emphasizing the need to untangle the cause-and-effect relationships between these surface durations and the other parameters with each dive cycle.

Path Analyses of Dive Cycles
Using path analyses, we investigated the effects of predive duration directly and indirectly on all other dive cycle parameters of descent duration, bottom duration and ascent duration, and thus all of their direct and indirect effects on postdive duration. The pathways based on the original log-transformed dive parameters (Fig. 2) therefore indicate the strengths of positive cause-and-effect relationships within the dive cycles separately for males and females, and for V-shaped dives versus U-shaped dives. For example, descent durations obviously explain much of the variation in ascent durations in both types of dives (adjusted r 2 values: V-shaped dives in females ¼ 0.53 and males ¼ 0.54; U-shaped dives in females ¼ 0.41 and males ¼ 0.39), due simply to both being linked through variation in dive depths. Predive durations also explain a lot of the variation in postdive durations in both types of dives (adjusted r 2 values: Vshaped dives in females ¼ 0.18 and males ¼ 0.35; U-shaped dives in females ¼ 0.36 and males ¼ 0.46), independent of variation in the durations of the various dive parameters (i.e. for comparison with univariate models showing separate path estimates not controlling for each other, see Appendix 4, Table A3, Fig. A2).
Positive predive duration 'preparation' effects on various V-shaped dive parameters were apparent for males, but not for females (Fig. 2a). Indeed, many of the effects throughout the V-shaped dive path analyses were lower for females than males, including some indications of 'recovery' in male postdive durations due to variation in descent durations. For U-shaped dives there were the same indications for both 'preparation' effects of predive durations and 'recovery' effects on postdive durations for males (and/or larger individuals), more so in females than in males (Fig. 2b). There was stronger evidence for 'recovery' in the effects of dive parameters on postdive durations in U-shaped dives compared with V-shaped dives for both sexes. Most significant effect sizes here were also biologically meaningful in scale when back transformed. For example, in female U-shaped dives the bottom duration effect on ascent duration of 0.02 corresponds to an increase of 1.02 s in ascent duration per s of bottom duration (Fig. 2b), and the largest effect of 0.99 for descent duration on postdive duration in male U-shaped dives corresponds to an increase of 2.69 s in postdive duration per s of descent duration. Note that in the univariate versions of these models (Appendix 4, Table A4) all paths had a strong effect on their own, but they were spread out more unevenly than in the full model path analysis controlling for all these different effects (Fig. 2).
One issue here is that the path analyses in Fig. 2 confound within-versus among-individual variation in dive behaviour, potentially obscuring some of the effects of interest here (see van De Pol & Wright, 2009). For example, the possibly sex-specific 'preparation' versus 'recovery' dive strategies could be explained by variation in behaviour across different individuals (e.g. due to The nonlinear term body mass squared was only included when significant. Fixed effects are given as effect sizes (for nonstandardized values) ± 95% confidence intervals in parentheses. Bold indicates significant values. Random effects of individual ID, location, year and day are given as proportions of total variation explained. See Appendix Table A1 for t and P values. body mass differences) or patterns of variation across dives within all individuals of the same sex. We therefore performed the same path analyses separately for mean-centered within-individual variation and for among-individual variation in the mean values.

Within-individual variation
The within-subject effects (Fig. 3) were notably different to the path analyses using the original dive parameters (Fig. 2). They were mostly significant and with similar effect sizes for both sexes, implying that any sex differences in the original analyses were driven by differences between individuals (Fig. 4) rather than sexspecific strategies within individuals. The effect of predive duration was much reduced in both V-shaped and U-shaped dives, suggesting that there was little 'preparation' going on in both sexes and dive types. In contrast, the within-individual effects of descent, bottom and ascent duration on postdive duration were much higher in V-shaped and U-shaped dives, suggesting 'recovery'-    Table 1 for definitions). Results are shown for (a) V-shaped no bottom duration dives and (b) U-shaped with bottom duration dives (see text for details). All parameters were log transformed. The strength of effect sizes (± 95% confidence intervals) are indicated as darkness of colours (lightest ¼ NS, light ¼ 1.00e-4e0.04, mid ¼ 0.05e0.25, mid-dark ¼ 0.25e0.55, dark ¼ 0.55e1.50), shown separately for females (red) and males (blue). See Appendix Table A3 for all of the statistical models involved and Table A6 for t and P values.
slight sex difference here, with larger 'recovery' effects in the different dive parameters on postdive duration in males, again because of their greater body sizes and deeper dives (see above).

Among-individual variation
There was some evidence for apparent 'preparation' in predive duration on descent in V-shaped dives in males but not females (Fig. 4a), and for both sexes in U-shaped dives (Fig. 4b). However, as we are considering among-individual variation here this cannot be 'preparation' as such and it is merely evidence that individuals with longer predive durations also had deeper longer dives (see Tables 3e6). This strongly suggests that any evidence in the original parameters path diagram (Fig. 2) for 'preparation' was due to the confound between within-and among-individual effects, because this effect is absent where we would expect to see evidence for it in the within-individual path analyses (Fig. 3). Interestingly, the among-individual effect of predive duration on postdive duration (Fig. 4) was greater than this effect for the original dive parameters (Fig. 2) and much greater than the effect for within-individual variation (Fig. 3).      Table A5). Results are shown for (a) V-shaped no bottom duration dives and (b) U-shaped with bottom duration dives (see text for details). All parameters were log transformed, except body mass. The strength of effect sizes (± 95% confidence intervals CIs) are indicated as darkness of colours (lightest ¼ NS, light ¼ 1.00e-4e0.04, mid ¼ 0.05e0.25, mid-dark ¼ 0.25e0.55, dark ¼ 0.55e1.50), shown separately for females (red) and males (blue). See Appendix Table A5 for the various statistical models and Table A6 for t and P values.
This suggests that any covariation here was driven by differences between individuals, with some individuals having longer surface durations than others in general. The same was also true for descent duration affecting ascent duration, which is much greater among (Fig. 4) than within individuals (Fig. 3), suggesting that the effect seen in the original data set (Fig. 2) was largely due to different individuals diving to consistently different depths.
Among-individual variation in dive behaviour was partly driven by differences in body mass within the sexes (see Appendix 4, Table A5). Fig. 4 therefore also includes all significant effect paths of body mass on each of these dive parameters. This neatly clarifies the apparently similar effects on descent and ascent durations in Tables 2 and 3, because in both V-shaped and U-shaped dives there is a direct effect of body mass on descent durations and any effect on ascent durations must have been indirect via the strong effect of descent durations on ascent durations. Male body mass had a negative effect on postdive duration, suggesting a mass-dependent advantage, such as increased buoyancy during ascent. If among-individual variation in body mass is driving the covariation of pre-and postdive durations, it is doing so mostly indirectly via total dive duration effects (Tables 3e6) on postdive duration. Notably, when selecting models for the amongindividual path analysis (Fig. 4), no interactions were significant (see Appendix 4, Table A5), meaning that body mass did not moderate any of the relationships between the other dive parameters.

Preparation-versus Recovery-Based Diving Strategies
We found no effect of within-individual dive-to-dive variation in predive surface durations in European shags, suggesting no evidence for a dive strategy based mainly on 'preparation', as has been found in other diving birds that seem to have anticipated dive durations (Butler & Woakes, 1979;Sato et al., 2002;Stephenson et al., 1986;Wilson, 2003;Elliott et al., 2008). Instead, our findings suggest that for shags any within-individual dive-to-dive variation in surface durations is likely to be due to 'recovery' from less predictable variation in dive durations, probably as a result of variation in local prey detectability, availability and capture rates. Although descent durations had the strongest effects, there were significant increases in within-individual postdive duration due to increases in the durations of all the different dive parameters, as expected if surface durations were responses to the need to recover the total cost of all aspects of the dive (Walton et al., 1998). Indeed, many of the same physiological abilities in respiration (see Frappell, Hinds, & Boggs, 2001;Lasiewski & Calder, 1971) that could be used to prepare avian divers for a lengthy dive would also allow them to be especially well equipped for efficient recovery after a dive (see above).
Given that shags search for and hunt freely moving fish prey in a highly variable foraging environment (Hillersøy & Lorentsen et al., 2012), anticipating dive durations may simply not be possible. The within-individual recovery effects of variation in descent, bottom and ascent durations on postdive surface durations were actually much stronger for U-shaped with bottom duration dives than for Vshaped no bottom duration dives (Fig. 3). This could be explained by greater dive-to-dive variation in U-shaped dive durations, as each of these dives will have involved active and unpredictable hunting pursuits (as seen in Weddell seals, Leptonychotes weddellii, see Kooyman et al., 1980), as compared to more standardized searching and sampling for possible prey during V-shaped dives. However, V-shaped dives were also variable in length within individuals, suggesting considerable dive-to-dive variation in the The nonlinear term body mass squared was only included in the models when significant. Fixed effects are given as effect sizes (for nonstandardized values) ± 95% confidence intervals in parentheses. Bold indicates significant values. Random effects of individual ID, location, year and day are given as proportions of total variation explained. See Appendix Table A1 for t and P values. depth of active searching and assessment of possible prey before surfacing without an active hunt. Therefore, V-and U-shaped dives might have shared similar characteristics earlier on in each dive prior to an actual pursuit taking U-shaped dives deeper and with a horizontal component, indicating a less clear initial biological division than our statistical categorization suggests (see also Cook et al., 2012). V-shaped dives may indeed just be aborted U-shaped dives when no prey was detected. Importantly, if the dive type is not decided upon before the dive is commenced, then the dive duration cannot be anticipated, and thus careful 'preparation' is unlikely to increase efficiency of foraging dives. For diving animals feeding on sedentary benthic species (e.g. diving ducks; Butler & Woakes, 1979;Stephenson et al., 1986), preparation may be more efficient as the individual can anticipate to some extent how long a dive should last, as information about the depth and the local prey availability would be easily accessible (Stephens et al., 2007). However, it would be far too much of a challenge for diving animals like shags foraging on solitary or shoaling pelagic prey to adjust predive durations to match dive-todive variation in time and energy costs, especially relative to the marginal benefits to be gained from such foraging (Green et al., 2005). Even if diving animals do not show dive-to-dive variation in the amount of preparation, it should still be expected that recovery results in individuals using increased respiration following a dive to return to their maximum sustainable O 2 storage more quickly (Kooyman, Kerem, Campbell, & Wright, 1973;Lasiewski & Calder, 1971) prior to the start of the next dive. Thus, any 'recovery' measured as the surface duration from the previous dive also represents the preparation needed to start any new dive at baseline levels of respiratory condition (Wilson & Quintana, 2004).
A nonmutually exclusive recovery effect could also have involved the need to restore body temperatures in between each dive by increasing metabolic rates for longer periods on the surface after increasing amounts of time submerged in cold waters (Enstipp et al., 2005;Gr emillet et al., 2001;Scholander, Hock, Walters, Johnson, & Irving, 1950). Both respiratory (O 2 storage and CO 2 release) and body temperature recovery provide broadly similar predictions, with escalating and possibly exponential curvilinear surface duration costs with increasing dive durations. Therefore, it is challenging to separate these two effects using the current data set (i.e. without accompanying body temperature data) as far as the within-individual dive-to-dive variation in dive cycle behaviour is concerned.
Whatever the physiological mechanisms involved were, our analyses emphasize the importance of distinguishing withinversus among-individual effects (see van De Pol & Wright, 2009). Indeed, many previous studies of this type in the literature carried out analyses similar to the original path analyses (Fig. 2), which gave the impression that there was preparation in predive surface durations determining dive durations (e.g. Lea et al., 1996). Once the analyses isolated the within-individual effects of pre-versus postdive surface duration (Fig. 3), in line with actual theoretical predictions, the effects of 'preparation' became difficult to discern. Likewise, when the different dive parameters were placed within a path analyses, it was possible to appropriately disentangle causeand-effect within each dive cycle. For example, what looked like similar effects of body mass on descent and ascent durations (Tables 3e6) turned out in the among-individual path analyses (Fig. 4) to be a direct effect of individual variation in body mass on descent duration that then itself influenced ascent duration independently of body mass (see below). Hence, any exploration of questions such as surface duration representing 'preparation' versus 'recovery' for dives, or of MVT optimal dive behaviour predictions, requires this type of statistical path analysis approach due to the covarying nature of dive parameters in a sequence of successive dives within and among individuals.

Among-Individual Effects and Variation in Body Mass
What appeared to be 'preparation' in predive surface duration in anticipation of total durations of dives in the original path analyses (Fig. 2) was completely due to among-individual variation (Fig. 4). Indeed, the random effect of individual ID in these mixed-effect models explained between 16 and 30% of the total variation in most dive parameters (see Table 2), which includes any variation in individual dive behavior due to body mass. We have shown that body mass had a strong but decelerating effect on dive parameters where almost all the effects of sex differences in our results overlapped entirely with effects of body mass, to the extent that the two factors were statistically inseparable for this data set (see Appendix 3, Table A2). Thus, as far as variation in foraging dive parameters are concerned, males simply represented a larger class of individuals than females.
The largest energy and thus O 2 drain for dives in shallow water is the sheer resistance in turbulent near-surface water due to movement (i.e. wave-induced hydrodynamic drag) and buoyancy (see Liu et al., 2017), which both decrease with depth and interact with the volume of the object down to a point of depth where they cease to have an effect (see Lovvorn & Jones, 1991). Increased buoyancy may thus explain the larger effect of descent on postdive duration in larger individuals, as they would have to work harder, possibly for a longer time, to descend. Indeed, it was the depth of a dive that seemed to be particularly affected by body size, and any subsequent dive parameters were only affected indirectly by mass (Fig. 4). Longer descent durations consequently led to the longer bottom, longer ascent and, together, to longer postdive surface durations for heavier individuals. This means that if among-individual variation in body mass was driving the covariation between pre-and postdive durations, it was doing so both directly and indirectly via descent duration effects on postdive duration.
The increased buoyancy of larger individuals may also have resulted in less energetically costly ascents based upon passive gliding (Lovvorn et al., 1999;Lovvorn & Jones, 1991), leading to relatively less energy and O 2 usage per s and thus quicker recovery once on the surface. This would explain why body mass had a negative effect on postdive surface duration after U-shaped dives. These results also argue against body temperature recovery as a major cost to diving for shags in these cold waters (see above). This is because heat loss and thermoregulation-related metabolic costs are body mass dependent in terms of surface-area-to-volume ratio, and so we should have seen similar effects of body mass on descent, bottom and ascent durations, and similar effects of each of these on among-individual variation in postsurface durations. The net result here in terms of travel (descent plus ascent) costs should therefore have been the same in terms of metabolic energy and O 2 use irrespective of body size, so long as individuals could choose their optimum dive depth for their body size. The deeper depths chosen by larger individuals thus imply that there was a reduced benefit and maybe even an increased energy or time cost for larger individuals when actively foraging (i.e. during bottom duration) at shallow depths. Such a cost might well have been due to increased hydrodynamic drag and buoyancy costs for larger bodies in water, due to larger air volumes in the respiratory air sacs and remaining plumage air layer (Lovvorn & Jones, 1991). This could then be avoided by diving deeper for larger (male) individuals.
Overall, these patterns of body mass-dependent dive depth variation among individuals, alongside the key role of descent duration on subsequent surface durations, support the notion that shags operate primarily under respiratory recovery from the previous dive, expelling CO 2 and loading O 2 ready for the next dive (Wilson & Quintana, 2004). However, the weak within-individual (Fig. 3) and strong among-individual (Fig. 4) link between preand postdive surface durations, independent of dive durations, would seem to indicate some form of temporal autocorrelations across dive cycles and surface durations. One cause of this might be accumulating effects, such as fatigue, within bouts that run across dive cycles depending upon body mass, challenging the assumption in most of the optimum diving strategy literature that each dive cycle is independent of the ones before and after it within a bout (see Carbone & Houston, 1996;Walton et al., 1998). Hence, surface 'recovery' in systems like diving shags is supposed to always be of sufficient duration to complete all CO 2 expelling and O 2 loading, and removal of any lactate build-up, before the start of the next dive (Lasiewski & Calder, 1971;Kooyman et al., 1973;Wilson & Quintana, 2004). Likewise, to avoid longer-term effects of repeated dives in cold water within bouts (Gr emillet et al., 2001;Scholander et al., 1950), body temperature should always be regained between dives (see above). To test this key assumption in studies of adaptive diving behavior and physiology, we are currently exploring the sources of any temporal autocorrelation to investigate any cumulative effects of fatigue or temperature on dive parameters across dive cycles within bouts in these data on European shags.

Conclusions
There was no evidence for 'preparatory' variation in predive surface duration that anticipates dive durations within individuals on a dive-to-dive basis. Instead, there seemed to be clear evidence of 'recovery' in almost all of the variation in postdive surface durations following considerable unexplained variation in the duration of both prey searching V-shaped dives and especially hunting U-shaped dives. This matches the stochastic foraging expected in such a pelagic hunter with unpredictable prey availabilities and dive durations. Our study therefore suggests that when applying MVT optimal foraging theories to within-individual variation in dive versus surface durations (see Carbone & Houston, 1996;Walton et al., 1998) it is critical to know the appropriate surface duration to compare with each dive.
There were no obvious sex differences in dive behaviour beyond the differences in body mass, which had its largest effect on descent duration, with all other dive parameters then being driven by descent duration, rather than directly by body mass. All of which suggests a critical role for wave drag and/or buoyancy costs in dive behavioural decision making in shags. It remains to be seen how much of the variation in surface times is explained by recovery of body temperature compared to O 2 loading and CO 2 expelling, and whether either of these is responsible for any cumulative effects that might exist across successive dives within a bout. Therefore, our analyses reveal important new details concerning adaptive diving behaviour in European shags, especially regarding tests by optimal MVT models of dive versus surface durations, but it also highlights important areas for future work in this area.

Appendix 2: t and P values for Tables 2e6 Appendix 3: Sex and Body Mass Effect Sizes
The effects of sex versus body mass (Table A2) on each dive parameter were determined by running linear mixed models separately for each factor. We then ran mixed models containing both sex and body mass, showing that all the effects of sex became nonsignificant when controlling for body mass, while body mass (being the more fine-grained covariate measure) explained the differences between individuals in all parameters. Body mass was able to explain the same variation as sex differences in all dive parameters. The effect of body mass was calculated from models for females only and multiplied by the mean weight difference between females and males of 300 g. This allowed us to compare the effect sizes for body mass and sex, which overlapped in size more or less completely, again suggesting that they represent the same biological effect and cannot be distinguished statistically in this data set. Effect sizes for models estimating: sex on its own; body mass on its own; both sex and body mass together; and quantitative comparisons of biological effect sizes assuming the usual 300 g difference between the sexes. Fixed effects are given as effect sizes (for nonstandardized values) ± 95% confidence intervals in parentheses. Bold indicates significant values. Random effects of individual ID, location, year and day are given as proportions of total variation explained.  Separated by sex for no bottom duration (V-shaped) and with bottom duration (U-shaped) dive types (see main text for further explanation). Mixed-effect models included individual ID, bout ID, location, year and dive date as random effects, which are given as proportions of total variation explained. Fixed effects are mean centered, given as log effect sizes (for nonstandardized values) ± 95% confidence intervals in parentheses. Bold indicates significant values except in intercepts. See Table A7 for t and P values.

Table A5
The among-individual effect of log dive parameters on log postdive duration, log ascent duration, log bottom duration and log descent duration. Separated by sex for no bottom duration (V-shaped) and with bottom duration (U-shaped) dive types (see main text for further explanation). Mixed effect models included individual ID, location, year and dive date as random effects, which are given as proportions of total variation explained. Fixed effects are mean values per individual, given as log effect sizes (for nonstandardized values) ± 95% confidence intervals in parentheses. Bold indicates significant values except in intercepts.     Figure A2. Path diagram summarizing the strength and direction of univariate parameters in dive cycles using log-transformed dive parameters. Results from univariate mixedeffect models are shown for (a) V-shaped no bottom duration dives and (b) U-shaped with bottom duration dives (see text for details). The strength of effect sizes (± 95% confidence intervals) are indicated as darkness of colours (lightest ¼ NS, light ¼ 1.00e-4e0.04, mid ¼ 0.05e0.25, mid-dark ¼ 0.25e0.55, dark ¼ 0.55e1.50), shown separately for females (red) and males (blue). Based on the full data set.