Ocean warming alters the distributional range, migratory timing, and spatial protections of an apex predator, the tiger shark (Galeocerdo cuvier).

Given climate change threats to ecosystems, it is critical to understand the responses of species to warming. This is especially important in the case of apex predators since they exhibit relatively high extinction risk, and changes to their distribution could impact predator-prey interactions that can initiate trophic cascades. Here we used a combined analysis of animal tracking, remotely sensed environmental data, habitat modeling, and capture data to evaluate the effects of climate variability and change on the distributional range and migratory phenology of an ectothermic apex predator, the tiger shark (Galeocerdo cuvier). Tiger sharks satellite tracked in the western North Atlantic between 2010 and 2019 revealed significant annual variability in the geographic extent and timing of their migrations to northern latitudes from ocean warming. Specifically, tiger shark migrations have extended farther poleward and arrival times to northern latitudes have occurred earlier in the year during periods with anomalously high sea-surface temperatures. A complementary analysis of nearly 40 years of tiger shark captures in the region revealed decadal-scale changes in the distribution and timing of shark captures in parallel with long-term ocean warming. Specifically, areas of highest catch densities have progressively increased poleward and catches have occurred earlier in the year off the North American shelf. During periods of anomalously high sea-surface temperatures, movements of tracked sharks shifted beyond spatial management zones that had been affording them protection from commercial fishing and bycatch. Taken together, these study results have implications for fisheries management, human-wildlife conflict, and ecosystem functioning.

INTRODUCTION

Climate change is now recognized as an environmental emergency, threatening species and ecosystems globally (Lenton et al., 2019; Ripple et al., 2019). In the oceans, many of the habitats facing the highest rates of climate‐driven change are productive areas with high economic importance (e.g., coastal shelves) and biodiversity hotspots (e.g., coral reefs) (Pörtner et al., 2019). The rate and magnitude of future changes are predicted to increase further from current levels (Cheng et al., 2019).

Ocean warming can be driven by both long‐term climate change and short‐term climate variability, such as marine heatwaves, and reports of species’ responses to these climatic phenomena are growing across marine ecosystems (Brown et al., 2016; Pinsky et al., 2013, 2020; Poloczanska et al., 2013, 2016; Scheffers et al., 2016). Observed and predicted responses of species to warming are varied, but commonly include shifts or expansions in their distributional range poleward (Poloczanska et al., 2016). For example, severe marine heatwaves in the northeast Pacific during 2014–2016 triggered poleward distributional expansions in a diverse range of species, including crustaceans, cnidarians, seabirds, and teleosts (Sanford et al., 2019). Climate variability and change is also altering seasonality in the ocean, shifting the annual cycle of surface temperatures toward earlier seasons (Stine et al., 2009). While associated observation of changes in the migratory timing or phenology of marine species is relatively limited (Poloczanska et al., 2016), the few studies investigating this phenomenon have found that seasonal migrations of highly mobile fishes to their northerly range are occurring earlier in the year compared to preceding decades (Dufour et al., 2010; Langan et al., 2021). For example, catches of tunas in the Bay of Biscay during their summer northerly migrations reflect albacore tuna (Thunnus alalunga) and bluefin tuna (T. thynnus) arriving about 8 and 14 days earlier compared to 40 and 25 years ago, respectively (Dufour et al., 2010). Accordingly, key research priorities in climate change ecology are to determine and predict the rate, direction, and timing of associated shifts in the space use and movements of species from climate variability and change.

Quantifying climate‐associated changes to the space use and movements of upper trophic‐level predators is particularly important given that they exhibit relatively high extinction risk (Myers & Ottensmeyer, 2005) and changes to their distributions could render them more vulnerable to exploitation, for example, resulting from shifting ranges outside of protected areas (Selden et al., 2020). Moreover, climate‐driven alteration to the movements of large predators, such as sharks, could change their likelihood of encounters with recreational water users (Chapman & McPhee, 2016) and and/or cause altering ecosystem dynamics through novel trophic cascades (Bastille‐Rousseau et al., 2019; Hammerschlag et al., 2019; Rosenblatt et al., 2017). Although shifts have been predicted (e.g., Birkmanis et al., 2020; Hazen et al., 2013; Niella et al., 2020, 2021), empirical evidence of climate‐driven shifts on the distribution or phenology of marine top predators is rare (c.f. Bangley et al., 2018; Dufour et al., 2010; Tanaka et al., 2021).

The tiger shark (Galeocerdo cuvier) is a globally distributed apex predator in tropical and subtropical seas, with a generalist diet inclusive of teleosts, elasmobranchs, sea turtles, seabirds, and marine mammals (Dicken et al., 2017; Ferreira et al., 2017). Tiger sharks exhibit both wide‐ranging movements in pelagic waters and periods of high residency to coastal habitats (e.g., Acuña‐Marrero et al., 2017; Hammerschlag et al., 2012; Lipscombe et al., 2020; Meyer et al., 2009), with larger individuals generally being more dispersive than smaller conspecifics (e.g., Afonso & Hazin, 2015; Ajemian et al., 2020; Lea et al., 2015). As an ectotherm, temperature is a key abiotic driver of tiger shark habitat use (Ferreira et al., 2015; Lea et al., 2018; Papastamatiou et al., 2013), regulating their coastal abundances and swimming activity levels (Payne et al., 2018). Accordingly, climate‐driven warming is predicted to shift their ranges (Niella et al., 2021). In the western North Atlantic, tiger sharks are distributed coastally from Massachusetts to South Florida and the Bahamas (Bigelow & Schroeder, 1948, 1953; Kohler & Turner, 2018). Here, they exhibit temperature‐driven seasonal migrations (Hammerschlag et al., 2015), moving by way of the Gulf Stream between a southerly cold season (November to April) home range, centered in The Bahamas, to a northerly warm season (May to October) home range that expands to include waters of the Mid‐Atlantic and Southern New England. Taken together, tiger sharks possess traits hypothesized to accelerate climate‐driven range shifts, namely large body size, high mobility, narrow thermal tolerance, negligible risk from predation, a highly generalist diet, and poleward migrations along prevailing currents (Pinsky et al., 2020). Understanding climate‐driven effects on the movement ecology of tiger sharks has been identified as a research priority for this apex predator (Holland et al., 2019). In the western North Atlantic, tiger sharks are relatively well protected from commercial fisheries due to high spatial overlap between their home range and management zones that prohibit commercial longline fishing (Calich et al., 2018; Graham et al., 2016; Queiroz et al., 2019). However, the waters encompassing the northerly portions of their home range are experiencing among the fastest rates of ocean warming globally (Cheng et al., 2019; Pershing et al., 2015), which could consequently shift parts, or all, of their home range outside of protected areas.

In the present study, we used a combined analysis of 9 years of satellite tracking data of tiger shark movements, remotely sensed environmental information, habitat modeling, and nearly 40 years of capture data from conventional tagging to evaluate potential changes in the distributional range, migratory timing, and spatial protections of tiger sharks in response to ocean warming in the western North Atlantic. To accomplish this, we sought to address the following five primary questions: (1) What is the preferred temperature range of the studied tiger shark population? (2) What is the influence of temperature on tiger shark space use relative to other environmental factors known to affect their movements (i.e., chlorophyll a and ocean depth; e.g., Calich et al., 2018)? (3) Does the distributional range of tiger sharks extend farther poleward in response to warming seas? (4) Do seasonal migrations of tiger sharks into their northerly range occur earlier in the year in response to ocean warming? (5) Have climate‐driven shifts in tiger shark space use altered their spatial overlap with management zones that afford them protection from capture in longline fisheries as target and/or bycatch? Addressing these questions will provide insights into the biological responses and associated conservation implications of ocean warming on a marine apex predator.

METHODS

Temperature preferences at swimming depth

To determine temperature and depth preferences of tiger sharks in the study region, pop‐off archival satellite tags were affixed to 10 female tiger sharks in the Northern Bahamas (Nassau and Little Bahama Bank; Table S1) in 2011, 2012, and 2019, permitting measurements of ambient temperatures experienced by sharks at their swimming depths. Both temperature and depth measurements were recorded continuously, but at different intervals depending on tag type, either at 4 min (Sea‐Tag MODS, Desert Star Systems LCC) or 10 s intervals (miniPAT tags, Wildlife Computers, USA or PSATGEOs, Lotek Wireless, Canada; Table S1). All pop‐off satellite tags were either recovered manually, permitting access to the full time‐series, or popped‐off and transmitted their data to an Earth‐orbiting Argos satellite, resulting in a subset of the full time‐series (transmission frequencies: 2.5 min [miniPAT], 10 min [PSATGEO], daily average [Sea‐Tag MOD]). Descriptive statistics of depths and seasonal temperatures experienced by sharks were generated and percent‐frequency histograms of temperature records were plotted using 2°C temperature bins. Additionally, we applied a resource selection analysis to identify the precise range of sea‐surface temperature (SST) where the probability of tiger shark spatial occurrence was highest (see section below on resource selection analysis).

Given known seasonal movement patterns of tiger sharks in the region (Hammerschlag et al., 2015; Kohler & Turner, 2018), data were evaluated with respect to the cold season (November–April) and warm season (May–October).

Tagging and data processing of spatial satellite‐linked data

Between May 2010 and January 2019, 69 tiger sharks were tagged off southeast Florida, southwest Florida, and the northern Bahamas with Smart Position and Temperature Transmitting tags (SPOT tag, Wildlife Computers) to quantify spatial movement patterns. At capture, sharks were sexed and measured for total length (TL). SPOT tags were affixed to the first dorsal fin (Hammerschlag et al., 2012), and tags were coated with antifouling materials (either Propspeed clearcoat, Micron 66, or IPM‐AST4 from Interphase Materials) to minimize biofouling. Prior to deployment, all SPOT tags were tested and confirmed for location accuracy at land‐based facilities.

The geographic location of each tagged shark was determined via Doppler‐shift calculations made by the Argos Data Collection and Location Service (www.argos‐system.org) whenever the shark's tag broke the water's surface and transmitted. Location accuracy was dependent on the number of tag transmissions received by Argos satellites. Argos provides location accuracy using location classes (LC) 3, 2, 1, 0, A, B, and Z (in decreasing accuracy), corresponding with the following error estimates: LC3 < 250 m, 250 m < LC2 < 500 m, and 500 m < LC1 < 1500 m. The error estimates associated with LC A and B are reported to be >1 km and <5 km, respectively (Tougaard et al., 2008). LC Z estimates are inaccurate or unreliable and were removed from the dataset prior to any analysis.

Due to irregular surfacing of sharks (and thus irregular transmission rates) and variation in satellite coverage at any given time, raw SPOT‐derived data are subject to autocorrelation and spatial biases. Therefore, prior to any analyses, positional data were interpolated and regularized to daily estimates using a Bayesian state‐space model that also accounts for Argos satellite telemetry precision (Graham et al., 2016) using the R package foieGras (Jonsen et al., 2019) in the R statistical software (v. 4.0.2.; R Core Team 2021). Individual tracks with data gaps >10 days between positions or tracks with <10 positions were not interpolated (Jonsen et al., 2019).

Given that we aimed to focus analysis on the same population or sub‐population of tiger sharks under study and test one type of climate response strategy (i.e., spatial range shifts), only regularized positions that occurred in the western North Atlantic Ocean were included in subsequent spatial analysis. Two individuals were tracked into the Gulf of Mexico, but these portions of their tracks were not included in the analysis because tiger sharks are unable to shift poleward in the region in response to warming waters. Previous work has shown that some species mitigate the effects of climate‐driven warming in the Gulf of Mexico by shifting deeper to colder water (Pinsky et al., 2013).

Following interpolation, regularization, and filtering of shark positions, 47 (68%) individuals were included in subsequent analyses which examined the spatial preferences and migratory tendencies of tiger sharks and how these patterns related to ocean warming (Table S2).

Resource selection analysis

While the focus of this investigation was to assess how ocean warming may effect tiger shark space use and migrations, previous research in the region has revealed that in addition to sea‐surface temperature (SST), both chlorophyll‐a concentration (ChlA) and ocean depth can influence tiger shark space use (Calich et al., 2018, Queiroz et al., 2016). Therefore, here we evaluated the degree to which these three variables influenced the probability of tiger shark spatial occurrence within each season using a population‐level resource selection function analysis (RSF; Manly et al., 2007). To accomplish this, we first calculated total seasonal home ranges for each individual using standard kernel density estimates based on the “adehabitatHR” R package (we set h=”href,” extent=4, default grid parameter) with the total home range being defined as the 95% isopleth. We then produced 100 random points (i.e., “available” locations) for each shark geolocation (i.e., “used” locations) and distributed random points evenly within each shark's estimated total home range to adequately capture availability (Fieberg et al., 2021). We obtained average seasonal rasters corresponding with each study year from Ocean Color L3 (https://oceancolor.gsfc.nasa.gov/l3/) for SST and ChlA, and an ocean depth raster from NOAA coastwatch (https://coastwatch.pfeg.noaa.gov/erddap/griddap/etopo180.html). We then subsequently extracted SST, ChlA, and ocean depth values at corresponding “available” and “used” locations for each season and year, and we assessed the effect of our variables of interest on resource selection using a generalized linear mixed‐effects model (GLMM) in the R package lme4 (Bates et al., 2014). This model included a random effect of shark nested within year to account for yearly variation between individuals, and a quadratic term for SST (i.e., SST2) to accommodate the possibility of an optimal temperature for tiger shark presence. We evaluated predictor independence using variance inflation factors (VIFs; Zuur et al., 2009). We detected no substantial variance inflation for the three primary predictors (SST, ChlA, and depth); however, we observed VIF>3 for SST2. This is to be expected since SST2 is derived directly from SST, and we retained this predictor to maintain the quadratic structure of the model. From the GLMM, we also calculated optimal seasonal temperature ranges for tiger sharks in the region given estimated model coefficients.

Space use and movement of tracked sharks

To examine the relationships between shark spatial occurrences and surface temperatures, satellite‐derived SST values were obtained from NOAA’s ERDAPP data server (https://coastwatch.pfeg.noaa.gov/), using the R package “rerddapXtracto” (V.0.4.5) (Mendelssohn, 2019) to match SST to each spatially and temporally explicit tiger shark position from SPOT‐tagged sharks. We used the NOAA 1/4° daily Optimum Interpolation Sea‐surface Temperature (OISST) product (https://coastwatch.pfeg.noaa.gov/erddap/griddap/ncdcOisst2Agg_LonPM180.html).

After confirming data‐met assumptions of normality, we applied linear models to test for the effects of SST on the mean latitudes occupied by individual sharks for each season. Each model included the observed individual mean latitudes by year as the response variable and both corresponding SSTs and shark total length as explanatory variables given the latter's previously reported influence on dispersion capacity (e.g., Ajemian et al., 2020; Lea et al., 2015).

To further explore whether and how climate may be driving response of sharks to ocean warming, we extracted SST anomaly values at each spatially and temporally explicit tiger shark position. Sea‐surface temperature anomalies are surface temperature departures from long‐term average temperature, and we derived anomalies from deviations from the SST average between 1971 and 2000 and extracted values for each shark location using the OISST product as described above. We then examined seasonal correlations between the yearly mean latitude occupied by each individual shark and the associated SST anomaly at each location using Spearman correlation tests.

To gain additional spatial insights into how ocean warming may be driving annual variability in tiger shark space use during the tracking period, we generated a composite map of averaged SST anomalies during all days when at least one individual shark was tracked south of 30°N latitude (i.e., within their south range) for each season. This was then compared to analogous composite maps of averaged SST anomalies during all days when at least one individual shark was tracked either north of 35°N or north of 40°N. To enable quantitative comparisons between composite maps, we calculated the area mean of SST anomalies of each composite map.

To evaluate the possible influence of ocean warming on tiger shark migration phenology, we grouped all shark locations by 6° latitudinal bins, encompassing the south (≤30°N), middle (31–35°N), and north (≥36°N) portions of their range. We then identified the minimum day of the calendar year (i.e., Julian date) each tracked shark occurred within each latitudinal bin and the associated SST anomaly at this location. We tested for any correlations between the minimum day of the calendar year an individual tagged shark was tracked within a latitudinal bin and the associated SST anomaly at the location by applying Spearman correlation tests. As all tiger sharks were tagged in the same latitudinal bin (i.e., within their south range), this analysis focused on the habitat selection northward without a bias in tagging location among individuals.

Space use and movement based on captured sharks

To explore the possible effects of longer‐term ocean warming on tiger shark spatial and temporal distributions in the western North Atlantic, we analyzed nearly 40 years of tiger shark capture data obtained from the U.S. National Marine Fisheries Service (NMFS) Cooperative Shark Tagging Program. Initiated in 1962, the Program primarily relies on volunteer fishers who tag and release sharks with provided conventional identification tags (Kohler & Turner, 2018). Most Program participants tagged sharks caught with a rod and reel while fishing recreationally from shore, while other participants included commercial fishers using longline and net gear, biologists, and NMFS observers. When a shark was tagged or recaptured, information on the species, date, location, measured or estimated size, and disposition was provided and placed in a Program‐maintained database (Kohler & Turner, 2018). We used information on the location and date of tag and recaptures of tiger sharks occurring within U.S. western North Atlantic Ocean, spanning 1980–2018 that spatially overlapped with our satellite tracking dataset. Given the reliance on volunteers reporting recaptures to the Program, only recaptures that occurred within waters of the U.S. exclusive economic zone were included in subsequent analyses, representing 90.3% of the data.

To evaluate consistency in thermal associations based on tracking data versus capture data, SST values at each tag‐recapture location were obtained from NOAA’s ERDAPP data server using 1‐day composites (Pathfinder Ver 5.3 [L3C]; https://coastwatch.pfeg.noaa.gov/erddap/griddap/nceiPH53sstd1day.html). Descriptive statistics and percent‐frequency scatter plots of these data were then generated using 2°C temperature bins and fitted with a curve using a combination of first, second, and third‐order models to determine the range, trend, and peak in SST values where tiger sharks were captured for comparison against thermal associations derived from pop‐off archival satellite tags and resource selection function analysis of tracking locations.

To help visualize decadal trends in climate change for comparison with potential decadal trends in shark catch data, we generated composite maps of averaged SST anomalies by decade for both cold and warm seasons. To enable quantitative comparisons between composite maps, we calculated the area mean of SST anomalies averaged by decade for the study area. Here, SST anomalies were based on deviations from historical SST averaged between 1850 and 1979, extracted from the COBE (Hirahara et al., 2014; Ishii et al., 2005). These SST data were a nearly spatially complete, interpolated 1°×1° product, spanning 1850 to present. It combined SSTs from ICOADS (Berry & Kent, 2009) release 2.0, the Japanese Kobe collection, and reports from ships and buoys. Data were gridded using optimal interpolation and bias‐adjusted as needed (i.e., SST from 1850 to 1941) using the “bucket correction” following the approach of Chan and Huybers (2019). Prior to interpolation analyses, data were also subject to quality control using a‐priori thresholds, and nearby observations were combined.

Potential decadal shifts in tiger shark capture locations were investigated by plotting catch densities by decade (1980–2018). To accomplish this, we applied Kernel density analysis to catch data using ArcGIS Pro (version 2.4.0, ESRI 2017), whereby kernel density estimates (KDEs) were generated using a 25 km smoothing parameter and a 1 km2 grid cell (Hammerschlag et al., 2016). We generated 50% density volume contours from KDEs (Worton, 1989), representing high catch density areas. To explore decadal shifts in spatial distribution of high catch density areas in association with long‐term warming trends, plots of high catch density areas were compared with corresponding plots of SSTs averaged by decade. To accomplish this, catch densities were overlaid on top of corresponding SSTs averaged by decade (and season) in ArcGIS Pro and visually compared. Second, we used ArcGIS Pro to manually measure decadal shifts (in km) of the leading (northern) edge of the high‐density catch area (50% KDE) and corresponding decadal shifts in the leading (northern) edge of the underlying preferred temperature isotherm (26–28ºC) of tiger sharks.

To evaluate the possible influence of climate on tiger shark phenology, we grouped all capture locations by the same 6° latitudinal bins as the satellite tracking data, encompassing the south (≤30°N), middle (31–35°N), and north (≥36°N) portions of tiger shark range. We then determined whether any correlations existed between year (1980–2018) and the Julian date of capture averaged across all sharks within each latitudinal bin by applying Spearman correlation tests.

Overlap with protected areas

We examined whether warming‐driven shifts in tiger shark distributions could have conservation implications by quantifying the extent to which variation in satellite tracked shark movements would affect their spatial overlap with protected areas that restrict commercial longline fishing, a source of mortality to tiger sharks in the region (NOAA, 2021). Protected areas had restrictions against bottom and/or pelagic longlines and restrictions were either year‐round or occurred during specific months (i.e., January through July; Table S3, Figure S1). By overlaying shark tracking points with shapefiles of spatial management zones (downloaded from /www.fisheries.noaa.gov) in ArcGIS Pro, we classified all shark tracking points as either “protected” or “unprotected” based on the time of year and the location where a shark was tracked relative to these management zones and also designated whether that protection was from bottom or pelagic longlines or both. SST anomaly values were then determined for each of these tracking positions as described earlier. To understand the potential influence of temperature on whether a shark occurred in a protected or unprotected area, we employed the use of a generalized linear mixed model with a binomial error distribution. The response variable was the probability of an individual shark position occurring within a protective zone and the fixed effects of the model included the SST anomalies with an interaction of season, and the individual shark ID was set as the random effect to account for individual variation in shark behavior. This model was performed for three separate datasets: when sharks were protected or unprotected from only pelagic longlines, only bottom longlines, and both types of longline fisheries. Significance was assessed via the p‐values of the predictor variables.

RESULTS

Temperature at depth preferences

Pop‐off archival satellite tags measured ambient temperature and depth experienced by 10 female tiger sharks (Table S1). Individual PSAT deployments lasted between 14 and 384 days, generating between 771 and 125,166 measurements of both temperature and depth, for a total of 348,938 measurements of each over 1201 tracking days (Figure 1a, Table S1). These data revealed general tiger shark preferences for relatively shallow (<15 m), warm waters above 22°C (78% and 91% of all depth and temperature records, respectively), despite differences in the size and date of individuals tagged (Figure 1a, Table S1). During the cold season, 98% of temperature records were between 22 and 30°C, and 82% of depth records were between 5 and 15 m. Similarly, during the warm season, 83% of temperature records were between 22 and 30°C, and 78% of depth records were between 5 and 15 m.

FIGURE 1

Temperature dependence of tiger sharks in the study area based on (a) temperature measured at swimming depth from sensors on pop‐off archival satellite tags, (b) population‐level resource selection functions using sea‐surface temperature (SSTs) at locations of SPOT‐tracked sharks; (c) SSTs at locations of tiger sharks captured between 1980 and 2018. Blue = cold season (November–April); red = warm season (May–October). In (a), data are from 10 sharks tracked in years 2011, 2012, and 2019 and temperature readings are binned by 2°C. In (b), values are relative log‐odds of selection for SST, where the vertical dashed line represents the SST at which selection is maximized in each season. In (c), values are percent of catches by decade within 2°C temperature bins. The fitted curve in the cold season is a second‐order polynomial: y = −6.2x2 + 41.82x − 36.14 (N = 2932, R 2 = 0.57, p < .001), whereas the fitted curve in the warm season is a third‐order polynomial: y = –2.392x3 + 21.34x2 – 45.39x + 27.76 (N = 5375, R 2 = 0.66, p < .001)

Tiger sharks satellite tracked with SPOT tags generated 5227 locations (Figure 2) across 9 years from 47 sharks (6 males and 41 females), ranging in size from 175 cm to 403 cm TL (mean ± SD = 316 ± 51.6 cm TL).

FIGURE 2

Filtered and interpolated daily position estimates (N = 5227) for 47 tracked tiger sharks. Points are color‐coded by season (blue = cold, red = warm) and graded by study year, starting 1 November 2010 and ending 31 October 2019. (a) Cold season: November–April; (b) warm season: May–October. Provided for spatial reference, two letter abbreviations for US states are FL = Florida, SC = South Carolina, NC = North Carolina, VA = Virginia, PA = Pennsylvania, NY = New York. Tagging locations indicated with yellow crosses

Population‐level resource selection models indicated, for the cold season, all the variables evaluated (SST, ocean depth, ChlA, SST2) significantly influenced the probability of tiger shark occurrence (Table 1). However, the relative effect of SST was the greatest (Table 1), further implicating SST as a key driver of tiger shark occurrence. The negative effect of SST2 (Table 1) indicated that tiger sharks selected for optimum temperatures during the cold season. For the warm season, ocean depth, ChlA, and SST2 also significantly influenced probability of tiger shark occurrence; however, the effect size of ocean depth and ChlA were relatively small and negative (Table 1). There was also a strong negative effect of SST2, indicating a temperature optimum was also selected by tiger sharks in the warm season (Table 1). Modeled thermal response curves based on predicted values suggested the highest probability of tiger shark presence at SSTs between ~25 and 31°C, peaking at 27–28°C, during both cold and warm seasons (Figure 1b), which agrees with shark temperature dependencies based on temperature measurements derived from pop‐off archival satellite tags (Figure 1a).

TABLE 1

Results of the population‐level resource selection functions by season

Variable	Cold season			Warm season
	β ± SE	z‐value	p‐value	β ± SE	z‐value	p‐value
Ocean depth	1.34 ± 0.06	23.83	<.001	−0.36 ± 0.03	−13.82	<.001
Chlorophyll‐a	−0.11 ± 0.04	−2.63	.009	−0.09 ± 0.03	−3	.003
Sea‐surface temperature	4.66 ± 0.51	8.97	<.001	0.13 ± 0.25	0.54	.59
Sea‐surface temperature2	−3.49 ± 0.46	−7.48	<.001	−0.9 ± 0.26	−3.46	<.001

β ± SE are slopes and associated standard errors. Relative differences in effect size among variables in the model can be assessed by comparing relative differences in the β values. Cold season = November through April; warm season = May through October.

Visual inspection of the raw tracking data revealed seasonal and annual variation in tiger shark movement patterns (Figure 2); however, mean latitudes occupied by sharks were inversely related to mean SSTs, demonstrating general preferences for warmer waters that occur in the south portions of their range (Figure 3a). Mean latitudes occupied by tiger sharks only extended north of ~35°N (the north portion of their range) when sharks experienced mean SSTs at or above ~22°C (Figure 3a). We found an effect of SST (p < .0001, coefficient estimate = −1.21), but not shark length (p = .91), on annual mean latitudes occupied by sharks during the cold season (R 2 = 0.39). We found an effect of SST (p < .0001, coefficient estimate = −1.55; Figure 3a) and shark total length (p = .03, coefficient estimate = 0.018) on annual mean latitudes occupied during the warm season (R 2 = 0.67). In the case of the former, for every 1°C increase in SST, the linear model predicted a 1.6° decrease in mean latitude occupied, whereas in the case of the latter, for every 1 m increase in shark length, the linear model predicted a 0.02° increase in mean latitude.

FIGURE 3

(a) Linear relationship between the yearly mean latitude occupied by tiger sharks and the associated sea‐surface temperature (SST). (b) Correlation between the yearly mean latitude occupied by tiger sharks and the associated SST anomaly. Cold season (November–April) = blue; warm season (May–October) = red. Data are means for each individual shark by year and season. Shaded area represents the 95% confidence intervals around a significant relationship

Mean latitudes occupied by sharks were positively related to SST anomalies. Mean latitudes occupied by sharks extended north of ~35°N (their north portion of their range) when SST anomalies exceeded mean values of ~1.5°C (Figure 3b). Annual variation in mean latitudes occupied by tiger sharks was not significantly correlated with SST anomalies in the cold season (r = 0.18, p = .24), but was positively correlated in the warm season (r = 0.48, p < .002; Figure 3b). A linear model plotted on these data predicted a 3.71°N increase in mean latitude occupied by tiger sharks for every 1°C increase in SST anomaly.

Consistent with these results (Figure 3b), the composite maps of average SST anomalies (Figure S2) revealed that tiger shark distributions extended successively farther north during periods of increasingly anomalous warm water. The area means of the SST anomalies averaged across days when at least one individual tiger shark was tracked south of 30°N (N = 2967 positions), north of 35°N (N = 998 positions), or north of 40°N (N = 290 positions) were 0.34, 0.73, and 0.92°C, respectively (Figure S2).

We found a negative correlation (r = −0.68, p = .01; Figure 4) between SST anomalies and the minimum Julian date individual sharks were tracked in northeast shelf waters (i.e., only in their north range: latitudinal bin ≥36°N). A linear model plotted to these data predicted that for every 1°C increase in SST anomaly, the minimum Julian date a tracked tiger sharks entered northeast shelf waters decreased by ~14 days. No other latitudinal bins exhibited a correlation (south range: r = −0.03, p = .84; middle range: r = −0.14, p = .45; Figure 4).

FIGURE 4

Tiger shark migration phenology in response to ocean warming based on SPOT satellite tracking locations between 2010 and 2019; plots are correlations between the minimum day of the calendar year (i.e., Julian date) an individual tagged shark was tracked within a latitudinal bin and the associated sea‐surface temperature (SST) anomaly at the location. Location data are grouped by 6° latitudinal bins, generally representing the south, middle, and north range of tiger sharks. Shaded area represents the 95% confidence intervals around a significant correlation

Space use and movement based on shark captures

Between 1980 and 2018, a total of 8764 tiger sharks were captured and tagged and/or recaptured as part of the Cooperative Shark Tagging Program (Figure S3). Analyses of these data (Figure 1c) also revealed a temperature dependency pattern in agreement with results from the tracking data (Figure 1a); specifically, tiger shark capture rates were also highest in coastal waters with SSTs between 22 and 28ºC. Such a temperature dependency was consistent across decades (Figure 1c).

The composite maps of decadal SST anomalies since the 1980s revealed increasing deviations in SSTs from historical (1850–1979) climatology (Figure S4). Across the cold seasons, the area means of SST anomalies averaged by decade progressively increased: 0.04ºC (1980s), 0.06ºC (1990s), 0.10ºC (2000s), and 0.29ºC (2010s), respectively (Figure S4). Across the warm seasons, the area means of SST anomalies averaged by decade also progressively increased: 0.01ºC (1980s), 0.16ºC (1990s), 0.32ºC (2000s), and 0.75ºC (2010s), respectively (Figure S4).

Kernel density analysis on tiger shark capture locations revealed a decadal poleward shift in the spatial distribution of high‐density captures consistent with decadal poleward shifts in preferred SSTs (Figure 5). Specifically, across decades in cold seasons, the northern edge of their optimal SST isotherm (26–28°C) shifted ~300 km poleward. In parallel, the leading edge of the high‐density catch area (50% KDE) also shifted poleward by ~300 km (Figure 5). Across the decades in warm seasons, the northern edge of their optimal temperature SST isotherm shifted ~400 km poleward. In parallel, the leading edge of the high‐density catch area shifted poleward by ~440 km (Figure 5).

FIGURE 5

Decadal patterns in high catch density areas of tiger sharks and associated decadal averages in sea‐surface temperatures (SSTs). High catch densities are based on 50% density volume contours from kernel density analysis on 8764 tiger shark captures between 1980 and 2018. SSTs are averaged by decade for both cold seasons (November–April; top row) and warm seasons (May–October; bottom row). Provided for spatial reference, two letter abbreviations for US states are FL = Florida, SC = South Carolina, NC = North Carolina, VA = Virginia, PA = Pennsylvania, NY = New York State. Underlying SSTs are on a gray scale; the vertical arrow in the legend indicates direction of temperature preference from lowest to highest (white to black)

We found a negative correlation (r = −0.40, p = .01; Figure 6) between year and the mean Julian date sharks was captured in northeast shelf waters (i.e., only in their north range: latitudinal bin ≥36°N); no other latitudinal bins exhibited a correlation (south range: r = −0.29, p = .07; middle range: r = −0.04, p = .82; Figure 6).

FIGURE 6

Tiger shark migration phenology derived from capture locations between 1980 and 2018. Plots are correlations between year and the mean calendar day of the year (i.e., Julian date) of shark capture within a latitudinal bin. Location data are grouped by 6° latitudinal bins, generally representing the southern, middle, and northern range of tiger sharks. Shaded area represents the 95% confidence intervals around a significant correlation

Overall, the probability of an individual shark position occurring within a protected area decreased with increasing SST anomalies (Table 2, Figure 7). With respect to management areas prohibiting both benthic and pelagic or only pelagic longline fishing, there was a significant interaction between SST anomaly and season on the probability of an individual shark position occurring within a protective zone (Table 2). Both models revealed a negative relationship between SST anomalies and the probability of protection, while the probability of a shark position occurring in a protective management zone was higher during the cold season compared to the warm season (Table 2, Figure 7a, b). With respect to management areas prohibiting only benthic longline fishing, there was a significant negative influence of SST anomalies on the probability of protection as well as significantly higher protection during the cold season (Figure 7c). However, the interaction between SST anomalies and season was not significant (Table 2).

TABLE 2

Results from generalized linear mixed models describing the relationship between sea‐surface temperature (SST) anomalies on the presence satellite tracked tiger sharks with an interaction of season (i.e., cold and warm) in management zones where pelagic longlines (PLL), bottom longlines (BLL), or both PLL and BLL were prohibited

Predictors	PLL and BLL			PLL			BLL
	Odds ratios	CI	p	Odds ratios	CI	p	Odds ratios	CI	p
(Intercept)	2.15	1.08–4.28	<.029	0.88	0.50–1.55	<.660	0.59	0.31–1.12	<.104
SST Anomaly	0.84	0.78–0.91	<.001	0.71	0.65–0.77	<.001	0.86	0.80–0.93	<.001
Season	21.89	17.15–27.96	<.001	19.16	15.35–23.93	<.001	8.30	6.92–9.95	<.001
SST Anomaly × Season	0.85	0.74–0.98	.023	0.69	0.60–0.79	<.001	0.92	0.82–1.03	.141

Cold season = November through April; warm season = May through October. ICC = interclass correlation coefficient, σ2 = variance of the random effect, N = sample size. For the seasonal comparison, the cold season is the reference.

FIGURE 7

Relationship between sea‐surface temperature anomalies and the probability of a satellite tagged tiger shark position occurring within a management zone prohibiting (a) both pelagic and bottom longline fishing, (b) only pelagic longlines, and (c) only bottom longlines, in the cold season (November–April; blue) and warm season (May–October; red). Shaded regions indicate 95% confidence intervals around a significant relationship as predicted by generalized linear mixed models. Shark tracking data from 2010 to 2019

DISCUSSION

The present study used multiple lines of parallel evidence to demonstrate alterations in the distributional range and migratory timing of an ectothermic apex predator, the tiger shark, from climate variability and change within the western North Atlantic Ocean. Moreover, changes in the space use of tiger sharks from ocean warming have subsequently decreased their spatial protections from commercial fisheries, increasing their vulnerability to potential exploitation. We discuss these lines of evidence in the sections that follow.

Temperature preference

Our combined analyses of satellite tracking, habitat modeling, and capture data suggest a strong preference of tiger sharks for waters between 26 and 28°C, with ~22°C representing a lower thermal optimal limit for tiger sharks tagged in the region (Figure 1). Interestingly, previous studies from Australia suggest that ~22°C is the thermal optimal for this species in the Pacific (Payne et al., 2018). Consistent with our data, other studies have reported that Atlantic tiger sharks have higher thermal performance optima than those from the Pacific (Lear et al., 2019). While beyond the scope of this study, it is plausible that tiger sharks exhibit geographic adaptation to temperature variation, as many species can locally adapt to climate (Dixon et al., 2015; Yampolsky et al., 2014).

Spanning the warmest decade on record for ocean temperatures (Osborn & Jones, 2019), satellite tracking of tiger sharks over 9 years revealed annual variation in space use and migratory patterns linked to ocean warming. Tiger sharks migrated earlier and farther in northern latitudes during years and seasons in which northeast shelf waters experienced high SST anomalies such that these waters warmed within preferred temperature ranges of tiger sharks. These data suggest that on average, for every 1°C increase in SST anomalies, tiger sharks have extended their movements farther north by nearly 4 degrees of latitude and arrived in northeast shelf waters about 14 days earlier.

In addition to the effects of temperature on tiger shark space use, we also found an effect of shark length such that larger individuals ranged farther north than smaller conspecifics during the warm season. This result is consistent with previous tracking studies demonstrating higher dispersive capabilities in larger tiger sharks (Afonso & Hazin, 2015; Ajemian et al., 2020; Lea et al., 2015), possibly related to combination of proposed factors such as higher movement capacity, relatively lower cost of transport, and lower metabolic rate per unit mass associated with larger body size (Lea et al., 2018; Sorte et al., 2013). However, shark length had no effect on mean latitudes occupied by sharks during the cold season and the effect size of SST on latitudes occupied by sharks during the warm season was nearly two orders of magnitude greater than that of shark length, implicating ocean warming as a greater driver of shark space use. It is plausible that other biological (e.g., reproductive state, health, energy levels; Hammerschlag et al., 2018; Rangel et al., 2021), demographic (e.g., age and sex; Niella et al., 2021), ecological (e.g., increased intraspecific competition for resources and mates; Papastamatiou et al., 2018; Ratcliffe et al., 2018), and environmental factors (e.g., oceanic fronts, dissolved oxygen; Lea et al., 2018; Queiroz et al., 2016) not evaluated here contributed to some of the variation observed in movements of tracked tiger sharks. That said, our interpretation of ocean warming as a key driver of variability in space use is strengthened by congruence among results of tracking data, habitat use measured via resource selection function, as well as multi‐decadal trends in tiger shark capture data.

It is also possible that tiger shark distributional patterns found here could be related in part to climate‐driven changes in their prey distribution, rather than ocean warming directly. Marine consumers, including sharks, have been shown to associate with areas of higher primary productivity (measured through ChlA), which directly supports local food webs and thus serves as proxy for prey availability and foraging opportunities (Block et al., 2011; Queiroz et al., 2016). While the conducted resource selection models revealed an effect of ChlA on utilization distributions, the effect size was relatively small and negative. Given water visibility is generally higher under lower ChlA conditions, visual predators such as tiger sharks may avoid ChlA‐rich waters to maximize visual acuity and consequently prey capture success (Ortega et al., 2020). Even if prey may not be the primary driver of the annual variability in space use of the tiger sharks tracked here, warming on the northeast shelf has already led to increases in species diversity and productivity (Friedland et al., 2020), which may offer increased foraging opportunities for tiger sharks that range there.

Capture locations are influenced by both the home range of sharks and the behavior of fishers, the latter of which includes the operational variables by which fisheries interactions occur. While catch data alone cannot be used to define shark home range, catch densities can provide insights into the distribution of tiger sharks within their range and can be used to subsequently detect climate effects on space use and movement (e.g., Dufour et al., 2010). Here, analysis of nearly 40 years of capture data suggests northward expansion in tiger shark distributions in relation to SST, matching observations of ocean warming (Cheng et al., 2019). Not only have high catch density areas shifted poleward in parallel with expansions in preferred water temperatures across decades, catches of tiger sharks off the northeast shelf have also progressively occurred earlier in the year, a pattern consistent with known effects of warming on the migration phenology of ectothermic teleosts in the region (Langan et al., 2021).

Space use shifts can place species at risk even when successful management approaches have been implemented to conserve them (Davies et al., 2019; Santora et al., 2020). In the western North Atlantic, no specific conservation or management measures exist for tiger sharks, but vulnerability of this species to fisheries has been relatively low in recent decades due to overlap between their highly suitable habitat areas and conservation management zones that prohibit commercial longline fishing (Calich et al., 2018; Graham et al., 2016; Queiroz et al., 2019), which has likely contributed to their population recovery and growth (Morgan et al., 2020). Here, we found that spatial overlap between tracked tiger sharks and these protective zones varied from as much as 100% to as little as only 5% in response to ocean warming over the past decade. Specifically, tiger sharks shifted poleward outside of protected areas during times of high SST anomalies. While individual variation was controlled for in our analysis as a random effect, individuals tracks varied in spatial distribution throughout the seascape, as would be expected. This contributed to relatively high intraspecific variability in the probability of shark positions overlapping with protected areas. Despite this individual variability in behavior and space use, it is clear that ocean warming decreased spatial protections from commercial longline fishing due to poleward expansion of shark movement.

While tiger shark populations in the study region are currently recovering from historical declines (Peterson et al., 2017), increased fisheries interactions could lead to declines as has been found for this species in other areas (Holmes et al., 2012; Roff et al., 2018). Indeed, previous modeling has predicted that tiger shark biomass in the western North Atlantic will be negatively affected under a scenario of decreased overlap between their highly suitable habitats and no‐fishing zones (Morgan et al., 2020).

CONCLUSION

Here we show congruence among multiple data sources, including satellite tracking data, habitat models incorporating remotely sensed environmental data, and multi‐decadal capture data. Our results provide evidence for alteration in the distributional range and migration phenology of tiger sharks in the western North Atlantic, a region strongly affected by climate variability and change. Given their role as apex predators in tropical and subtropical seas, these climate‐driven changes in tiger shark space use and migratory patterns could lead to shifts in ecological interactions through alterations in predator–prey dynamics, and may also redistribute sharks into areas where human–wildlife conflict is likely to occur. Therefore, understanding, predicting, monitoring, and ultimately mitigating such impacts is both a fundamental research and socio‐economic priority. Of conservation concern, our data indicate that ocean warming has already shifted tiger shark movements outside of management zones that protect them from capture in target and bycatch fisheries. Ensuring management approaches are responsive to changing environmental conditions, species distribution, and fishing effort is required to ensure sustainable fisheries in a changing climate (Hazen et al., 2018; Holsman et al., 2019; Pinsky & Mantua, 2014).

CONFLICT OF INTEREST

The authors declare no competing interests.

AUTHOR CONTRIBUTIONS

NH conceived the study, conducted the fieldwork, analyzed the data and wrote the manuscript. LHM, MJR, MRB, and GMS also contributed to data analysis; LHM contributed to figure preparations; AJG contributed to fieldwork, data analysis, and figure preparation; BK contributed to data analysis; LJN and CMT provided capture data from the NMFS Cooperative Shark Tagging Program. All authors contributed significant input and revisions to the manuscript.

Supplementary Material

Click here for additional data file.