Human use of oceans has dramatically increased in the 21^st century. Sea turtles are vulnerable to anthropogenic stressors in the marine environment because of lengthy migrations between foraging and breeding sites, often along coastal migration corridors. Little is known about how movement and threat interact specifically for male sea turtles. To better understand male sea turtle movement, and the threats they encounter, we satellite-tagged 40 adult male sea turtles of four different species. We calculated movement patterns using state-space modeling (SSM), and quantified threats in seven unique categories; shipping, fishing, light pollution, oil rigs, proximity to coast, marine protected area (MPA) status, and location within or outside of the U.S. Exclusive Economic Zone (EEZ). We found significantly higher threat severity in northern and southern latitudes for green turtles (Chelonia mydas) and Kemp’s ridleys (Lepidochelys kempii) in our study area. Those threats were pervasive, with only 35.9% of SSM points encountering no high threat exposure, of which 47% belong to just two individuals. Kemp’s ridleys were most exposed to high threats among tested species. Lastly, turtles within MPA boundaries face significantly lower threat exposure, indicating MPAs could be a useful conservation tool.

2.1 Study Area/Species Collection

We captured turtles as in Hart et al. [100] from 2009 – 2019. Forty adult male sea turtles of four different species were captured from four locations using a boat (jumping from a boat, snorkeling) or net capture via trawler. Sample sizes are as follows: Kemp’s ridley = 6, hawksbill = 1, loggerhead = 8, green = 25. Capture location sample sizes are as follows: Dry Tortugas National Park = 24, Florida Keys National Marine Sanctuary = 6, Northern Gulf of Mexico = 9, Buck Island National Reef Monument = 1. We followed standard morphometric data collection, and attached platform transmitter terminals (PTT) to each turtle carapace using slow-curing epoxy (two-part Superbond epoxy; see Hart et al. [100]). Turtles were tracked using Wildlife Computers (Redland, Washington, U.S.A.) SPOT or SPLASH10 transmitters. Tracking data ranged from 8 June 2009 to 7 August 2020 [100, 101].

2.2 Collection and Calculation of Threats/State-Space Modelling

We performed a switching state space model (SSM) on the raw location data in order to estimate each turtle’s true locations at regular time intervals due to significant positional uncertainty in the raw satellite data [101]. Briefly, we used a Bayesian hierarchical movement model implemented in the R package ‘bsam’, using the ‘hDCRWS’ model specification and a time step of 1 day [103, 104, 105].  We set the Markov Chain Monte Carlo (MCMC) parameters following Roberts et al. [106], which used adaptive sampling for 7000 draws, taking 10,000 samples from the posterior distribution, and then thinning by five to reduce MCMC autocorrelation, resulting in 2000 posterior samples from which to make inference. This process ultimately resulted in an improved dataset by eliminating location errors and provided one location point for each turtle per day.

We collected a total of 8875 SSM points for threat analysis [101]. Through review of scientific literature and professional consultation, we collected data for seven primary threats to male sea turtles (Fishing, Shipping, Drilling Platforms, Light Pollution, MPA boundaries, located within or outside the U.S. Exclusive Economic Zone (EEZ), and coastal threat (within 10 km of a coastline) [8, 35, 36, 103, 107, 108, 109]. Raw location data have spatial accuracy ranges that vary between 500 m to 1.5 km. Raw tracking data were therefore fit to a hierarchical, behavior-switching state-space model (SSM), which was then used to increase the accuracy of tracking data and to determine home ranges of each individual [103]. This allowed for estimation of the behavioral modes of individual turtles (unique behavioral patterns), regularization of the locations in time, and accounting for location error in the raw satellite data. In order to accurately depict the threats within the area of each SSM point (1.5 km), we created a 2-km radius buffer around each SSM turtle point using R [105], within which threats were assessed. The threat data were collected and prepared as described in the following sections.

2.3 Fishing Data

Threats from fishing can come from a variety of sources (artisanal, longlining, commercial, nets and trawlers, etc.). Since 2016, all commercial fishing vessels within U.S. waters over 65 feet in length are required to have an AIS (Automatic Identification System) transponder tag attached, which tracks the ships every hour via satellite global positioning system (GPS) and ground-based receivers placed by the U.S. Coast Guard [111]. At present, only 2% of the world’s fishing vessels have AIS tags on board, but these ships account for more than 50% of total fishing effort [111]. We used a fishing density raster layer of ground-based, AIS tracked ships as a representative subsample of fishing fleet intensity from marinecadastre.gov, a joint collaborative data repository for marine-related research by the National Oceanic Atmospheric Administration (NOAA), and the Bureau of Ocean Energy Management (BOEM) [112]. This layer includes tracks of fishing vessels that leave U.S. waters in the Gulf of Mexico and therefore provides a representative subsample for turtles that move beyond the U.S. EEZ boundary. 

Fishing intensity was measured in grid cells 1 hectare in size, with each cell representing the total number of fishing craft that passed through that cell with an AIS transponder onboard within a given year. We added and averaged the total number of fishing vessels per cell from 2015–2017 to get the mean fishing intensity per cell. We then created a single raster layer for analysis using ArcGIS Pro ver. 2.5 [113]. Although only three years of fishing data were available, we assume fishing density was similar enough in previous years that the average cell values of the three years represent past fishing seasons. We averaged the fishing intensity score of all raster cells within each SSM turtle buffer and assigned that value as the fishing threat score for that point in ArcGIS Pro ver. 2.5 [113].

2.4 Shipping Data

Shipping data were also obtained using AIS tagged ships, and were downloaded in their raw format courtesy of the U.S. Coast Guard with certain identifiers scrubbed for privacy. We were able to obtain data for the entire study period within from the Marine Cadastre data repository [112]. These data cover the entire study area and are a suitable, representative sub-sample of shipping data, clearly showing all shipping lanes within our study area.

     AIS-tagged ships in the United States account for 50–60% of shipping activity [110]. We clipped all AIS data to each 2-km turtle buffer by date and then merged the data into a single vector shapefile to create a layer of shipping points that coincide with the presence of each 2-km SSM turtle buffer. Because of the large size of the AIS data files, we ran an RStudio instance on Google’s Cloud Computing Engine [105]. The total number of shipping points within each 2-km buffer was then added and assigned as the threat score for that SSM point.

2.5 Drilling Platform Data

We downloaded drilling platforms point data, also referred to as oil rigs, oil platforms or drilling rigs to represent oil derived threats from the Marine Cadastre data repository [112]. We calculated the number of platforms within each turtle buffer by clipping the oil rig layer to the turtle layer and merging the data into a single Vector shapefile to create a layer of oil rig points within each SSM turtle buffer using ArcGIS Pro ver. 2.5 [113]. Drilling platforms were corrected by date to ensure they were in use during the date associated with the date of the 2-km turtle buffer.

     Marine Cadastre, although very useful in acquiring data for the United States, is missing drilling platform data for other parts of our study area, specifically Mexico and Cuba. In order to understand if turtles that left the U.S. EEZ encountered oil threats, we used a world, oil exploration shapefile called PETRODATA, that covers oil drilling hotspots around the world [114]. Upon comparison with our existing dataset, 99.7% of oil rig points fall within the PETRODATA polygon for the United States, Gulf of Mexico, oil exploration polygon. Therefore, we felt it was comparable because no public oil rig data are available for Mexico at the time of writing this manuscript. No international turtle points fell within the confines of oil polygons so further calculations were not necessary.

2.6 Light Pollution Data

In 2011, the SUOMI VIIRS (Visible Infrared Imaging Radiometer Suite) satellite was launched to track multiple spatial data, such as snow and sea ice cover, active wildfires, sea and ice surface temperatures, and day/night light reflectance and radiance at high resolution [115]. We created our light pollution threat layer by combining all available light radiance raster data from NOAA’s Earth Observation Group public download domain and averaging the total radiance for each pixel [115]. In total, 54 raster files, ranging from January 13^th – March 8^th, 2021, were combined by taking the average light radiance of each pixel, and then recording the average value of all pixels within each 2-km turtle buffer user R [105]. We assumed that the light radiance during the study period did not vary annually and that the light data we collected are representative of all years for which we have tracking data.

2.7 MPA, EEZ, and Proximity to Coast Data

We downloaded both the MPA Layer and EEZ Layer as vector layers from the Marine Cadastre data repository [112]. We created the coastal threat vector layer in ArcGIS Pro version 2.5 by making a 10-km buffer around all available land within the study region [113]. The SSM points that were > 10 km from coastline, or within the U.S. EEZ or an MPA boundary were assigned a value of “0” to indicate threat absence, whereas points that were < 10 km from the coast or outside of the U.S. EEZ or an MPA boundary were assigned a value of “1” to indicate threat presence.

     Because Marine Cadastre focuses on primarily U.S. waters, we needed data on international MPAs, specifically for Cuba and Mexico. Those data were downloaded from the IUCN’s World Database on Protected Areas website [116]. We followed the same format of assigning values of “0” for turtles that were within the confines of those MPAs, and values of “1” for turtles outside the confines of those MPAs.

2.8 Statistical Analyses

To directly add and compare the effects of individual threats, we standardized all threat categories to a mean of zero and standard deviation of one, which allowed us to take the sum of all threats directly and create a cumulative threat score for each SSM point. Through preliminary data exploration, we discovered the data were substantially non normal, and greatly spatially autocorrelated. As a result, we removed SSM points that were within 4 km of one another. Data removal reduced the number of SSM points from 8875 to 474. Despite removal, spatial autocorrelation still existed, but the degree to which it existed was reduced. Moran’s I of cumulative impact scores changed from 0.458 (p < 0.001) at distance class I to Moran’s I of 0.313 (p < 0.001). Complete removal of spatial autocorrelation from our dataset would have thinned the data to too few points to be able to run an analysis on; therefore, we decided to strike a balance between reducing the data, yet also minimizing the degree to which spatial autocorrelation existed in our dataset.

To test the prediction that individual or combined threats varied with species along a latitudinal gradient, we tested our thinned data using PERMANOVAs, a permutation-based test for significance as the data were still very non normal. PERMANOVAs were run using the R package ‘vegan’ [117]. We tested this using the interactive effects of species and latitude responding to threat. We included the threat x species interaction because we expected that species may vary in their response to spatially varying threats. Because of the potential for areas with variables of high threat to be clustered, data may not show a direct, linear relationship with latitude. To better understand any latitudinal relationships present in our data, we ran breakpoint regressions between latitude and threat (individual threats vs. latitude and cumulative threat scores vs. latitude) to see if the breakpoint model was a better fit to the data. Due to the large spatial gap present and low sample size, we removed turtle 14, the lone hawksbill of this study from this portion of the analysis. All data were given an alpha of 0.95 for detecting statistical significance. All data were analyzed using R [105].    

In order to test the prediction that threats vary by species, we first calculated median values (due to the presence of outliers) of threat scores for each of the four numerical threats (Shipping = 0, Light = 1.3, Fishing = 1.2, Oil = 0) and used presence scores for the remaining three (Coast = 1, MPA = 1, EEZ = 1). Values above the median value for numerical threats or that had a score of 1 for categorical threats were categorized as ”high threat,” whereas values below the median or a score of zero, respectively, were considered ”low threat.” We then calculated the percentage of days during the study period an individual turtle encountered high threats by dividing the number of days a high threat was encountered by the sum of their SSM points. We calculated average values for each species from these percentages to understand how often sea turtle species were exposed to each threat during the study period. Preliminary data analysis discovered our data were very non normal. Therefore, we ran PERMANOVAs on each threat category percent by species. Turtle 14, the single hawksbill captured for our study, was removed from this part of the analysis.

To test the prediction that turtles within MPA boundaries experienced lower threat, we recorded the mean time each individual turtle spent outside of an MPA using their SSM points. If individuals spent more than the mean value (23.6%) outside of an MPA, they were counted as a “non-MPA” turtle. Individuals that spent more time than the mean value within an MPA were counted as “MPA” turtles. Three threat variables of continuous data (Light, Shipping, Fishing) that were categorized as high threat were compared between non-MPA and MPA turtles for statistical significance with a Welch’s T-Test in R [105]. T-Test analyses were modified for unequal variance and if the equal variance assumption was violated.

To test the prediction that concentration of turtle SSM points is lowest in areas of high threat, we created a 10x10 km grid cell fishnet over the study region and then using the ‘Spatial Join’ tool to merge all points within each grid of the fishnet in ArcGIS Pro version 2.5 [111]. The number of turtle points within each grid cell was treated as a form of density for that specific area. Threat values of each grid cell were averaged if more than one turtle point was present. We then ran linear regression analyses to test for relationships between density and each individual threat, and several combinations of threat layers: All Combined Threats, Oil Threat (Oil, Shipping, and Coastal layers), Boat Threat (Fishing, Shipping, and Coastal, EEZ layers), and Fishing Threat (Fishing, Coastal, EEZ, and MPA layers). We additionally ran breakpoint regressions on our data to determine if density responded to threat nonlinearly or in linear segments using the R package ‘segmented’ [118].

To better understand threat interactions through time, we added the total number of high threat categories encountered on a given day for each SSM point and then plotted the threat on a color gradient (from 0–6) by month of the year and individual turtle, which clearly displays the daily number of high threats each turtle was exposed to during the tracking period. Turtles were sorted by capture location, MPA status, and species. Plotting data in such a manner can provide a visual interpretation of threats by individual, as well as show the timing of threats by month. All data were analyzed using R version 4.1.0 [105].

Please see the attached metadata release for navigating our dataset. Attached is all the transformed data, locations where raw data was accessed, and code for figures and statistical analysis. Coordinates have been removed to protect sensitive species locations.Therefore some of these data transformations in the code cannot be run. Requests for coordinate data can be made to Dr. Kristen Hart: kristen_hart@usgs.gov