Diving into the vertical dimension of elasmobranch movement ecology.

Knowledge of the three-dimensional movement patterns of elasmobranchs is vital to understand their ecological roles and exposure to anthropogenic pressures. To date, comparative studies among species at global scales have mostly focused on horizontal movements. Our study addresses the knowledge gap of vertical movements by compiling the first global synthesis of vertical habitat use by elasmobranchs from data obtained by deployment of 989 biotelemetry tags on 38 elasmobranch species. Elasmobranchs displayed high intra- and interspecific variability in vertical movement patterns. Substantial vertical overlap was observed for many epipelagic elasmobranchs, indicating an increased likelihood to display spatial overlap, biologically interact, and share similar risk to anthropogenic threats that vary on a vertical gradient. We highlight the critical next steps toward incorporating vertical movement into global management and monitoring strategies for elasmobranchs, emphasizing the need to address geographic and taxonomic biases in deployments and to concurrently consider both horizontal and vertical movements.

INTRODUCTION

The aquatic realm is defined by its three-dimensional volumetric scale, facilitating vertical and horizontal movements of species on a daily basis compared to the more constrained two-dimensional space inhabited by terrestrial species [with the exception of flying and select climbing and burrowing fauna; moving with three versus six degrees of freedom; ()]. Selection of optimal thermal habitat, dynamic movements between temporally productive foraging patches, and predatory and competitive interactions occur across both vertical and horizontal planes. Defining vertical mobility in tangent with horizontal displacements of species is therefore a key factor to consider in understanding ecological-anthropogenic interactions in support of conservation and management of aquatic species. Ultimately, the vertical dimension of movement is as, if not more, important as horizontal movement for most aquatic species, with depth use influencing survival, fitness, and resilience to ongoing climatic impacts ().

While the horizontal distributions and movement patterns of marine megafauna (large fishes, cetaceans, pinnipeds, etc.) have been examined at a diverse range of temporal and spatial scales, facilitating global-scale comparative studies (, ), these studies in the vertical plane are limited. This is understandable given the very large size of datasets that would be needed for analyzing both horizontal and vertical space use at the global scale within a single study. These comparative studies, however, are particularly relevant to establish baselines of habitat use and generate hypotheses for threatened megafauna and taxa experiencing global population declines, such as sharks, rays, and skates (hereafter “elasmobranchs”) (). From an ecological perspective, vertical movement data can further our knowledge of the relative use of vertical bathomes and spatial and temporal overlap in habitat use within and among species, revealing the likelihood of intra- and interspecific interactions (). Given their mobility, elasmobranchs can connect disparate deep-sea and shallow water ecosystems in tangent with high- to low-latitude movements, regulating trophic interactions and nutrient cycling on localized to ocean basin scales (). When considering management and conservation, data on vertical distributions can inform us of a species’ vulnerability to anthropogenic threats, which is particularly relevant for elasmobranchs given that they can be subject to both targeted fisheries and bycatch worldwide (, ). For example, capture probabilities and levels of mortality are influenced by depth overlap between elasmobranchs and fishing gears (), as well as gear selectivity () and species physiology (). Data on vertical distributions can consequently inform encounter rates and catchability with various fishing gears and inform potential mitigation measures to reduce bycatch of vulnerable species while maintaining capture of fishery targets (). Similarly, data on vertical movements can facilitate and improve our understanding of the biases of both fisheries-dependent and fisheries-independent techniques for estimating species’ abundance, including fishery observer records, and aerial, scuba, camera [e.g., baited remote underwater video (BRUV)], and environmental DNA surveys. Knowledge of the vertical range of aquatic species is critical for informing dynamic spatial management planning efforts, such as the designation of protected areas where multiuse activities may occur, by ensuring that representative vertical, as well as horizontal, habitats are protected (, ).

Elasmobranchs are often tracked through a combination of biotelemetry and biologging devices. These electronic tags can record the depths encountered by tagged individuals, and, depending on the tag type, data are either directly downloaded once a tag is physically recovered (archival tags or recovered pop-up satellite archival tags; hereafter “archival tags”) or remotely transmitted via satellite when an individual surfaces or the tag detaches and floats to the surface at a preprogrammed date (“satellite tags”) (). For satellite tags, depth data can be transmitted as coarse time series and/or in a summarized binned histogram, where the sampling rate and histogram bin sizes are selected by the user a priori, whereas archival tags typically record high-resolution, continuous time-series data. To date, studies on elasmobranch depth distributions have generally been conducted on a species-specific basis (, , ). Yet, an early study () and recent reviews have identified consistent vertical patterns among certain taxa, including diel vertical movement (DVM) and deep-diving behaviors (, ). Furthermore, a recent synthesis of a regional predator assemblage in the South Atlantic revealed that dynamic species-specific space use was tied to water column structure (). Direct comparisons of vertical space use among elasmobranchs, however, are often hampered by inconsistent reporting of depth data. For example, the temporal period over which summary statistics are calculated varies, as do the metrics themselves (e.g., reporting mean versus median depth), while the selection of depth bins, particularly those used in older tags, can be highly variable across users tied with the species’ ecology. To facilitate meaningful inter- and intraspecies comparisons of elasmobranch vertical habitat use at a global scale, data harmonization is required across deployments and tag types.

In this study, we characterize standardized vertical habitat use patterns of elasmobranch species on a global scale. We synthesize depth data from 989 biotelemetry tags deployed across 38 species, encompassing eight orders (14 families), spanning those inhabiting coastal and oceanic habitats from polar to tropical latitudes, and including a wide range of foraging and movement modalities. Specifically, we (i) characterize the vertical distributions and key vertical movement metrics for each species and test whether these were maintained across broad geographic regions and unique ocean basins (i.e., assessment of intraspecific variation), (ii) quantify the level of interspecific overlap in vertical distributions, (iii) assess the relative occurrence of species’ diel vertical behaviors [nDVM (normal DVM), rDVM (reverse DVM), and neutral], and (iv) examine the intrinsic and extrinsic factors driving variation in vertical habitat use across global elasmobranch species. In doing so, we provide a valuable reference point for hypothesis generation to further our ecological knowledge of elasmobranchs as well as informing future monitoring and management efforts.

RESULTS

Data synthesis

We compiled 96,169 days of vertical movement data from 989 individuals, encompassing 38 species that were representatives of 14 families of elasmobranch. Most depth data were obtained from tag deployments on representatives of the families Lamnidae and Carcharhinidae, accounting for 43.7 and 21.7% of the data days and 37.4 and 28.9% of the individuals, respectively. The number of tagged individuals varied across species, ranging from 1 (Munk’s pygmy devil ray, Mobula munkiana; pelagic stingray, Pteroplatytrygon violacea; and Cuban dogfish, Squalus cubensis) to 187 (white shark, Carcharodon carcharias) (Table 1). The number of days of data recorded by the tags per individual ranged from 7 to 784, with a mean (± SD) of 98.4 ± 80.9 data days per tag deployment across all species. Collectively, archived and transmitted time-series data were available for 59.7% of tags (n = 590), with the remaining tags providing data in the form of summarized histograms. Habitat type, assigned on the basis of the presence of a species on the continental shelf as found in previous studies, was classified as “coastal” for 11 species, “transient” for 18 species, and “oceanic” for 9 species (table S1).

Table 1.

Vertical reference metrics for 38 elasmobranch species.

TS, time series; IQR, interquartile range. Mean skewness and percentage metrics are calculated from time-series data only. Maturity status (I, immature; M, mature; and U, unknown maturity status) was estimated on the basis of length measurements taken during tagging activities and assessed against published maturity lengths (table S1). Calculations of mean depth, median depth, and mean skewness were only available for species where time-series datasets were available. NA, not applicable.

Species	N (TS/total)	Maturity status I/M/U	Mean data days ± SD	Total data days	Depth range (min-max)	Mean max depth ± SD (m)	Mean depth ± SD (m)	Median depth	Mean skewness ± SD
								(IQR) (m)	(range)
Pelagic thresher shark Alopias pelagicus	6/6	6/0/0	76.3 ± 62.5	458	0–584	452 ± 159.9	115.3 ± 39.1	82 (42.4–177.3)	0.94 ± 0.23 (0.54–1.2)
Bigeye thresher shark Alopias superciliosus	3/3	0/1/2	48 ± 21.7	144	5–543	518.2 ± 24.3	219.9 ± 31	211.6 (84.3–342.5)	0.22 ± 0.29 (−0.07–0.51)
Common thresher shark Alopias vulpinus	0/8	1/2/5	68.5 ± 36.7	548	0–572	261 ± 191.6	NA	NA	NA
Arctic skate Amblyraja hyperborea	5/5	5/0/0	38.6 ± 20	193	317–1400	1306 ± 110.5	944.4 ± 155.8	952.3 (835.3–1061.5)	0.02 ± 0.71 (−0.7–0.81)
Big skate Beringraja binoculata	2/5	2/3/0	189 ± 106.3	945	1–500	377.8 ± 172.1	81.9 ± 62.2	81.8 (54.1–97.3)	0.62 ± 0.15 (0.52–0.73)
Silvertip shark Carcharhinus albimarginatus	7/11	8/2/1	93 ± 41.6	1023	0–792	441.4 ± 220.3	42.6 ± 3.8	41 (31.4–52)	2.05 ± 1.02 (0.52–3.64)
Grey reef shark Carcharhinus amblyrhynchos	5/5	0/5/0	25.2 ± 26.5	126	0–147	99.8 ± 28.5	27.5 ± 6.4	23.1 (13–42.1)	0.53 ± 0.36 (0–0.86)
Bronze whaler shark Carcharhinus brachyurus	7/7	1/6/0	126.1 ± 74.1	883	0–129	96.9 ± 26.3	25.9 ± 10	24 (10.5–37.6)	0.6 ± 0.82 (−0.17–2.16)
Silky shark Carcharhinus falciformis	32/37	32/2/3	51.2 ± 48.2	1893	0–1112	294.5 ± 187.2	41.2 ± 14.3	40 (22–57.5)	0.75 ± 0.8 (−0.73–2.35)
Bull shark Carcharhinus leucas	17/18	1/17/0	89.8 ± 67.4	1616	0–256	148.1 ± 66.5	24.7 ± 11.4	22.7 (12.9–33.7)	0.97 ± 0.93 (−0.55–3.17)
Blacktip shark Carcharhinus limbatus	5/10	0/10/0	63.4 ± 50.1	634	0–132	60.3 ± 47.9	11.1 ± 4.3	10.7 (5.5–14.5)	0.76 ± 1.1 (−0.35–2.22)
Oceanic whitetip shark Carcharhinus longimanus	19/22	12/5/5	88.6 ± 41.4	1950	0–659	253.2 ± 127.7	32.5 ± 6.8	25.4 (10.2–48.6)	1.34 ± 0.56 (0.49–2.98)
Galapagos shark Carcharhinus galapagensis	9/10	4/6/0	63.4 ± 40.7	634	0–528	317.4 ± 126.0	53.2 ± 7.2	53.2 (32.4–74.2)	0.36 ± 0.74 (−0.89–1.58)
Caribbean reef shark Carcharhinus perezi	10/10	0/10/0	118 ± 97.4	1180	0–697	311.4 ± 161.9	29.5 ± 11.7	26.9 (18.6–33.4)	4.17 ± 2.26 (0.15–8.06)
White shark Carcharodon carcharias	93/187	126/41/20	113.5 ± 80.2	21,220	0–1277	541.5 ± 339.5	48.9 ± 38.5	21.3 (6.4–61.9)	2.49 ± 2.23 (−4.43–10.01)
Basking shark Cetorhinus maximus	31/66	7/51/8	121.6 ± 77.9	8025	0–1504	782.6 ± 426.5	177.8 ± 140.1	123.9 (56.2–273.1)	1.68 ± 1.74 (−0.79–6.36)
Tiger shark Galeocerdo cuvier	46/55	17/38/0	54.9 ± 45.9	3019	0–1275	519.7 ± 241.7	40.1 ± 19.9	26.5 (10–53.5)	2.82 ± 1.75 (−0.17–8.85)
School shark Galeorhinus galeus	16/17	0/17/0	78.4 ± 47.3	1332	0–696	234.4 ± 226.0	68.5 ± 46	56.6 (38–86.8)	0.53 ± 1.03 (−1.32–2.82)
Bluntnose sixgill shark Hexanchus griseus	2/2	2/0/0	201 ± 0	402	159–904	891.8 ± 16.6	449 ± 1.7	528 (267.5–590.8)	−0.19 ± 0.01 (−0.2 - −0.18)
Lutz’s stingray Hypanus berthalutzae	1/2	0/2/0	17.5 ± 14.8	35	0–78	51.8 ± 36.4	3.1	1 (0–2.5)	4.94
Shortfin mako shark Isurus oxyrinchus	9/57	54/2/1	94.9 ± 44.5	5409	0–1888	481.1 ± 346.3	55.5 ± 33.9	42.4 (23.1–75.2)	2.34 ± 1.67 (0.88–5.3)
Longfin mako shark Isurus paucus	1/3	0/3/0	94.7 ± 30.9	284	0–1752	1482.7 ± 459.6	222.5	239 (77–340.5)	0.21
Salmon shark Lamna ditropis	11/59	11/48/0	117.2 ± 68.3	6917	0–968	411.6 ± 193.6	49.6 ± 32.6	31.9 (10.4–71.9)	2.19 ± 1.66 (0.73–6.27)
Porbeagle shark Lamna nasus	42/64	48/16/0	128.5 ± 78.7	8224	0–1305	601.7 ± 312.9	93.2 ± 73.7	66.7 (23–138.4)	1.71 ± 1.6 (−0.04–7.78)
Reef manta ray Mobula alfredi	35/64	35/13/16	76.9 ± 52	4923	0–711	275.4 ± 161.2	25.2 ± 11.5	20.6 (8.6–38.1)	1.3 ± 1.03 (−0.48–4.85)
Oceanic manta ray Mobula birostris	11/11	1/7/3	57.7 ± 22.5	635	0–1246	442.5 ± 325.3	38.2 ± 21.4	31.5 (14–54.8)	3.21 ± 3.47 (0.37–11.15)
Spinetail devil ray Mobula mobular	1/14	14/0/0	40.4 ± 26.2	565	0–476	340 ± 105.2	11.7	1 (0–7.5)	4.68
Munk’s pygmy devil ray Mobula munkiana	1/1	0/0/1	28	28	0–126	126	11.7	9 (5.5–15)	2.79
Sicklefin devil ray Mobula tarapacana	5/6	5/1/0	48 ± 40	288	0–1637	1208.1 ± 547.2	89.6 ± 11.3	59.4 (26.8–107.6)	3.96 ± 2.01 (1.48–6.8)
Starry smooth-hound Mustelus asterias	7/7	5/2/0	129.3 ± 79.6	905	0–118	83.5 ± 35.9	30.5 ± 9.8	28 (18.4–43)	0.26 ± 0.27 (−0.01–0.69)
Broadnose sevengill shark Notorynchus cepedianus	2/5	0/5/0	71.4 ± 51.1	357	0–222	194.8 ± 38.5	40.2 ± 7.3	26.5 (20.5–55.5)	1.52 ± 0.24 (1.35–1.68)
Blue shark Prionace glauca	52/101	10/86/5	77.9 ± 50.8	7871	0–1792	702.6 ± 411.9	86.5 ± 34.2	51.8 (14.4–125.1)	1.72 ± 0.85 (−0.05–4.32)
Common sawshark Pristiophorus cirratus	3/3	0/3/0	14.7 ± 7.5	44	5–121	105.4 ± 17.6	79.7 ± 10.1	81.7 (63.8–95.1)	−0.45 ± 0.27 (−0.73 - −0.21)
Pelagic stingray Pteroplatytrygon violacea	1/1	0/1/0	58	58	3–428	428	104	100.5 (68–129)	0.93
Whale shark Rhincodon typus	48/61	12/48/1	147.5 ± 161	8997	0–1896	1055.6 ± 537.2	40.8 ± 17	24 (6.8–51.1)	5.02 ± 3.3 (0.14–16.24)
Greenland shark Somniosus microcephalus	27/28	8/20/0	138.1 ± 126.7	3867	0–1547	969.2 ± 290.1	379.9 ± 155.8	379.7 (308.9–440.1)	0.66 ± 1.06 (−1.15–3.48)
Scalloped hammerhead Sphyrna lewini	16/17	0/17/0	30.8 ± 42.8	523	0–971	555.3 ± 263.8	58.9 ± 33.8	43.3 (26.3–64.5)	2.83 ± 1.77 (0.14–5.63)
Cuban dogfish Squalus cubensis	1/1	0/1/0	14	14	324–710	710.1	463	441.1 (418.3–485.5)	1.42

Tagging locations were distributed across 17 of the 31 marine biogeographic realms defined by Costello et al. () (Fig. 1). One-quarter (25.0%) of all tag deployments were in the North Pacific marine biogeographic realm, followed by the Tropical Western Atlantic (15.3%), and the Tropical Indo-Pacific and Coastal Indian Ocean (9.6%) (Fig. 1 and table S3). Tag detachment locations were available for 845 individuals and spanned 21 different marine biogeographic realms, with 70.3% of the tags detaching within the same marine biogeographic realm they were deployed in.

Fig. 1.

Deployment and pop-up and/or recapture locations of tracked elasmobranchs.

Yellow triangles indicate deployment and red circles indicate pop-up and/or recapture of the 989 elasmobranchs included within the analysis for this study. Numbers refer to the ocean biogeographic realms as defined by Costello et al. () (see table S3). Pop-up locations were not available for 144 tags.

Vertical habitat use and reference metrics

A whale shark, Rhincodon typus, recorded the deepest depth in our dataset with a maximum depth of 1896 m, which approached the physical limit of the pressure sensors of electronic tags. Depths greater than 200 and 1000 m were attained by 31 and 13 of the 38 species, respectively (Table 1). A total of 26 species spent over 95% of their time within the top 250 m of the water column (Fig. 2), while 4 species spent less than 50% of their time in the top 250 m (Fig. 2; table S4). Only 8 species spent more than 25% of their time in the top 5 m (a zone reflecting surface-oriented behaviors and where animals are potentially detected by aerial surveys), and 13 species spent less than 5% of their time within this zone (Fig. 2 and table S4). Coastal species spent the highest proportion of time in the top 100 m of the water column (98.9 ± 1.6%), followed by transient (70.5 ± 30.8%) and oceanic (54.4 ± 36.2%) species. The range of mean sea surface temperature (SST) (°C), the difference between the minimum and maximum mean SST at a location of deployment for each species, was narrower for coastal species (3.1 ± 2.6) compared to transient (5.2 ± 4.0) and oceanic (4.1 ± 5.1) species (fig. S2).

Fig. 2.

Comparison of epipelagic water occupation by each tagged elasmobranch species.

Mean percentage of time at liberty spent by tagged elasmobranchs within the (A) top 5, (B) top 50, (C) top 100, and (D) top 250 m of the water column. Error bars represent ±1 SD and are truncated at 0 and 100%. Exact values can be extracted from table S4, along with the mean percentage of time spent in the top 10 m. Species are sorted from top to bottom by lowest to highest use of the top 250 m to ease interpretability. Cuban dogfish and Arctic skate spent all their time deeper than 250 m.

Median depths varied among species, ranging from the shallow occupancy of 1 m [interquartile range (IQR) = 0 to 2.5 m] for Lutz’s stingray (n = 1 with time-series data) and 11 m (IQR = 5.5 to 14.5 m) for blacktip sharks, Carcharhinus limbatus (n = 5) to deeper waters occupied by Arctic skate, Amblyraja hyperborea (892.5 m; IQR = 835.3 to 1061.5 m; n = 5) (Table 1 and fig. S3). Distributions of depth data for 34 of the 37 species where time-series data were available were positively skewed (i.e., higher use of shallower depths with tails into deeper waters; Fig. 3, Table 1, and fig. S3). In contrast, depth data for the common sawshark, Pristiophorus cirratus, (skewness = −0.45 ± 0.27; n = 3) and bluntnose sixgill shark, Hexanchus griseus (−0.19 ± 0.01; n = 2), were negatively skewed toward deeper waters, whereas the depth distribution for the Arctic skate was symmetric (0.02 ± 0.71, n = 5) (Table 1).

Fig. 3.

Vertical distributions and diel behavior of 15 elasmobranch species.

The hourly median depth distributions of 15 elasmobranch species determined from hourly median depths from each satellite-tagged individual within each species. Only species with >1000 days of depth time-series data were incorporated into this figure (fig. S3 shows a corresponding figure with all available species). Violin plots represent the full distribution of the data, with colors relating to family. Boxplots depict the lower quartile, upper quartile (and thus the interquartile range), and median within the data, with whiskers extending from the shallowest to the deepest depth observed within each species. Whiskers are capped to 1200 m to improve visual interpretation, with the maximum depths of species that exceed this threshold stated at the bottom of the whisker. Bars represent the estimated detection zones of aerial surveys (top 5 m; drone icon), scuba-diving surveys (top 50 m; diver icon), and longline fishing (top 250 m; fish and hook icon) used within this study. Pie charts represent the proportion of individuals within each species that primarily exhibited nDVM, rDVM, or no clear evidence of DVM (neutral) as determined by nonparametric Wilcoxon signed-rank tests applied to time-series data. Species are ordered by habitat type, moving from oceanic to transient to coastal species from left to right.

Regional variability

Tag data were available from more than one Food and Agriculture Organization (FAO) region for 18 species, revealing patterns of intraspecific variation in depth use among regions (figs. S6 to S23). Six of the 12 species with enough individuals tagged among regions to enable statistical comparison displayed significant variation in depth use (differences among FAO areas were significant for >2 metrics). Vertical distribution patterns of the other six species were consistent among regions (figs. S6 to S23).

Interspecific overlap in use of vertical habitats

Coefficients of similarity for time-series depth data binned and averaged at 10-m intervals for 28 elasmobranch species had a mean of 0.77 ± 0.27, ranging from 0 (between five coastal species and Arctic skates) to 0.996 (between whale sharks and tiger sharks, Galeocerdo cuvier) (Fig. 4). Hierarchical clustering based on dissimilarity revealed four discrete clusters (Fig. 5 and fig. S4). Cluster 1 contained seven of the eight coastal species and the oceanic whitetip shark, Carcharhinus longimanus, with six of these species being of the family Carcharhinidae. Species in this cluster had relatively high occupancy of the top 50 m of the water column, and very little use of depths >100 m (Fig. 5). Cluster 2 contained more than half of the total species included in this analysis (n = 16) and was distinguished by depth distributions spanning surface waters to depths >100 m (Fig. 5 and fig. S4). Cluster 2 could be further broken up into three subclusters. The first contained the silvertip shark, Carcharhinus albimarginatus, Galapagos shark, Carcharhinus galapagensis, scalloped hammerhead shark, Sphyrna lewini, and silky shark, Carcharhinus falciformis, and was characterized by a relatively even depth distribution across the top 100 m. The second subcluster contained species with the most similar depth distributions, with a mean coefficient of similarity of 0.97 ± 0.02, and included white sharks, oceanic manta rays, Mobula birostris, tiger sharks, whale sharks, salmon sharks, Lamna ditropis, and shortfin mako sharks, Isurus oxyrinchus. These species had depth distributions that were strongly right skewed, with high occupation of the top 10 m. The final subcluster contained the basking shark, Cetorhinus maximus, sicklefin devil ray, Mobula tarapacana, pelagic thresher shark, Alopias pelagicus, school shark, Galeorhinus galeus, blue shark, Prionace glauca, and porbeagle shark, Lamna nasus, and was distinguished by having relatively higher occupation of deeper depths (Fig. 5 and fig. S4). Clusters 3 and 4 consisted of the Arctic skate and Greenland shark, Somniosus microcephalus, respectively, that resided in deeper waters (Fig. 5 and fig. S3).

Fig. 4.

Vertical habitat overlap between elasmobranch species.

Matrix of vertical habitat overlap (Bhattacharyya coefficient) among species, where zero indicates no overlap between depth distributions and one indicates identical depth distributions. Calculations were based on time-series depth data binned at 10-m intervals for each individual and averaged across a species. Only species with five or more individual depth time-series datasets were incorporated into this analysis (n = 26).

Fig. 5.

Clustered depth distributions of elasmobranchs.

Depth distributions for 26 elasmobranch species with time-series data binned at 10-m intervals. Note that the plot has been limited to the top 150 m to ease interpretation but extends to 1850 m (see fig. S4 for the full plot). Italicized lettering next to each species name indicates the habitat type of each species (c = coastal, t = transient, and o = oceanic). The dendrogram and clusters on the right side of the figure resulted from hierarchical cluster analysis performed on dissimilarity of the Bhattacharyya coefficient. Numbered clusters represent species grouped according to similarity in vertical habitat use.

Diel variation in vertical distributions

Time series of depth data enabled daytime (10:00 to 14:00) and nighttime (22:00 to 02:00) median depths to be compared for 577 tagged individuals of 37 species (Fig. 3 and Table 2). nDVM (diving deeper during the day than night) was the most common strategy across individuals and species (found in 46% of individuals and 89% of species), followed by neutral patterns (no difference between day and night; 41% of individuals and 73% of species) and rDVM (deeper during the night than day; 13% of individuals and 38% of species). However, diel patterns of movement were not consistent within species. Twelve species had at least one individual present in each of these diel movement categories (Table 2), and three species did not have any individuals exhibiting patterns of nDVM or rDVM [Arctic skate (n = 5), blacktip shark (n = 5) and spinetail devil ray, Mobula mobular (n = 1)]. Of the 16 species with 10 or more individuals tracked, 10 (63%) displayed all three behavioral modes (Table 2). Intraspecific variability in diel patterns tended to be more common in oceanic and transient species than in coastal species. The reef manta ray, Mobula alfredi, was the only coastal species to display all three patterns (nDVM = 36.1%, rDVM = 16.7%, neutral = 47.2%, n = 36), compared to 58 and 50% of transient and oceanic species, respectively.

Table 2.

DVM patterns for elasmobranch species.

Only species with time-series data available were used, with data split into local day (10:00 to 14:00) and night (22:00 to 02:00) periods. Percentage values represent the proportion of individuals of a species displaying nDVM, rDVM, or neutral patterns (no difference between day and night). Day and night counts display the number of days of data available for each respective diel period.

Species	Number of tags	% nDVM	% rDVM	% Neutral	Day count	Night count
Pelagic thresher shark Alopias pelagicus	5	100	0	0	368	422
Bigeye thresher shark Alopias superciliosus	3	100	0	0	144	144
Arctic skate Amblyraja hyperborea	5	0	0	100	143	133
Big skate Beringraja binoculata	2	50	0	50	387	407
Silvertip shark Carcharhinus albimarginatus	7	71.4	0	28.6	510	516
Grey reef shark Carcharhinus amblyrhynchos	3	33.3	33.3	33.3	92	92
Bronze whaler shark Carcharhinus brachyurus	4	75	0	25	531	526
White shark Carcharodon carcharias	91	38.5	12.1	49.5	11,612	11,527
Silky shark Carcharhinus falciformis	32	50	12.5	37.5	1323	1282
Galapagos shark Carcharhinus galapagensis	9	66.7	0	33.3	443	459
Bull shark Carcharhinus leucas	16	62.5	0	37.5	1415	1432
Blacktip shark Carcharhinus limbatus	5	0	0	100	457	452
Oceanic whitetip shark Carcharhinus longimanus	19	26.3	31.6	42.1	1676	1674
Caribbean reef shark Carcharhinus perezi	10	80	20	0	1103	1103
Basking shark Cetorhinus maximus	30	26.7	13.3	60	3960	3760
Tiger shark Galeocerdo cuvier	45	44.4	6.7	48.9	2069	2046
School shark Galeorhinus galeus	14	64.3	0	35.7	1016	967
Bluntnose sixgill shark Hexanchus griseus	2	100	0	0	400	398
Lutz’s stingray Hypanus berthalutzae	1	0	100	0	13	11
Shortfin mako shark Isurus oxyrinchus	9	55.6	22.2	22.2	1034	1032
Longfin mako shark Isurus paucus	1	100	0	0	115	94
Salmon shark Lamna ditropis	10	70	0	30	1653	1655
Porbeagle shark Lamna nasus	42	81	0	19	5066	5074
Broadnose sevengill shark Notorynchus cepedianus	2	50	0	50	202	201
Reef manta ray Mobula alfredi	36	36.1	16.7	47.2	2273	2246
Oceanic manta ray Mobula birostris	11	9.1	54.5	36.4	510	582
Spinetail devil ray Mobula mobular	1	0	0	100	67	68
Munk’s pygmy devil ray Mobula munkiana	1	100	0	0	27	27
Sicklefin devil ray Mobula tarapacana	5	100	0	0	189	177
Starry smooth-hound Mustelus asterias	7	71.4	0	28.6	900	905
Blue shark Prionace glauca	54	59.3	7.4	33.3	4176	4132
Common sawshark Pristophorus cirratus	3	33.3	0	66.7	40	44
Pelagic stingray Pteroplatytrygon violacea	1	100	0	0	47	50
Whale shark Rhincodon typus	48	12.5	45.8	41.7	7902	7922
Greenland shark Somniosus microcephalus	27	37	0	63	3380	2886
Scalloped hammerhead Sphyrna lewini	15	20	26.7	53.3	441	456
Cuban dogfish Squalus cubensis	1	100	0	0	13	13

Drivers of elasmobranch depth preferences

Our best model explaining variation in elasmobranch median depths (see table S5 for the full list of candidate models) indicated that spatial and species terms were strong predictors of median depth [Best Integrated Nested Laplace Approximation (INLA) regression model, Table 3]. Interaction terms between species and realm and sex had strong effects on median depths, whereas the interaction between species and SST had a weak positive effect. The fixed effects of sex, SST, and trophic level alone had no effect on median depth. Median depth also decreased with maximum size across species (posterior mean, −0.21; CrI, −0.32, −0.03) and for species primarily associated with coastal habitats (posterior mean, − 0.21; CrI, −0.38, −0.04) (Table 3). We found no evidence of an effect of maturity status on median depth, even when included as an interaction term with species.

Table 3.

Drivers of elasmobranch depth preferences.

Mean of the posterior coefficient estimates and the quantiles of the posterior distribution (0.025 and 0.975) from the best Bayesian regression model examining the median depths of tagged elasmobranchs from 38 species within the first 7 days of tracking (excluding the first day of deployment) using INLA. Trophic level, maximum size, and primary habitat for each species were determined from existing literature (table S1). The maturity status for each individual was determined from length measurements taken when the animal was tagged and compared to published lengths at maturity (table S1). SST at the tagging location obtained from the National Oceanic and Atmospheric Administration’s Multi-Scale Ultra-High Resolution level 4 analysis on a 0.01° spatial resolution and averaged across the 7 days following deployment. Spatial terms are derived from the latitude and longitude of the tagging location, and the realm is the biogeographic realm () where the tagging occurred. Terms for which the 95% credible interval of the posterior does not overlap zero are italicized. Asterisk denotes the interaction term.

Variable	Mean	0.025 quantile	0.975 quantile
Intercept	1.91	0.29	3.42
Habitat: coastal	-0.21	−0.38	−0.04
Habitat: coastal transient	0.16	0.00	0.32
Habitat: oceanic	0.05	−0.13	0.22
Maturity: immature	−0.02	−0.08	0.04
Maturity: mature	0.03	−0.03	0.09
Maturity: unknown	0.00	−0.10	0.09
Max size	−0.17	−0.32	−0.03
Sex: F	−0.04	−0.12	0.05
Sex: M	0.02	−0.06	0.11
Sex: U	0.01	−0.12	0.15
SST	0.00	−0.01	0.01
Trophic level	0.00	−0.34	0.38
Range for space	18.80	8.20	37.30
Stdev for space	0.17	0.11	0.24
Precision for species	0.02	0.20	0.01
Precision for species*realm	0.14	0.27	0.09
Precision for sex*species	0.10	0.20	0.06
Precision for species*SST	0.01	0.02	0.01

DISCUSSION

The vertical dimension of movement is required to gain a complete understanding of aquatic species’ three-dimensional space use to inform ecology and management. Compiled vertical movement data for elasmobranchs derived from telemetry devices revealed high intra- and interspecific variability in vertical movement patterns, across both geographical and diel scales, highlighting the complex behaviors adopted by this diverse taxonomic group. Yet, substantial vertical overlap occurred between many epipelagic elasmobranchs, and vertical habitat use of certain threatened species was consistent across global oceans, indicating convergent strategies and potential for standardized management practices to mitigate bycatch and enact effective conservation measures. We discuss the ecological and applied implications of observed patterns and identify future research priorities.

Ecological insights

Vertical habitat use and reference metrics

The majority of the 38 species we examined had epipelagic distributions (n = 31), with median depths in the top 200 m of the water column. However, many of these species frequented deeper bathomes, with 31 recorded at mesopelagic depths (>200 m) and 13 at bathypelagic depths (>1000 m). The motivation for irregular movements into deep, dark, cold waters for epipelagic species is difficult to determine but could be related to some combination of foraging opportunities, thermoregulation, bioenergetics, reproduction, predator avoidance, and/or navigation through magnetic, chemical, topographic, electric, and light cues (). Entering deep waters for sustained periods, however, requires a capacity to endure cooler temperatures and potentially more hypoxic conditions that are likely to be outside the “normal” environmental niche of many epipelagic species. Seven of the 10 epipelagic species found to access depths >200 m have either morphological [i.e., large size and therefore thermal inertia; ()] and/or physiological adaptations [i.e., regional endothermy; ()] that slow the rate of body temperature cooling at depth to facilitate access to deeper maximum depths. Individuals may also use behavioral adaptations to cope with these unfavorable conditions, for instance, by bounce diving between shallow and deep waters to thermoregulate and/or reoxygenate ().

Regional variability

Several species, including oceanic whitetip, tiger, scalloped hammerhead, school, and silky sharks, though capable of undertaking deep dives to meso- and bathypelagic depths (, ), recorded consistent use of relatively shallow habitats across biogeographic realms. This uniform behavior identifies that generalized fisheries management approaches could be adopted among FAO regions, particularly where regional management is hindered by data limitations for a given species. However, high intraspecific variability in vertical movement patterns of species both within and between regions was also observed. While regional-scale studies have also recorded such variation (), our global-scale study provides strong evidence that this phenomenon is common across species and regions. Drivers of intraspecific differences can be difficult to unravel without fine-scale, species-level analyses but most likely relate to the combined effects of local prey composition, abundance and distribution, discrete oceanographic conditions that dictate water column structure, and bathymetry.

Interspecific overlap in use of vertical habitats

Vertical habitat overlap among species analyzed was generally high, potentially reflecting the bias toward tracking studies of more accessible, epipelagic species, as well as the current pressure limitations of available tagging technologies. High overlap in use of vertical habitats among species indicates the increased likelihood of ecological interactions, such as intraguild predation and competition. The cluster analysis revealed high overlap between more closely related species and/or species inhabiting similar habitat types (i.e., coastal or oceanic habitats). Depth distributions primarily limited to the top 50 m were recorded for species in cluster 1, reflecting their coastal distribution and therefore bathymetric limitations, with the notable exclusion of the oceanic whitetip shark, which inhabits the open ocean. Cluster 2 grouped 16 species with relatively broad epipelagic distributions that was then broken into three distinct subclusters. The first of these contained three carcharhinids and the scalloped hammerhead shark, all of which are primarily piscivorous species found in warm temperate to tropical waters. The second contained three lamnids (white shark, salmon shark, and shortfin mako shark), the tiger shark, and the filter-feeding whale shark and oceanic manta ray, characterized by high surface use with occasional excursions into deeper waters. The latter three species contained the most similar vertical distributions of our analyses. Notably, these species also share other comparable movement behaviors, including distributions dominated by residency in coastal sites in warm temperate and tropical regions tied with preferential prey (i.e., zooplankton and sea turtles, respectively), where bathymetry may limit vertical movements, and offshore migrations through deeper waters with records of deep-diving behavior (–). In the case of the filter-feeding oceanic manta ray and whale shark, spatial overlap may result in competition for prey, whereas overlap of these taxa with tiger sharks may increase the likelihood of predator-prey interactions (). The final subcluster recorded relatively deeper distributions at epi- and mesopelagic depths, grouping oceanic and transient species, all of which (excluding the relatively understudied pelagic thresher shark) have been tracked on long-distance movements and/or migrations [e.g., blue shark (), porbeagle shark (), school shark (), and sicklefin devil ray ()]. These broad distributions could require changes in vertical behaviors with migration phase and/or horizontal position, as has been previously recorded for latitudinal migrations in basking sharks () and blue sharks (). Although overlap in vertical habitats does not necessarily indicate that species are simultaneously present in space and time, similarities in habitat and horizontal distributions, and evidence for intraguild predation from diet studies () for species that clustered together in our analyses, support spatial and temporal overlap. Future work could consider range distributions and both horizontal and vertical movements concurrently and over time to quantify the full extent of species overlap for elasmobranchs.

Diel variation in vertical distributions

Diel changes in vertical patterns were exhibited with nDVM the dominant behavior. Prevalence of nDVM patterns is expected, given the daily migration of global zooplankton, forage fish, and associated predators toward the surface at dusk and toward the mesopelagic at dawn driven by avoidance of visual predators in well-lit surface waters (). While elasmobranch dive behavior is commonly linked to foraging (), their time at depth can be limited by exposure to unfavorable physical conditions, such as reduced temperatures and dissolved oxygen concentrations, resulting in deviations to the model DVM pattern. Spatial and temporal plasticity in DVM patterns by zooplankton and other prey species as a result of “landscapes of fear” and lunar illumination can also modify systematic nDVM ().

Many elasmobranch species showed high intraspecific variation and plasticity in DVMs, displaying various combinations of nDVM, rDVM, and neutral patterns. Variation in DVM can result from prevailing environmental conditions (), bathymetry, spatial and temporal variation in prey distribution (), and/or sample size limitations (). Reef manta rays, for example, displayed nDVM patterns around the Chagos Archipelago (), rDVM around the Seychelles (), and both patterns off eastern Australia (), presumably due to site-specific vertical distributions of zooplanktonic prey and baseline productivity. Individuals may also display both patterns of rDVM and nDVM tied with habitat type, as observed for basking sharks moving between shallow, inner shelf habitats and deep, shelf-edge habitats (). Evidence for intraindividual flexibility in DVM is potentially favorable in the context of climate change, where plankton distributions are predicted to shift unpredictably with prevailing environmental conditions (). The depth behavior of each individual was summarized across its full time at liberty, potentially masking discrete shifts in DVM throughout the tag deployment tied with specific locations.

Drivers of vertical distributions

Fine-scale vertical behaviors of elasmobranchs are thought to be driven by a need to optimize the energy costs of locomotion and foraging, while remaining within physiological limits imposed by abiotic factors, such as ambient water temperature and dissolved oxygen levels (, ). It was not possible to incorporate these fine-scale variables into our analysis given data aggregation; however, key behavioral drivers emerged at the spatial-temporal range examined.

As would be expected, tagging location was a significant predictor of median depth, a result of the heterogeneous environmental, oceanographic, and bathymetric conditions that drive the distribution of fishes () and are tied to the specific region of capture (i.e., shelf or open ocean). Our best model, however, included an interaction between species and marine biogeographic realm, indicating that the unique environmental conditions that influence the horizontal distributions of species and ultimately underpin the classification of realms () also drive variability in vertical habitat use among elasmobranchs.

While interspecies variation was a strong predictor of median depth for the species examined, phylogeny contributed little to explain variation in our models, raising the question over whether vertical ecology could be reliably predicted across the breadth of the elasmobranch phylogenetic tree. Limited evidence of relationships between phylogeny and functional diversity within sharks was reported by Cachera and Le Loc’h () in agreement with our data for a subset of species. It is possible that distantly related species could have independently evolved similar patterns of depth use, driven by similar habitat and environmental conditions at discrete locations, such as those found in the deep ocean (). Given that our global data were limited to only 38 of the 1143 recorded elasmobranch species and 15 of the 57 extant families (), phylogenetic patterns cannot be ruled out and warrant further investigation. Furthermore, taxonomic, demographic, and geographic biases inherent in these data highlight the need for focused studies on underrepresented elasmobranch taxa. This is vital for effective conservation and management of elasmobranchs, given that more than one-third of all sharks and rays are threatened with extinction, according to the IUCN (International Union for Conservation of Nature) Red List of Threatened Species ().

Trophic level and maximum size were both included in our best model, but only maximum size was a significant driver of median depth, with larger species residing in shallower waters. This may be explained by the fact that both large zooplanktivores (i.e., whale sharks, basking sharks, and manta rays) and large higher-order predatory species (i.e., white sharks and tiger sharks) were tagged during their surface-based foraging phase [see, e.g., ()]. Models predicted that median depths were shallower for species associated with coasts than transient or oceanic taxa, which contain those feeding across a broad range of trophic levels and of various size classes. This is expected, given the greater habitat depth available to oceanic species and the oligotrophic environment that they inhabit compared to increased productivity in shallow water around islands and coastlines (). Previous studies have also shown marked ontogenetic shifts in horizontal habitat occupied [e.g., silky sharks in the Gulf of Mexico () and whale sharks in the northeast Pacific ()], and species-specific shifts in vertical habitat use with age (). We found no evidence, however, for a relationship between median depth and maturity stage, even when maturity was incorporated as an interaction term with species (which was not included in our best model). This lack of relationship highlights diverse strategies across the species studied, but we acknowledge the inherent limitation of including only certain size classes for species within the analysis, likely constrained by tag size (i.e., ability to tag juveniles or small bodied species) and ease of capture (i.e., cryptic habitat use of certain life stages).

Unexpectedly, we found limited evidence that SST at the tagging location was a significant driver of median depth, even when considered as an interaction with species. Previous single-species studies have reported temperature as a strong regulator of distribution, especially for ectotherms with individuals moving into deeper, cooler waters when surface waters were warmer [e.g., basking sharks () and oceanic whitetip sharks ()]. It is possible that the remotely sensed measure of SST used in the models was too coarse to describe the fine-scale effects on individual vertical movements. In addition, the effect of SST is modulated by subsurface physical variables, such as mixed-layer depth (), which can be highly variable but could not be estimated here. Aligning vertical data with archived temperature records provides an alternative avenue to assess how elasmobranch movement in this plane will be influenced by changing climates.

The mating habits of most elasmobranchs remain poorly characterized, yet many species are known to sexually segregate on horizontal axes because of habitat and prey preferences, and for females to reduce negative interactions during pregnancy and enhance offspring survival (, ). Yet, elasmobranchs are also known to sexually segregate in the vertical plane [i.e., white sharks () and flapper skate, Dipturus intermedius ()]. While sex on its own was not important, the interaction between sex and species had an effect within our model. However, while this interaction effect improved model fit, no species-specific effects were confidently (95% significance level) estimated as being either positive or negative. It is unclear whether this is a result of small effect sizes or data limitations for each of the species, or that our metric of central tendency masked variation in habitat use, particularly at fine temporal and spatial scales. Given the importance of understanding reproduction dynamics to conservation and management, more detailed investigations in this area should be a priority for future research.

Implications for management and monitoring

Assessing the efficacy of methods for detecting elasmobranchs (i.e., to assess biodiversity or abundance) and to mitigate against the likelihood of individuals encountering major threats (i.e., bycatch in fisheries) is key to ensure that approaches are effective and relevant for management. As a first step to incorporating vertical habitat use into these assessments, we generated the vertical metrics calculated for each species and estimated the proportion of time elasmobranchs spent within depth bins that overlapped with common detection methods and fishing efforts using pelagic longlines, providing an important first step for incorporating vertical habitat use into these assessments. In fisheries management, metrics of vertical habitat use can inform estimates of catchability to improve key fishery-derived abundance indices (), bycatch mitigation strategies (), and protected area designation and management (, ). Our estimates of vertical overlap with longline fishing are based on published hook depths (). However, they should be interpreted as indicative, rather than conclusive, given the high intraspecific variability in vertical movements displayed by many species, variation in hook depth of pelagic longlines, and catch rates of elasmobranchs that depend on the target species, season, habitat being fished (), and vertical movements of animals attracted to bait.

Monitoring techniques

The proportion of time that a species spends near the surface dictates suitability for several techniques used to sample abundance. For example, aerial surveys are increasingly used in elasmobranch research to monitor occupancy and population size () but are only effective for species that spend a significant proportion of their time in the top few meters of the water column during daylight hours. This metric will therefore dictate confidence in estimated abundance values. Similar considerations of species depth use should be applied to each sampling approach on a case-by-case basis (i.e., detection bias), acknowledging that extrapolations on the estimates provided here should be undertaken with caution, as the maximum depth that an animal can be detected is likely spatially and temporally variable and influenced by prevailing ocean conditions, overhead weather conditions, turbidity, the size and shape of the animal, and the underlying habitat (e.g., sandy bottom versus deep water) (–). The proportion of time that a species spends at the surface is also important for determining the utility of location-transmitting tags, such as Smart Position or Temperature Transmitting tags, in tracking the horizontal positions of a species. Consequently, while the thresholds presented here and the proportions of time spent within each band by each species should be considered as conservative, they serve as useful baselines for future monitoring activities, when accounting for local conditions and context.

Generalized management approaches

The identification of consistent patterns of vertical habitat use, both intra- and interspecifically, is important from a management perspective, when considering generalized management approaches for data-deficient regions. For example, species that displayed consistent patterns of DVM among individuals (e.g., porbeagle sharks) would benefit from management strategies where fishing gears are only deployed during one diel period (i.e., day or night), or at target depths modified over a diel cycle, to reduce interaction rates (). Similarly, for species that recorded similar vertical distributions among FAO areas (e.g., the “Critically Endangered” oceanic whitetip shark), depth-tied management approaches could be consistently applied among regions as a first step to proactive conservation. High overlap in the depth distributions among species indicates a similar risk of exposure to threats that vary on a vertical gradient (e.g., fishing gears), identifying that coordinated management approaches will be effective, especially in the short term while data acquisition is ongoing. The next crucial step to inform spatial and temporal management strategies will be three-dimensional overlap analysis of elasmobranch species and their prey with fisheries and considering dynamic fisheries management approaches ().

Future directions

Through a global synthesis of elasmobranch vertical movement data collected by biologging and telemetry tags, it was possible to explore underlying relationships, commonalities, and variation among species and highlight opportunities for improved conservation and management. Yet, many questions remain. The investigation of the fine-scale physical and ecological drivers of vertical distributions at the assemblage level in particular should be prioritized. Given that temperature data are routinely recorded concurrently with pressure (i.e., depth), a pertinent next step would be to analyze the thermal niches of these species in relation to their vertical distribution and investigate potential thermal limits and thermoregulatory behaviors that may be driving observed vertical movement patterns among species with and without regional endothermy (). Environmental modeling could also be used to investigate how reduced oxygen availability may physically limit vertical movements () to predict climate change effects. Fortunately, continued technological advancement and innovation in the field promises to provide the tools needed to address these questions. Improvements to the diversity and transmission capability of sensors in biologging and biotelemetry tags () are also allowing us to contextualize further the vertical movement patterns of elasmobranchs, enabling comparisons of the energetics and foraging ecology among species and the investigation of the fine-scale mechanisms and processes underlying movement behaviors.

A wider range of species targeted by tagging studies is required to address spatial and taxonomic biases present in current global data. The creation of more affordable electronic tags will improve accessibility and equity in data-poor developing nations, whereas tag miniaturization will facilitate deployment on smaller-bodied species and younger-size classes that are currently underrepresented. Last, the development of tags that can withstand pressures >2000 m will be essential in extending movement studies to deep-water sharks that remain among the most understudied elasmobranch taxa.

MATERIALS AND METHODS

Data collection

Depth data were collated from either archived time series from recovered archival tags or transmitted time series and/or time-at-depth binned histogram data from satellite tags, for 1038 individuals across 38 species of elasmobranch tagged between 2000 and 2019. Deployment locations and, where possible, end-of-track positions were also reported for each tag deployment. All data were quality- and suitability-checked before inclusion in the study (see the “Data processing” section below). All animal handling procedures were approved by the respective ethical review committees at the data owners’ institutions. Elasmobranchs were tagged throughout the Arctic, Atlantic, Indian, and Pacific Oceans (Fig. 1), and tags were deployed on animals captured using either baited lines, during commercial fishing operations (longlining or purse seine) or while the animals were free-swimming.

For each species with corresponding tag data, we synthesized information on principal biological traits. Estimates of trophic level, maximum reported size, primary habitat type occupied for each species, and maturity status of each individual were compiled from the published literature (table S1). Primary habitat type was based on the presence of species in continental shelf habitats (), with species distributed exclusively on the continental shelf assigned as coastal, those off the continental shelf as oceanic, and those moving between the two and/or resident in both as transient (see table S1 for classifications and associated citations). Sex-specific maturity status (“mature” or “immature”) was estimated on the basis of length measurements taken during tagging; all length measurements were converted to total length if an alternate measurement was provided (i.e., fork length) via established length-length conversions (table S2) before being referenced against published length-at-maturity data (table S1). Where length data were available, but the sex was not specified, we conservatively assessed maturity status on the basis that the individual was female (i.e., longer to mature). Where length data were lacking (n = 53), maturity was defined as “unknown.” All data processing, analysis, and visualization took place in the R Statistical Environment.

Data processing

Data from tagged individuals providing less than 7 days of depth data were excluded. Data were trimmed to the time at liberty by identifying the point of initial submersion (i.e., tag deployment), and the point at which a constant depth was recorded for >24 hours (i.e., tag detached from animal and floating at the surface and/or mortality whereby a constant depth was observed). The first day of data (24 hours) was removed from each dataset to reduce biases due to potential atypical behaviors associated with capture and tag deployment (). Each tag dataset was also checked for drifting pressure sensors (i.e., a consistent and unnatural shift in the minimum depths being recorded), and where drift was apparent, the dataset was either cut to the period before drift [n = 5; e.g., ()] or removed from the database if the remaining period was less than 7 days (n = 1), resulting in a final dataset of 989 usable tags. Data were also queried for values >0 m [i.e., above the surface, which can occur in instances of exposure of tag at the air-water interface such as when an individual breaches ()] and where present, converted to zero values. To account for the influence of biogeography on vertical habitat use within and among species, we assigned marine biogeographic realms [as defined in ()] to each individual on the basis of the start and end points of the tag deployment (see table S3 for breakdown by species). Beyond latitude and longitude, these realm assignments account for influence of environmental heterogeneity and land barriers on species distributions and movements (). Additional positional information for individual sharks was not included because of the variable accuracy of light-based geolocation [error can be in the order of hundreds of kilometers; ()] and the frequency of geolocation data collection among datasets, which varies with tag type and time spent in the photic zone (). SST (°C) at the location of tag deployment was obtained from the National Oceanic and Atmospheric Administration’s Multi-Scale Ultra-High Resolution level 4 analysis on a 0.01° spatial resolution and averaged across the 7 days following deployment to correspond with INLA models (see the “Data analysis” section below).

Data analysis

For each species, we calculated the mean and total data days (i.e., tracking days with vertical data available, this may not represent the total days of a tag deployment because of transmission gaps in uplinks to satellites) and the absolute maximum and mean maximum depths across all individuals. Each dataset was summarized into the percentage of time that a species occupied a depth bin. For time-series data, bin sizes were set at the top 5, 10, 50, 100, and 250 m. Where possible, these depth bins were also used to summarize histogram data, combining existing bins where necessary (e.g., combining 0- to 5-m and 5- to 10-m bins into a single 0- to 10-m bin). Where this was not possible for a particular bin size, a tag was excluded from the respective calculation. Bin sizes were selected to coarsely correspond with the depth ranges represented in common monitoring methods and those widely reported in the literature. We therefore considered time spent within the following bands: (i) top 5 m of the water column for surface-oriented behavior and potentially detectable by aerial surveys (depending on turbidity) (); (ii) the top 10 m (); (iii) the top 50 m, which is commonly reported and represents a reasonable maximum detectable depth (considering dive depth and downward visibility) by scientific scuba surveys () and where the majority of BRUV deployments occur (); (iv) the top 100 m; and (v) the top 250 m to be representative of pelagic longline set depths ().

For individuals with time-series data available, we calculated the mean depth, mean skewness, and median depth and then averaged these for each species and determined the associated SD and/or IQR. This approach weighted each individual equally (i.e., did not bias to individuals with longer tracking durations). Given relatively high skewness values in many species, median rather than mean depth values were used in further visualization and analysis. To visualize core vertical habitat use, hourly median depths were calculated for each individual. These were then combined for each species and overall (grand) medians and IQRs were calculated. The full distributions of the data were visualized using violin plots for (i) species with >1000 days of combined time-series data available and (ii) all species with time-series data.

We also investigated intraspecific variation in vertical movements across the major FAO fishing areas. For species tagged in multiple FAO areas, we visualized vertical metrics (median depth, maximum depth, % top 5 m, % top 100 m, and % top 250 m) by deployment area, and statistical comparisons were made when sample sizes within FAO areas consisted of at least five individuals of a given species. Nonparametric Kruskal-Wallis tests were used to test for differences among FAO areas, given the data heteroscedasticity, with significance accepted at P < 0.05.

Overlap metrics and cluster analysis were used to compare the global vertical habitat use between species with five or more individuals with time-series datasets (n = 28 species). Specifically, depth data were binned into 10-m intervals for each individual and averaged, summarizing the proportion of time spent in each depth bin to a maximum of 1850 m (the maximum bin observed in time-series data). For each species, these proportions were averaged across individuals to obtain one distribution of proportion time at depth. Then, for each pair of species, the Bhattacharyya coefficient was computed using the proportion of time at depth. This coefficient ranges from 0 to 1 and measures the similarity between two discrete probability distributions, where 0 indicates no overlap between distributions and 1 indicates identical depth distributions (). Last, hierarchical clustering analysis was performed, using the “hclust” function in R, on a dissimilarity measure (1 − Bhattacharyya coefficient) using Ward’s minimum variance method to group species with relatively high overlap in vertical distributions into clusters. The optimum number of clusters was chosen by plotting the within-cluster sum of squares against the number of clusters (fig. S1).

Patterns of DVM (i.e., where vertical movements change across a daily cycle) are ubiquitous across elasmobranch species (, ). Variation in quantitative definitions, however, has precluded quantitative comparisons among species. To compare standardized diel variation in vertical distributions, we calculated the median daytime and nighttime depths for each individual from available time-series data. To ensure that estimated medians were not biased by crepuscular activity and/or movements through regions of varying day length, daytime depths were calculated from data between 10:00 and 14:00 local time and for nighttime between 22:00 and 02:00 local time. We note that, for all Arctic skates and some Greenland sharks, deployment of tags in the northern hemisphere summer resulted in almost constant daylight for tagged individuals. As data were often skewed, two-sample Wilcoxon tests with significance set at P < 0.05 were used to compare median daytime and nighttime depth [as per ()] within species and across biogeographic realms for each species. Movements were designated as nDVM, where depth was significantly deeper during the day than at night, rDVM, where depth was significantly deeper during the night than the day, and neutral, where no significant difference was found between daytime and nighttime depths.

Drivers of elasmobranch depth preferences

To investigate drivers of depth preferences of elasmobranchs, we fitted hierarchical Bayesian models with INLA using the “R-INLA” package (). The INLA methods determine an approximation to the posterior marginal distribution for parameters, which affords a large reduction in computational time compared to simulation-based (e.g., Markov chain Monte Carlo) methods ().

Within our initial full model, our response variable—log(median depth + 1)—was fitted within a normal likelihood. After removing data from the initial 24 hours after deployment, the next 7 days of track data per individual (or available data during this set period if data gaps occurred because of nonreceived transmissions) were used to ensure that the calculated median depths were closely tied to the spatial locations (i.e., the tagging location). This accepts varying levels of mobility among species, i.e., swim speeds that could lead to variable straight-line distances traveled over this set time period from the tagging location. We included an intercept term and a suite of fixed effects (“trophic level,” “maturity status,” “mean SST,” “max species size,” “primary habitat,” “marine biogeographic realm,” and “sex”) that were assigned conservative priors with a mean of 0 and SD of 10, which, given the scale of the continuous covariates, were sufficiently wide to allow these covariates to explain the full range of the response data. The categorical variables (“maturity status” and “primary habitat”) were given a sum-to-zero constraint to avoid identifiability issues with the intercept. Therefore, there was one parameter for each class in the variable and the parameter was interpreted as the difference from the mean over all classes (i.e., there was no reference class). We also included interaction terms between “species” (categorical) and “starting marine biogeographic realm” (i.e., tagging location), “maturity status” and “mean SST,” and “sex.” These were all included as hierarchical iid effects (i.e., random effects). As these interaction parameters could reasonably explain the full range of the data, we applied penalized complexity priors () of P(SD > 0.5) = 0.01, which favor parsimony (i.e., an SD of 0 for the random effect) and so act to regularize the model. As intrinsic and ecological traits are likely to be more similar between closely related species, we also included a phylogenetic term, based on the one fossil calibration full elasmobranch tree in ()—accessible via www.sharktree.org—and an iid species term. This was implemented as a correlated random effect using the covariance matrix that described the phylogeny. As with other random effects, this model component was hierarchical, with the strength of the effect being determined by a hyperparameter that was learned along with the rest of the model. When using structured random effects such as this, including an unstructured term is important to prevent forcing variance into the structured component and thus confounding inferences. Penalized complexity priors were assigned to both the phylogenetic and species terms with P(SD > 0.5) = 0.01. Last, we included a full hierarchical spatial term with a spatial Gaussian process over the latitude and longitude of the tagging location (decimal degrees). Here, we included penalized complexity priors () on the range [P(range < 10) = 0.01] and the marginal SD [P(SD > 0.1) = 0.01]. From this initial full model, we performed model selection by removing terms and comparing Watanabe–Akaike information criteria (WAIC) scores to compare the predictive fit of each model. We then selected the most parsimonious model (based on the number of effective parameters) within a delta-WAIC of two of the lowest WAIC as the best model.