Environmental cycles and individual variation in the vertical movements of a benthic elasmobranch

Trends in depth and vertical activity reflect the behaviour, habitat use and habitat preferences of marine organisms. However, among elasmobranchs, research has focused heavily on pelagic sharks, while the vertical movements of benthic elasmobranchs, such as skate (Rajidae), remain understudied. In this study, the vertical movements of the Critically Endangered flapper skate (Dipturus intermedius) were investigated using archival depth data collected at 2 min intervals from 21 individuals off the west coast of Scotland (56.5°N, −5.5°W) in 2016–17. Depth records comprised nearly four million observations and included eight time series longer than 1 year, forming one of the most comprehensive datasets collected on the movement of any skate to date. Additive modelling and functional data analysis were used to investigate vertical movements in relation to environmental cycles and individual characteristics. Vertical movements were dominated by individual variation but included prolonged periods of limited activity and more extensive movements that were associated with tidal, diel, lunar and seasonal cycles. Diel patterns were strongest, with irregular but frequent movements into shallower water at night, especially in autumn and winter. This research strengthens the evidence for vertical movements in relation to environmental cycles in benthic species and demonstrates a widely applicable flexible regression framework for movement research that recognises the importance of both individual-specific and group-level variation.


Introduction
Trends in the vertical movements of mobile species are widespread in marine ecosystems. They extend across ecological timescales, from the cycles of the tides (Gibson 2003;Krumme 2009) to those of the seasons (Hyndes et al. 1999;Milligan et al. 2020) and those exhibited over ontogenetic development (Grubbs 2010;Frank et al. 2018). They are also prevalent across trophic levels, from the diel vertical migration (DVM) of plankton (Hays 2003;Bianchi and Mislan 2016) to the diving patterns of cetaceans (Keen et al. 2019;Caruso et al. 2020) and pelagic sharks (Papastamatiou and Lowe 2012;Schlaff et al. 2014). An important area of research in marine ecology is to understand the drivers of these trends.
Among elasmobranchs, research on vertical movement has focused on pelagic sharks and their diving patterns using animal-borne data loggers (Speed et al. 2010;Papastamatiou and Lowe 2012;Schlaff et al. 2014;Andrzejaczek et al. 2019). This work has demonstrated associations between Responsible Editor: J. Carlson.
Sophie Smout and James Thorburn should be considered joint senior author.
164 Page 2 of 18 vertical movement and periodic environmental cycles, individual characteristics and ecological interactions. Associations between vertical movement and environmental cycles are termed depth-specific periodic behaviours . Over short timescales, many elasmobranchs exhibit tidal-driven movement, moving into shallower waters with incoming tides to minimise energy expenditure or predation risk or to forage (Ackerman et al. 2000;Wetherbee et al. 2007;Starr 2009, 2010). Over diel scales, solar light levels are often linked to vertical movement via DVM, with crepuscular or nocturnal movement into shallower or deeper water (normal and reverse DVM, respectively), and changes in vertical activity (Carey et al. 1990;Andrews et al. 2009;Doherty et al. 2019;Arostegui et al. 2020). In many cases, these movements are associated with prey availability, but links to thermoregulation and bioenergetic efficiency are also documented Sims et al. 2006;Papastamatiou et al. 2015). Over the lunar cycle, movements can be associated with changes in tidal range or moonlight for similar reasons (West and Stevens 2001;Graham et al. 2006;Weng et al. 2007). Likewise, longer-term seasonal variation in solar light levels through increases or decreases in day length (photoperiod) can be important, ultimately through the regulation or entrainment of physiological processes such as reproductive cycles (Carey et al. 1990; Kneebone et al. 2012;Nosal et al. 2013).
Alongside environmental cycles, biological characteristics are frequently associated with vertical movement (Papastamatiou and Lowe 2012;Schlaff et al. 2014). There is increasing interest in the role of individual variation and behavioural plasticity, but this remains poorly studied (Jacoby et al. 2014;Hertel et al. 2020;Shaw 2020). Sexual segregation in elasmobranchs is also understudied but has the potential to give rise to differences in vertical movement, for instance in association with the movement of gravid females into shallow, sheltered, inshore areas for parturition (Wearmouth and Sims 2010). Likewise, ontogenetic differences in vertical movements have been documented, especially in association with ontogenetic dietary shifts (Grubbs 2010; Thorburn et al. 2019). However, in many cases, the links between environmental cycles, individual characteristics and ecological interactions are difficult to demonstrate (Papastamatiou and Lowe 2012;Schlaff et al. 2014).
The vertical movements of benthic elasmobranchs, such as skate (Rajidae), have received limited attention and generalising patterns from pelagic species is challenging (Humphries et al. 2017;Siskey et al. 2019). Available evidence indicates that in some instances skate undergo DVM in association with diel changes in solar light levels (Peklova et al. 2014;Humphries et al. 2017), perhaps in response to diel rhythms in prey such as Norway lobster (Nephrops norvegicus) and teleosts such as hake (Merluccius merluccius) in the North Atlantic (Aguzzi and Sardà 2008;De Pontual et al. 2012;Brown-Vuillemin et al. 2020). In turn, the apparent influence of solar light suggests that lunar phase and photoperiod may be important, as reported for other elasmobranchs (Schlaff et al. 2014) and in marine ecosystems more widely (Migaud et al. 2010;Kronfeld-Schor et al. 2013). However, in general, the prevalence of DVM in skate and the influence of other environmental conditions remain poorly understood.
The flapper skate (Dipturus intermedius) is a large, benthic rajid, recently recognised as a distinct species within the common skate (Dipturus batis) species complex (Iglésias et al. 2010). According to the International Union for Conservation of Nature Red List of Threatened Species, this species complex is Critically Endangered (Dulvy et al. 2006). Once widespread in the North East Atlantic, the flapper skate has been strongly affected by overfishing (Brander 1981;Dulvy and Reynolds 2002) and is now principally found in small areas within its former range, including off the west coast of Scotland where the Loch Sunart to the Sound of Jura Marine Protected Area (LStSJ MPA) has been established for its conservation . Juveniles hatch at about 20 cm in length (Dulvy et al. 2006) and males and females reach maturity at 185.5-197.5 cm in length (Iglésias et al. 2010) and 14-21 years in age (Régnier et al. 2021), respectively. Mating is thought to occur inshore in early spring, after which point females move into shallow (25-50 m) water to lay eggs and males may move offshore (Day 1884;Little 1995;NatureScot 2021). Individuals grow throughout life, with females reaching more than 250 cm in length. During this time, individuals are thought to follow a predominately benthic lifestyle, feeding on a variety of prey, including crustaceans, teleosts and elasmobranchs (Steven 1947;Wheeler 1969).
The movement ecology of flapper skate is poorly understood, but several studies have examined their vertical distribution and depth preferences Thorburn et al. 2021). Modelling suggests a preference for depths from 100-300 m in close proximity to the coast . Within the LStSJ MPA, Neat et al. (2015) reported a daily depth range of 50-180 m for three tagged individuals. More recently, Thorburn et al. (2021) aggregated data from multiple studies to examine long-term (seasonal and ontogenetic) variation in depth use in terms of how well the MPA covers skate depth preferences. They showed that skate used depths from 1-312 m but core depth ranges (20-225 m) varied seasonally and across size classes, with shallow-water (25-75 m) use increasing in winter, especially by large females. This pattern did not appear to be driven by temperature, but other drivers of depth use were not investigated. However, only two studies have used longitudinal vertical movement time series to examine individual-specific trends in vertical movement and their drivers (Wearmouth and Sims 2009;Pinto and Spezia 2016). These studies suggest that skate may exhibit depth-specific periodic trends in relation to tidal, diel and lunar cycles, but the small number of analysed individuals (N = 4-6) limits current understanding of individual-specific trends in vertical movement across ecological timescales.
The aim of this study was to characterise the vertical movements of flapper skate in relation to environmental cycles and individual characteristics. Following previous research (Wearmouth and Sims 2009;Pinto and Spezia 2016), it was hypothesised that skate vertical movements would be associated with (1) tidal, (2) diel, (3) lunar and (4) seasonal cycles. Archival tag records collected from the LStSJ MPA were used to examine evidence for these environmental cycles and to quantify the relationships between vertical movements and environmental cycles. A flexible regression method was used for this analysis that enabled both individual-specific responses to environmental cycles and similarities in the responses of different individuals with shared characteristics to be investigated.

Study site
The LStSJ MPA occupies 741 km 2 on the west coast of Scotland (Fig. 1). The coastline is punctuated by sea lochs, peninsulas and narrow channels. This surrounds a complex bathymetric environment that encompasses shallow areas (< 50 m deep), glacially over-deepened basins and channels up to 290 m deep (Howe et al. 2014). In surveyed regions, muddy sediments are widespread, while coarser, rockier sediments are found over smaller areas (Boswarva et al. 2018).
Hydrodynamic conditions across the MPA are resolved hourly by the West Scotland Coastal Ocean Modelling System (Aleynik et al. 2016). Tidal ranges vary between ± 0.5 m during spring tides to ± 5 m during neap tides. At any one time, spatial variation in tidal heights is typically less than 1 m but can reach up to 2.5 m in specific areas. Daylight hours are very similar across the area at any one time, extending from a minimum of eight hours in winter to nearly 20 h in summer. Peak variation in sea water temperature occurs between March (~ 7.5 °C) and August (~ 16.0 °C), with the upper 100 m experiencing + 1-2 °C of vertical stratification during the summer and autumn. Salinity in the upper layers (0-34.5 psu) is strongly associated with proximity to sources of freshwater input, which skate seem to avoid (Lavender et al. in review), but remains relatively stable elsewhere. The flow regime is complex, dominated by winds and tidal currents that interact with the bathymetry and incised coastline.

Tagging
Forty-five skate were captured, tagged and released in three locations within the MPA in 2016-17, following the methods of Neat et al. (2015) (Table S1; Fig. 1). Skate were captured using baited angling lines with barbless hooks. The total length (snout tip to tail tip) and disc width (wing tip to wing tip) were measured over the dorsal surface of each skate. Star Oddi milli-TD archival tags were programmed to record pressure (depth), to a resolution of 0.24 m and an accuracy of 4.77 m, and temperature (°C), to a resolution of 0.032 °C and an accuracy of 0.1 °C, every 2 minutes. Tags were attached externally, via two stainless steel pins through the leading edge of the left wing, and secured on a padded housing. Data from archival tags were recovered for recaptured skate whose tags were removed and returned. For these individuals, maturation status (immature, mature) was defined using a model for maturation with total length (Iglésias et al. 2010). In all cases, the predicted maturation status was unambiguous.

Data processing
Data processing, visualisation and modelling were implemented in R, version 4.0.2 (R Core Team 2020). For recovered time series, depth was calculated from pressure via pressure/1.01626. Vertical activity was calculated as the difference in depth between sequential observations. To focus on undisturbed vertical movements, data in a window around capture events that lasted until midnight on the day of capture were filtered from the time series (see Supplementary Information §1).

Environmental cycles
Tidal, diel, lunar and seasonal cycles were considered as putative periodic environmental drivers of vertical movement. Tidal heights (m) were extracted from predictions made every 15 min for Oban (56.41535°N, −5.47184°W) by POLTIPS-3 v. 3.6.0.0/13 and linearly interpolated to 2 min time steps. Diel cycles were parameterised with a metric of solar light levels or as diel period (day/night), depending on the model formulation (see below). For the former, sun angle (degrees above the horizon) was calculated for Oban at each time step using the suncalc package (Thieurmel and Elmarhraoui 2019). For the latter, 'day' was defined as the time between sunrise and sunset, and 'night' as the time between sunset and sunrise, evaluated at Oban, using suncalc. Thermal stratification was considered as an alternate predictor for diel cycles but not implemented because data exploration indicated diel cycles in vertical movement over winter (when thermal stratification was minimal). Lunar phase (rad) was calculated for each time step using the lunar package (Lazaridis 2014). Seasonal trends were parameterised in terms of photoperiod (hours between dawn and dusk on each day). Lengthening and shortening days with the same photoperiod were distinguished using a factor ('photoperiod direction') for the direction of change in photoperiod. While there are other correlated drivers of seasonal trends, such as temperature, photoperiod is the least variable across the study area and the most strongly associated with diel cycles through solar light levels. Furthermore, recent work suggests that depth use is not driven by temperature over seasonal timescales (Thorburn et al. 2021). The expression of these covariates at Oban assumes that skate remained around this area over their time at liberty and that the values of the covariates remained relatively stable across this area at any one time, which is supported by analyses of passive acoustic telemetry data (Lavender et al. in review) and environmental conditions (see 'Study site'). However, other meteorological and hydrodynamic variables that were considered as candidate drivers of vertical movement were not included in models given more substantial variability over space and the uncertain locations of individuals.

Exploratory data analysis
Depth time series for each individual were examined visually at multiple temporal scales and in relation to covariates to examine evidence for any relationships between vertical movement and environmental cycles. Functional principal components analysis (FPCA) was used to compare the structure and shape of observed time series (Ullah and Finch 2013;Wang et al. 2016), via the fda package (Ramsay et al. 2020). Under this approach, smooth functions are used to describe the shape of individual time series and FPCA is applied to these functions to examine their similarities and differences. This approach was applied to all individuals with sufficient data for each month from 17th March 2016 until 17th February 2017, during which time observations were always available for at least five individuals. For each 1-month period, individual time series were represented using 4th order B-splines with 1000 basis functions. FPCA was then applied to the smoothed time series and PC scores were visualised to examine how time series differed within and between life-history categories (immature/mature males/ females) (see Supplementary Information §2).

Depth models
Additive models were fitted to depth time series for each individual to investigate relationships between depth and environmental cycles. The female (1547) with the longest time series (772 days) was taken as an example individual for model development. A five-stage protocol was implemented: A. Depth time series were thinned to reduce serial autocorrelation.
B. An initial model of (thinned) depth observations in relation to environmental cycles was fitted. This model did not include an autocorrelation structure and is termed the 'AR0 DP model'. C. The autoregression order 1 (AR1) parameter was estimated from the autocorrelation function (ACF) of the model's residuals. D. The initial model was refitted using the estimated AR1 parameter to capture residual serial autocorrelation. This model is termed the 'AR1 DP model'. E. Using the AR1 DP model, model smooths, predictions and residual diagnostics were examined.
The first step was to thin the individual's depth time series to reduce serial autocorrelation (A). Thinning was implemented by selecting every 30 th observation (accounting for breaks in the time series due to the removal of data around recapture events), since this degree of thinning was sufficiently large to reduce autocorrelation while sufficiently small to not influence model smooths.
An initial Gaussian additive model (the 'AR0 DP model') was fitted to the thinned depth observations in relation to smooth functions of tidal (height), diel (sun angle), lunar (phase) and seasonal (photoperiod) cycles (B). Tidal height was included to account for temporal changes in water depth and to test the hypothesis that vertical movement differs either side of slack tide. A three-way interaction between sun angle, photoperiod and photoperiod direction was included to examine DVM in relation to solar light levels with changes in day length. Following previous results (Wearmouth and Sims 2009;Peklova et al. 2014;Pinto and Spezia 2016;Humphries et al. 2017), it was hypothesised that skate would occupy deeper depths during the day when solar light levels are higher. Weaker, more concentrated DVM was expected in summer when nights are shorter and brighter. The incorporation of photoperiod direction in this term allowed the interaction to differ depending on the direction of change in photoperiod. An interaction between sun angle and lunar phase was also included to examine lunar cycles and test the prediction of this hypothesis that DVM will weaken with increasing moonlight.
The model was implemented using the bam function in the mgcv package (Wood et al. 2015) (see Supplementary Information §3.1). Tidal height, sun angle and photoperiod were implemented with cubic regression splines. Lunar phase was implemented using a cyclic spline, with boundary knots at (0, 2π). Interactions were modelled using tensor smooths. The parameter was set to 1.4 to limit overfitting. A basis dimension of k = 10 was sufficient (according to the k-index diagnostic) for most model terms, while k = 5 × 5 was sufficient for interacting terms.
An AR1 correlation structure was used to capture residual serial autocorrelation (see Supplementary Information §3.1). To incorporate an AR1 structure, the AR1 parameter was estimated from the ACF of the residuals from the AR0 DP model (C). The model was refitted with the estimated AR1 parameter (becoming the 'AR1 DP model') (D). The final model of depth (DP) was: where is the model intercept, s is a smooth function, te are tensor product interactions, p is the AR1 parameter, e refers to residuals and t indexes observations.
Following model fitting, model smooths, predictions and diagnostics were analysed for the AR1 DP model using standard mgcv functions. Posterior simulation was used to quantify the magnitude of covariate effects ) (see Supplementary Information §3.2). Expected values of the response (depth) were simulated from the model at contrasting values for covariates of interest (such as low versus high values of tidal height), with other covariates fixed at appropriate values. Posterior distributions generated for contrasting values were compared, with the median difference in simulated depths between these values taken as a measure of effect size.
For all subsequent individuals with sufficient data (N = 17), the same protocol (A-E) was followed to generate a set of 'AR1 DP models' that were analysed in the same way. Models with an identical structure were fitted to enable comparisons among individuals, although the term for photoperiod direction was dropped if there were insufficient data. Using the set of models, ensemble-average effect photoperiod t , photoperiod direction t + te sun angle t , lunar phase t +̂e t−1 sizes (for example, the average expected change in depth between low and high values of tidal height across all individuals) were estimated for each covariate, by summarising pooled posterior simulations from all individuals, to support graphical interpretation of model predictions (see Supplementary Information §3.3). In all cases, ensemble averages overlapped with the posterior distributions from multiple models. FPCA was used to examine how estimated smooths varied in shape among individuals, by representing smooths as B-splines and then applying FPCA to those functions, as described for the observed time series (see Supplementary Information §3.4).

Vertical activity models
For modelling, the total, absolute vertical activity in day/ night, averaged over the duration (hours) of each period, was taken as a response variable (VA, mhr −1 ). This representation of vertical activity increased the signal to noise ratio in each time interval, improved the normality of the response and reduced data volume, while accounting for the effect of diel duration on total vertical activity. For each individual, VA was modelled in relation to diel period (day/ night), lunar phase, photoperiod and photoperiod direction, by fitting an 'AR0 VA model' to the VA time series, estimating the AR1 parameter and then fitting the 'AR1 VA model', defined below: where period=day is the model intercept; period=night is the difference in intercept at night; ∅ t is an indicator variable which takes a value of 1 when period t is night and 0 otherwise; and period × photoperiod direction is a factor that distinguishes both diel period and photoperiod direction. As for the depth time series, 'AR1 VA models' fitted to each individual were then analysed using standard mgcv functions, posterior simulation and FPCA (see Supplementary Information §4).

Archival time series
Archival time series, comprising 4,337,193 observations, were successfully recovered from 21 skate (six immature females, nine mature females and six mature males;

Exploratory data analysis: visual inspection of the time series
Individuals used a wide depth range within 0-311.82 (median = 115.47) m ( Fig. 2 and S1-3). Vertical activity was generally low (median = 0.004 ms −1 ) but ranged from 0.00-1.13 ms −1 . Over seasonal timescales, individuals exhibited marked variation in vertical movement and consistent patterns were difficult to identify ( Fig. 2 and S1-3). Almost all individuals ranged extensively in depth ( Fig. 2a) but most also used narrow depth ranges over prolonged periods lasting from weeks to months (Fig. 2b). There was clear evidence for repeated movements to/from a narrow range of depths. Repetition was most noticeable over daily-weekly timescales and lasted for prolonged periods in four immature females (1538, 1509, 1522 and 1558), two mature females (1512 and 1507) and three mature males (1511, 1520 and 1523) (Figs. 2c and S1-3). In general, vertical movement was more restricted in spring and summer, when individuals tended to use deeper water, compared to autumn and winter, when higher vertical movement was associated with increased use of shallower water (Fig. 2a-c). Over shorter timescales, individual time series were extremely variable, with periods of relative stasis, sometimes lasting more than one week, punctuated by high vertical activity (Fig. 2d). For most individuals, there were times when vertical movement coincided with diel cycles, but this pattern appeared to be irregular through time and varied within and among individuals (Fig. 2e).
Individual time series consistently emerged as distinct in monthly FPCAs (Fig. S4). In each case, the variation among time series was relatively difficult to characterise effectively with few dimensions, with first and second PCs explaining 37.29-56.74 and 13.76-34.42% of variation, respectively. There was no evidence for strong similarities within lifehistory categories at any time; however, there were relatively few individuals within each group.

Exploratory data analysis: correlations between depth and covariates
Exploratory analysis suggested relationships between depth, sun angle and photoperiod (Fig. S5) but not tidal height or lunar phase. Most individuals used depths from near the surface (< 25 m) to greater than 200 m in both day and night throughout the year. However, depth at night was shallower on average, especially in autumn and winter when the median difference in depth between day and night was 75 m. This seasonal difference appeared to be driven by fewer visits to the shallowest depths (< 25 m) during the middle of the day, alongside fewer visits to depths beyond 200 m, in autumn and winter.

Depth models
An additive model of depth in relation to tidal, diel, lunar and seasonal cycles for individual 1547 estimated statistically significant effects for all smooth terms (Table S2; Fig. 3 and S6). The small, positive effect of tidal height on depth was consistent with the rise and fall of the tides (Fig. 3a). Sun angle had a clear influence on both tensor smooths (Fig. 3b-d). The sun angle × lunar phase interaction was driven by a strong diel pattern throughout the lunar cycle (Fig. 3c). Even at night, lunar phase did not appear to have influenced depth. However, the estimated effect of sun angle changed over the course of the year (Fig. 3c-d).
In winter, around the winter solstice (photoperiod = 8.46 h), there was a clear diel pattern with movement from shallower water (~ 75 m) at night to deeper water (~ 120 m) during the day (Fig. 3c). The diel cycle weakened over spring and broke down almost entirely in summer. By the summer solstice (photoperiod = 19.89 h), the expected depth was deepest and remained deep throughout most of the day (Fig. 3c). After the summer solstice, the expected depth remained deep (~ 120 m) over both day and night (Fig. 3d). In autumn, the diel pattern resumed relatively rapidly, strengthening in winter.
The relationships between depth and environmental cycles for individual 1547 were borne out for other individuals, despite substantial intra-specific variation (Fig. 4). There was a consistently small, positive effect of tidal height. Across all posterior simulations, the median difference in depth between the lowest and highest tides was 5.05 m with a median absolute deviation (MAD) of 7.01 m; however, differences from −3.60-26.26 m were supported (Fig. 4a). In general, sun angle was positively associated with depth ( Fig. 4b-c). This effect was strongest in winter (low photoperiod) when the median difference between selected low and high sun angles (−40.42 and 26.93°, respectively) was 24.46 ± 26.99 m. However, there was considerable variation among individuals in the strength and shape of this pattern, with limited associations (< 10 m change) between sun angle and depth for four individuals (e.g., 1509 and 1520), moderate associations (10-50 m change) for 11 individuals (e.g., 1502, 1511) and strong associations (> 50 m change) for two mature females and one mature male (1507, 1509 and 1536) ( Fig. 4b and S7). The diel pattern was weaker in summer when the median difference in depth between low and high sun angles (−10.14 and 57.01°) was 9.37 ± 12.34 m and there were no strong associations between depth and sun angle ( Fig. 4c and S7). There was no evidence that lunar phase substantially influenced the depth of any individual over extended periods during the night or day (Fig. 4d-e). Similarly, depth was only weakly associated with photoperiod on average, with depths on the longest day 8.73 ± 20.79 m and 22.42 ± 36.04 m deeper than when the photoperiod was shorter (13.22 h) earlier or later in the year, respectively ( Fig. 4e-f). However, while the effect of lengthening photoperiod was small (< 10 m change) for ten individuals (e.g., 1502 and 1507), for five individuals (e.g., 1511 and 1522) the effect was moderate (10-50 m) and for one mature female and two males (1507, 1536 and 1518) the shift in the mean depth with photoperiod exceeded 50 m. Individual 1536 exhibited the strongest contrast in the expected depth between low and high photoperiods (Fig. S8). For 13/18 individuals with sufficient data, these effects mirrored the effect of declining photoperiod as the days shortened, although there were two exceptions (female 1519 and male 1536) (Fig. S8). In all cases, variation in the shape of these smooths appeared principally attributable to individual-level, rather than group-level, variation, with no well-defined clusters attributable to life-history category visually apparent from FPCA scores (Fig. S9). The deviance explained by individual models ranged from 3.74-56.45 (median = 29.15) %. There was a moderate correlation between deviance explained and the duration of individuals' time series (Spearman's rank correlation, rs = 0.49, N = 18, p = 0.04).

Vertical activity models
VA ranged from 2.05-177.08 (median = 74.79) mhr −1 . For individual 1547, the estimated mean VA during the day was 69.42 ± 1.27 mhr −1 (standard error) ( Table S3). The model of VA (Table S3, Fig. 5 and S10) for this individual suggested that VA was marginally higher at night (difference = 3.80 ± 0. 87 mhr −1 , t = 4.37; Fig. 5a; Table S3). During the day, VA was approximately 2-5 mhr −1 higher during the new and full moon, around spring and neap high tides, but there was no clear effect of lunar phase at night (Fig. 5b). In both day and night, posterior simulation showed that VA was 9.09-57.75 mhr −1 (5th-95th quantiles) lower at higher photoperiods ( Fig. 5c-d). For example, with increasing photoperiod from low values when the individual was tagged (8.46 h) to the summer solstice (19.89 h), the median decrease in VA during the day was −47.93 ± 6.66 mhr −1 MAD-a 69.05% decrease relative to the mean VA during the day across the time series. Moving towards the winter solstice, VA during the day partially increased, but less substantially (17.70 ± 6.72 mhr −1 ).
However, relationships between VA and environmental cycles varied among individuals (Fig. 6). The distribution of estimated differences between VA at night versus the day was right-skewed, ranging from −1.58-10.93 (median = 3.66) mhr −1 or −5.37-32.65% of daytime levels (Fig. 6a). Half of the individuals showed little (≤ 5 mhr −1 ) difference in VA between night and day, but higher (> 5 mhr −1 ) VA at night was clearly evident for others, especially females 1509 and 1538 for which VA at night was > 25% higher than daytime levels. Lunar phase was not strongly related to VA (Fig. 6b-c). Nine individuals showed weakly cyclical relationships between lunar phase and VA, for specific combinations of explanatory variables is shown, with other variables fixed where necessary at: tidal height = 2.50 m, lunar phase = 0π radians, photoperiod = 13.22 h and photoperiod direction = 'lengthening'. In b, lunar phase is shown in radians; 0π, π/2, π and 3π/2 refer to new moon, first quarter, full moon and last quarter (see inset). In c-d, expected depths are defined by lower colour bar and marked by 5 m contours. Diagonal stripes mark observed combinations of sun angle and photoperiod with VA approximately 2-5 mhr −1 deeper at new moon and/ or full moon, but these effects were small. The effect of photoperiod differed among individuals (Fig. 6d-g). For most (13/18) individuals (e.g. 1502 and 1507), daytime VA declined with photoperiod as the days lengthened; in nine cases (e.g., 1509 and 1519), only weakly (< 10 mhr −1 ); in eight cases (e.g., 1502 and 1507) moderately (10-50 mhr −1 ); and in one case (1536) substantially (> 50 mhr −1 ) ( Fig. 6d and S11). In contrast, two individuals (1523 and 1533) increased VA by more than 10 mhr −1 with photoperiod, while three remaining individuals (1520, 1539, 1538) did not show consistent patterns (Fig. S11). In general, these patterns mirrored the expected changes in depth after the summer solstice as photoperiod declined, although there were three exceptions (1502, 1536 and 1548) (Fig. S11). In all cases, trends in day and night were similar (Fig. 6e, g). There was no clear evidence for clusters in the change in VA with environmental variables by life-history category from visual inspection of FPCA scores (Fig. S12). Deviance explained ranged from 0.48-77.23 (median = 20.50) % and

Discussion
This study describes the vertical movements of a critically endangered benthic elasmobranch and demonstrates the influence of some environmental cycles alongside individual variation for understanding these movements over daily and seasonal timescales. Movement patterns span periods of low activity that are sometimes restricted around core depths for weeks or months as well as periods of higher activity, providing additional evidence that benthic species may switch between periods of low and high activity over extended timeframes (Kawabe et al. 2004;Wearmouth and Sims 2009;Humphries et al. 2017). The results further develop the evidence base for vertical movements in relation to environmental cycles in benthic species (Wearmouth and Sims 2009;Peklova et al. 2014;Humphries et al. 2017), with a significant portion of the variation in flapper skate movement structured over daily and seasonal scales. For flapper skate, this study provides the first evidence that interactions between environmental cycles affect vertical movement, extending previous research that has mainly focused on the independent effects of specific predictors (Wearmouth and Sims 2009;Pinto and Spezia, 2016). In line with recently documented seasonal trends in average depth use (Thorburn et al. 2021), there is a substantial seasonal shift in diel vertical movement patterns, with normal DVM and elevated nocturnal activity weakening in summer when skate tend to be less active and spend more time in deeper water. This points towards solar light as an important correlate of vertical movements and may have implications for their vulnerability as bycatch in other parts of their range (Bendall et al.

Fig. 5
Parametric terms and smooths from additive model of total vertical activity per hour (VA) in relation to environmental cycles for female 1547. a Diel phase, b lunar phase; c diel phase x photoperiod (lengthening); and d, diel phase x photoperiod (shortening). In a mean expected VA in day/night is shown, surrounded 95% confidence intervals. In b-d smooths, which depict expected change in VA with covariates, are shown for day (light grey) and night (dark grey) and surrounded by 95% pointwise confidence bands. (Smooths centred so that confidence bands solely reflect uncertainty in the smooths and not overall mean.) In b lunar phase is shown in radians: 0π, π/2, π and 3π/2 refer to new moon, first quarter, full moon and last quarter, respectively (see inset)

Fig. 6
Predicted changes in total vertical activity per hour (VA) across environmental variables for all individuals. Each plot shows expected VA ( VA ) of each individual given: a effect of diel phase, b-c effect of lunar phase during b day versus c night, d-e effect of photoperiod, when the day length is increasing, during d day versus e night, and f-g, effect of photoperiod when the day length is shortening, during f day versus g night. In a, connected points and error bars mark mean expected VA ± 95% confidence intervals for each individual. In b-g, lines represent expected VA for each individual (centred, to enable comparisons, by mean of the function, ̂V A , for each individual) in grey (immature females), black (mature females) or blue (mature males) across the range of values for that variable experienced by each individual. Unless otherwise stated, predictions are shown for diel phase = 'day', lunar phase = 0π radians; and photoperiod = 13.22 h 2017), as noted for Arctic skate (Amblyraja hyperborea) (Peklova et al. 2014) and other elasmobranchs (Gilman et al. 2019;Siskey et al. 2019;Arostegui et al. 2020). There is also evidence for a role for individual variation in movement, alongside ontogenetic segregation in depth use documented previously (Thorburn et al. 2021).
At short timescales, tidal-associated movement patterns are associated with vertical movement in many elasmobranchs (e.g., Ackerman et al. 2000;Wetherbee et al. 2007;Starr 2009, 2010). Yet proposed explanations for this behaviour typically hinge on the benefits of movement into shallow water with the incoming tide, either for foraging (Ackerman et al. 2000;Starr 2009, 2010) or predator avoidance (Wetherbee et al. 2007;Knip et al. 2011;Guttridge et al. 2012), which are unlikely to apply to flapper skate that tend to occupy deeper water (Wearmouth and Sims 2009;Neat et al. 2015;Thorburn et al. 2021). In agreement with this view, this study indicates that the vertical movements of tagged skate were not strongly influenced by tidal cycles at the resolution considered here. Nevertheless, tidal cycles may influence skate movement in other ways. For example, reports from recreational anglers of higher catches around slack tide suggest that tidal cycles may be associated with fine-scale changes in foraging behaviour. Alongside evidence for tidal influences on the movements of other elasmobranchs (Ackerman et al. 2000;Whitty et al. 2009;Campbell et al. 2012) and benthic flatfish (Metcalfe et al. 1990;Scott et al. 2016), this suggests that the influence of the tides would benefit from further investigation, particularly with accelerometry (Gleiss et al. 2013;Coffey et al. 2020).
DVM is widely documented among elasmobranchs, especially sharks (Carey et al. 1990;Andrews et al. 2009;Arostegui et al. 2020). While flapper skate occupied a range of depths over the diel cycle, this study adds to the evidence base for DVM in flapper skate (Wearmouth and Sims 2009;Pinto and Spezia 2016) and benthic predators more widely (Sims et al. 2006;Gorman et al. 2012;Harrison et al. 2013;Peklova et al. 2014;Cott et al. 2015;Humphries et al. 2017). Notwithstanding variation, tagged individuals tended to occupy shallower depths and exhibit slightly higher vertical activity levels at night, especially in the autumn and winter when depths were 75 m shallower on average at night. However, at short temporal scales the timing of DVM often appeared irregular, as noticed for other skate species (Peklova et al. 2014;Humphries et al. 2017).
In other elasmobranchs, including batoids, DVM has been attributed to thermoregulation and bioenergetic efficiency Sims et al. 2006;Papastamatiou et al. 2015) and foraging (Andrews et al. 2009;Arostegui et al. 2020). For flapper skate, given elevated DVM in winter when vertical temperature gradients are minimal, a role for foraging seems most likely. Studies of common skate suggest a broad, largely benthophagous diet comprising benthic invertebrates, cephalopods, teleosts and elasmobranchs (Steven 1947;Wheeler 1969). Switches between different foraging modes for these prey over a diel cycle are one potential driver of DVM. Low activity, opportunistic predation is common in benthic fish (Smale et al. 1995(Smale et al. , 2001Fouts and Nelson 1999). For flapper skate, depth and vertical activity time series suggest that this behaviour may predominate in deeper water in daytime, perhaps in a central place away from disturbances. However, the consumption of demersal and pelagic teleosts, especially by larger individuals, points towards the use of more active predation strategies too (Brown-Vuillemin et al. 2020). Coupled with lower activity in deeper water during the day, an increase in the prevalence of more active predatory behaviour at night could contribute towards DVM since nocturnal movement from deep habitats will naturally cause movement into shallower water. Cyclical trends in prey availability are likely to influence these movements, as noted for other elasmobranchs (Carey et al. 1990). For example, within the LStSJ MPA it is likely that skate respond to the diel rhythms of prey such as Nephrops (Aguzzi and Sardà 2008;Thorburn et al. 2021) and sandeels (Ammodytes spp.) (Winslade 1974). For example, in shallower (30-70 m) parts of the MPA, the increased availability of sandeels to flapper skate at night in winter, when they remain buried in sand, may contribute towards the emergence of DVM at this time. However, while these examples illustrate some of the potential ecological dynamics that may affect skate, further data on skate diets, prey availability and movement are required to examine fully the links between foraging and movement.
Alongside foraging, the avoidance of unfavourable nearsurface conditions that prevail in the day and/or summer may influence DVM. There is some support for the hypothesis that skate avoid higher light levels near the surface through seasonal changes in DVM. The absence of appreciable effects of moonlight appears to conflict with this hypothesis, but may be partly attributable to the impact of cloud cover. Other factors, such as seasonal sexual activity, have been postulated as explanations for seasonal trends in DVM (Carey et al. 1990), but their influence remains unknown for flapper skate.
In addition to seasonal patterns in DVM, there was evidence for repeated movements to and from a central depth over daily-weekly timescales for prolonged periods. These kinds of movements have been documented in several skate species (Humphries et al. 2017), but the longevity of the repeated movements illustrated here is exceptional. The pattern is reminiscent of central foraging or refuging behaviour (Humphries et al. 2017). In other taxa, this is typically associated with movement to and from resting, recovery or nesting sites (Boyd et al. 2014;Papastamatiou et al. 2018;Brost et al. 2020). In flapper skate, this behaviour is particularly obvious in immature females, which suggests that it is not associated with reproduction. More reasonable hypotheses include a spatiotemporal separation in the times or areas that are most suitable for resting versus foraging, which may result from the need to avoid intra-specific competition with large individuals (Humphries et al. 2016;Papastamatiou et al. 2018); territorial behaviour (Martin 2007); or personality traits, with 'shy' individuals remaining in familiar territory (Patrick and Weimerskirch 2014).
Individual characteristics, such as personality (Jacoby et al. 2014), age (Grubbs 2010) and sex (Wearmouth and Sims 2010), can certainly be important drivers of movement. The results presented here demonstrate a significant role for individual variation, with distinctive patterns such as the extensive vertical movements of female 1547 and the prolonged period of low activity of male 1548 not shared by other tagged individuals. One emerging explanation for this kind of variation is personality (Jacoby et al. 2014;Patrick and Weimerskirch 2014). Individual variation in small samples may also be attributable to characteristics that appear to be associated with specific individuals but are in fact shared more widely among individuals in a population. For example, in this study, there was no clear evidence for clustering in individuals' responses to environmental cycles by lifehistory category, which may partly reflect relatively small sample sizes for each life-history category, especially if differences only emerge at specific times (e.g., in association with reproduction). For this reason, future research would benefit from more data for all life-history categories, especially juveniles, for which data remain sparse.
Despite evidence for structure in the vertical movements of flapper skate, unsurprisingly a notable portion of variation remained unexplained by models. The absence of space from models is a significant limitation for benthic organisms: while models implicitly assume that individuals are 'free' to respond to changes in the value of covariates, benthic movements are restricted by bathymetry (Humphries et al. 2017). Building on previous research of flapper skate movement inferred from passive acoustic telemetry (Lavender et al. in review) and the analyses of vertical activity developed here, approaches that integrate these two sources of information to reconstruct possible movement pathways over the seabed would further advance our understanding of the ways in which bathymetry shapes vertical movement. The roles of oceanographic variables and anthropogenic factors, such as aquaculture farms, that are associated with the movements of other elasmobranchs (Dempster et al. 2005;Schlaff et al. 2014) could be also investigated with improved localisations of individuals, in tandem with new tag deployments and parameters derived from regional hydrodynamic models (Aleynik et al. 2016). However, biotic variables, such as prey movement and competition, are likely to be the primary drivers of movement (Papastamatiou and Lowe 2012;Schlaff et al. 2014). The lack of information on foraging in particular is a major ecological and conservation knowledge gap, not just for flapper skate but also more broadly (Papastamatiou and Lowe 2012). Coupled with improved estimates of individual space use, new data on these biotic processes should help to explain the complex movement patterns exhibited by flapper skate and the role of environmental conditions as indirect or direct drivers of these patterns.
The study of vertical movement has implications for skate management. For many skate species, bycatch is a major conservation concern (Dulvy et al. 2014;Oliver et al. 2015;Bendall et al. 2017). Most bycatch occurs in bottom trawls with a relatively low headline (Silva et al. 2012). The benthic habit of skate makes them particularly vulnerable to this gear, but vertical movements substantially affect their exposure (Bizzarro et al. 2014;Gilman et al. 2019) and catchability (Kynoch et al. 2015;Bayse et al. 2016). The addition of 'tickler' chains to the front of trawls to startle fish from the seabed into nets appears to increase skate bycatch in bottom-trawl fisheries (Kynoch et al. 2015). Consequently, in the LStSJ MPA, their use has been prohibited. The results of this study suggest that this measure is likely to have the biggest impact during the day and in summer, when vertical activity is lower. Later in the year, when skate are more active, tickler chain removal may have less of an impact. However, if higher activity is associated with movement off the seabed, the vulnerability of skate to trawling in autumn and winter may actually be lower, as noted for other species such as spurdog (Squalus acanthias) (Ellis et al. 2015). This highlights the need for high-resolution data on the vertical movements of skate in relation to the seabed, which could be provided by accelerometers and sonar tags (Gleiss et al. 2013;Lawson et al. 2015;Coffey et al. 2020). This research is essential for the effective implementation of fisheries management measures and spatial management approaches, such as marine protected areas, designed to achieve reductions in bycatch and improve skate conservation (Siskey et al. 2019).

Supplementary Information
The online version contains supplementary material available at https:// doi. org/ 10. 1007/ s00227-021-03973-1. skippering angling vessels for tagging work. Clinton Blight provided tidal predictions from POLTIPS-3. This is a single-user coastal tidal prediction software. The licence is sold by the National Oceanography Centre to end-users and universities. Samuel Iglésias kindly shared data and analyses of length-maturity relationships; David Miller provided statistical advice; and Francis Neat provided additional support and assistance with this research.

Availability of data and material
The data that support the findings of this study are available from Marine Scotland and NatureScot but restrictions apply to the availability of these data, which were used with permission for the current study, and so are not publicly available.
Code availability All data processing and analyses were implemented in R (R Core Team 2020). Code is provided in the Supplementary Information.

Conflict of interest
The authors declare that they have no conflict of interest.
Ethical approval The research was approved by the ethics committee of the University of St Andrews (number SEC21024). All regulated procedures involving animals were carried out in compliance with The Animals (Scientific Procedures) Act 1986 under the Home Office Project License number 60/4411 by competent Personal License holders.

Consent for publication All authors and tag data sources consent for publication.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.