Environmental and Biological Factors Influencing Dispersal of Neonate Leatherback Turtles (Dermochelys coriacea) From an Endangered Costa Rican Nesting Population

Quantifying early life movements is essential to understanding migratory pathways and habitat use that can impact individuals’ success later in life. To gauge how neonatal movements set the stage for later habitat use, we tracked neonate leatherback turtles (n = 94) with acoustic tags from Pacuare, Costa Rica, in 2016 and 2018. We analyzed movements using a first passage time analysis and random walk models, the results of which indicated neonates followed a fixed compass direction as they traveled away from shore and that strong currents in these areas resulted in advection. We combined the tracking data with concurrent environmental variables in a generalized additive mixed model framework. Our results showed the south-east current flow in this area has spatial and temporal structure consistent with large-scale geostrophic currents and not tidal current or local wind speed influences. After accounting for advection by currents, true neonate swimming speed was significantly related to current speed, first passage time, and the year. Neonates had three main response strategies to currents above 0.5 m s–1, with most increasing their swimming speed and the rest maintaining either a constant or decreased swimming speed. Neonates were significantly larger in 2018 than in 2016 but their average swimming speed was not significantly related to body size, indicating that environmental factors were more important contributors to their dispersal. We conclude that abiotic factors, including the strength and direction of the currents, significantly affect the swimming and dispersal strategy of neonate leatherback turtles and these results can help to inform strategies for releases of neonate turtles from hatcheries, future tracking studies, and conservation efforts.


INTRODUCTION
Advancements in biotelemetry have helped bridge the gap between animals' fine and broad-scale movement patterns, opening up new opportunities for the movement ecology of different life stages to be quantified and understood (Nathan et al., 2008). Biotelemetric monitoring is especially important for understanding individuals' movements in settings where direct observations are complicated or obscured, such as the marine environment. However, the need to attach biotelemetric instruments either directly to, or in, an organism has created limitations for studying the movements of young life stages, which may experience larger tag-to-body size ratios and higher natural mortality rates Shillinger et al., 2012a). Although satellite tags have become progressively smaller in size and capable of attachment to the larger juvenile life stages for animals such as birds, sharks, sea turtles, and marine mammals (Kjellén et al., 2001;Vincent et al., 2002;Hsu et al., 2007;Weng et al., 2007;Mansfield et al., 2014), their application on the much smaller, dispersing neonate phase is limited. Tagging of this vulnerable life stage has most often been applied to the fledgling stage in birds Hake et al., 2003;Kooyman and Ponganis, 2007;Cadahía et al., 2010;Péron and Grémillet, 2013;Vega et al., 2016). This leaves a knowledge gap for many other migratory species, wherein the links between neonate fine-scale distribution and adult long-term movement paths remain unknown (Hays et al., 2010). Characterizing this early life stage is important, as the movement and behavioral responses of neonates to different environmental conditions could potentially have carryover effects on the distribution, energetics, or survival of later life stages (O'Connor et al., 2014).
The habitat uses of adult marine animals has been characterized using biotelemetry in combination with advanced statistical models (Block et al., 2011;Hazen et al., 2012a) but there is a need for the application of these methods to other life stages. Statistical tools such as first passage time analyses can allow for the detection of changes in speed in response to environmental conditions at different spatial scales (Fauchald and Tveraa, 2003) and mechanistic movement models such as random walk models (Turchin, 1998), can provide further insight into behavioral processes and states. The application of such approaches to the movement trajectories of marine species has allowed for migratory and foraging behavior (Bailey and Thompson, 2006;Hurme et al., 2019) and navigational abilities (Girard et al., 2006) to be inferred for adult life stages.
Sea turtles are highly migratory species and research on their dispersal has been biased toward adult breeding females, since they can be satellite tagged on the beach and tracked globally for relatively long periods of time (Hays et al., 2006;Shillinger et al., 2008;Shillinger et al., 2011;Fossette et al., 2014). The behavior and habitat associations of neonate turtles during their early dispersal has been characterized for some species using tracking technology (Okuyama et al., 2009;Gearheart et al., 2011;Thums et al., 2013;Scott et al., 2014;Thums et al., 2016), field observations (Frick, 1976;Lohmann and Lohmann, 1992), and laboratory studies Wyneken et al., 1990;Goff et al., 1998;Smith and Salmon, 2009). Sea turtles have been shown to disperse from the nesting beach in a non-random fashion (Pilcher et al., 2000), using a mixture of geomagnetic information (Light et al., 1993;Lohmann and Lohmann, 1993;Lohmann and Lohmann, 1994), and environmental cues, such as the direction of surface waves, to set a fixed compass direction into deeper waters (Frick, 1976;Lohmann et al., 1990;Wyneken et al., 1990;Lohmann and Lohmann, 1996;Pilcher et al., 2000). Neonate sea turtles undergo a "frenzy" period of continuous and highly energetic swimming lasting around 24 h during their dispersal from the nesting beach that serves to quickly transpose them into deeper waters (Wyneken and Salmon, 1992). In the small number of studies conducted on neonate leatherbacks, they have been shown to be efficient dispersers during their frenzy stage, with a smaller aerobic scope than other species (Jones et al., 2007), steady swimming speeds, and moderate metabolic rates (Wyneken and Salmon, 1992;Jones et al., 2007;Hoover et al., in press).
Consequently, variability in ocean currents that result in displacement of neonates from their intended path may require different behavioral strategies in response, potentially having long-term carry over effects on their movement and distribution in older life stages (Mansfield et al., 2009;Benson et al., 2011;Shillinger et al., 2012b;O'Connor et al., 2014). Two main hypotheses on long-term carry over effects with respect to neonate and juvenile turtle dispersal have been considered previously: (1) successful dispersal of neonates/juveniles by ocean currents determines where they will forage as adults (Hays et al., 2010;Gaspar et al., 2012;Scott et al., 2014;Lalire and Gaspar, 2019), and (2) successful dispersal of neonates/juveniles into ocean currents determines population abundance at the natal beaches, as seen with regional variation in nesting densities on beaches close to favorable ocean currents (Okuyama et al., 2010;Putman et al., 2010aPutman et al., , 2011Shillinger et al., 2012b;Putman, 2018;DuBois et al., 2020). The impact of local temporally and/or spatially heterogeneous current structures on the behavioral swimming response of dispersing neonates has not been empirically researched for leatherback sea turtles.
In this study, we collected and analyzed empirical data on Caribbean leatherback neonate dispersal during two years in relation to environmental conditions to address the following hypotheses: (1) there are inter-annual differences in current speeds and thus, neonate swimming speeds, (2) there is spatiotemporal variability in environmental drivers of current speeds and neonate swimming speeds, and (3) that neonate body sizes impact their current-corrected or "true" swimming speeds. Specifically, we acoustically tracked and quantitatively analyzed Caribbean leatherback neonate movement tracks using first passage time and random walk models. We then modeled the effect of the tidal and lunar cycle and changes in current speeds over different times and years in relation to neonate true speeds and tested for significant relationships between neonate true swimming speeds and their body size. These results provide important information for designing future neonate tracking studies, distribution modeling, and informing coastal conservation efforts.

Neonate Collection and Study Site
Fieldwork was conducted in waters directly off the shore from Pacuare Nature Reserve in the Limón Province, Costa Rica ( Figure 1A). Data from Hoover et al. (in press) collected in 2016 was combined with new data collected in 2018. Atlantic leatherback neonates (n = 94) were obtained from eggs that had been collected from multiple in situ nests in 2016 (n = 42) and 2018 (n = 52). These eggs were then either reburied in a hatchery (n = 60) in Pacuare Nature Reserve or hatched in incubators near the hatchery (n = 34), following the protocol outlined in Williamson (2018). Before trials began, neonates were weighed and standard morphometric measurements (standard carapace length and width, flipper length, and head width) were taken. Neonates were chosen based upon availability and body condition, with only individuals showing healthy body conditions and activity used. Neonates were then tagged an FIGURE 1 | Maps of neonate tracks off the shore of Pacuare Nature Reserve, Costa Rica (A) for 2016 (B) and 2018 (C), with the reserve location indicated with a red box. Maps of true neonate tracks (red), were plotted using their heading and true swimming speeds, and raw neonate tracks (black). The black star shows the hatchery where turtles would emerge naturally from. The basemaps for this figure were created using ArcGIS R software by Esri. ArcGIS R and ArcMAP TM are the intellectual property of Esri and used under license, copyright © Esri, all rights reserved. hour prior to release, using Vemco V5-180 kHz transmitter tags that were attached to a 2 m line-float tether set-up and temporarily fixed to the neonates' carapaces with Vetbond adhesive, following the methods outlined in Hoover et al. (2017) (Supplementary Figure S1).

Acoustic Tracking
Acoustic tracking was conducted from a small (6 m) powerboat from August to September in 2016 and from July to August in 2018, time periods which reflected peak emergence for neonates and were chosen based on accessibility to the field site. Tracking began around 400 m from the shore to be outside of the surf zone, where neonates could not be tracked with the boat and where high predation occurs (Gyuris, 1994;Pilcher et al., 2000;Wyneken et al., 2000;Wilson et al., 2019). Tracking took place in front of the hatchery to mimic as realistically as possible the path the neonates would have taken if dispersing naturally from the hatchery. Tracking was conducted during daylight hours (0600-1600 h), beginning at different times and tidal states. Although neonates often emerge at night, they are also known to depart during the daytime and previous neonate tracking and observational studies have been performed successfully during the day (Frick, 1976;Lohmann and Lohmann, 1992;Okuyama et al., 2009;Scott et al., 2014). Additionally, due to neonates being in their continuous swimming "frenzy" stage during the entire period of our study, we did not anticipate any difference in swimming behavior during the day versus at night. After tagging, each neonate was brought onto the boat and released from the same location (∼400 m offshore, outside of the surf zone, and directly in front of the hatchery, with some variation depending on where the hatchlings were placed in the water relative to the boat) and tracked individually with a Vemco VR100 acoustic receiver and a VH180-D-10M directional hydrophone for as long as weather conditions and neonate availability permitted, generally 30-120 min. Releasing the neonates past the surf zone ensured that the swimming trajectories of the neonates were not directly influenced by nearshore waves and the bottom topography, but instead by surface currents. The boat remained behind neonate turtles at a distance of approximately 20 m, and all positions were recorded with the receiver. At 15-minute intervals along each track, a Holdpeak 866B hand-held digital anemometer was used to collect wind speed. A Barska 7 × 42 WP range-finder monocular and Plastimo Iris 50 hand-held compass were used to determine the distance and bearing of the neonates from the boat to ensure accuracy of recorded neonate positions. At the end of the tracking period for each neonate, the turtle was recovered, the acoustic transmitter and tether removed, and the turtle was released.

Surface Drifters
Surface drifters were constructed following a similar method to Hoover et al. (in press). These drifters had over two meters of line connected to a metal drogue and a bright orange flag on the top that increased visibility but was also heavy enough to lay horizontally on the surface; the drifters were tracked using a GPS-enabled phone (Supplementary Figure S2). To minimize the wind-induced movement of the drifters (slippage), the drifters had a large submerged to unsubmerged ratio (Niiler and Paduan, 1995;Lumpkin and Pazos, 2007). These surface drifters were deployed every 20-30 min along the tracks, when the weather permitted, to determine surface current speeds.

Processing of Acoustic Tracking Data
All data processing and analysis was conducted using R statistical software (R Core Team, 2018). Each neonate turtle and drifter track were divided into regular 5-minute position intervals. Neonate and drifter data were then merged on their matching time stamps for each position. Neonate positions were corrected using bearing and distance of the neonates from the boat for each 5-minute time point of the tracks. Raw track speed of both the neonates and drifters was found by taking the distance traveled between each point divided by the time duration between points. The "geosphere" R package (Hijmans et al., 2019) was used to find the bearing of each neonate and drifter point.
True neonate swimming speed (termed "in-water swimming speed" by Hoover et al., in press), which accounts for the effect of the currents on the neonate swimming speed, was found by using the raw track speed of the neonates (S N ) and the heading in radians ( ) to determine the neonate east-west (u N ) and north-south (v N ) speed components (1-2) (Bailey and Thompson, 2010;Hoover et al., in press). The drifter eastwest (u D ) and north-south (v D ) speed components were found using similar equations, but using the raw track speed and direction of the drifters instead of the raw track speed and direction of the neonates. The final true speed (V N ) for each neonate position point was found using equation (3). Neonate true headings ( T ) were derived from taking the arc tangent of the absolute difference in neonate east-west and northsouth speed components with drifter north-south and east-west components (4).
It is worth noting that in computing the true headings of the neonates, wave, wind, and inertial influences could differentially impact the neonates and the drifters. However, inertial influences are unlikely to be important, as the mean size of the neonates gives them an effective radius that minimizes inertial effects to less than 5 • of the deflection angle (Beron-Vera et al., 2015;Brooks et al., 2019) and at the modest wind speeds measured (1-5 m s −1 ), assuming a 0.27% windage effect calculated from a similar drifter (Meyerjürgens et al., 2019), the error associated with the windage would be (0.0027-0.0135 m s −1 ), which is orders of magnitude smaller than the neonate swimming speeds. Wave influences, or Stokes drift, would be aligned onshore in these nearshore waters and perpendicular to the isobaths and would potentially hinder the neonate's offshore trajectory and slow their swimming speed. However, recent work in the Gulf of Mexico suggests that the summertime Stokes drift has a sharp maximum in probability distribution in inshore waters of depths less than 100 m at 0.01 m s −1 (Clark, 2015), and even in hurricane conditions the median drift was only 0.03 m s −1 . Thus, over the short time scales of these experiments, the drifter and neonate motions are only very weakly influenced by these effects and this influence is within measurement error.

First Passage Time Analysis
To determine where neonates were spending more time along their track within a given spatial radius, a first passage time (FPT) analysis was conducted. The R package "adehabitat" (Calenge, 2006) was used to calculate the FPT in hours for each neonate position. After accounting for a 5-minute acclimation period at the beginning of each track, the minimum and maximum step length for each of the tracks was used to develop a range of radii. For each neonate track, the log-transformed variance in FPT was then calculated and used to identify the peak in variance and associated radius of interest. The mean radius was calculated by averaging the radii of interest across tracks, and this overall spatial radius was then used to determine the FPT for every position within each individual track.

Random Walk Models
Random walk models can be used as a null hypothesis of movement, as they assume an individual is moving through its environment randomly with independent distributions of move lengths and angles and no correlation between current and previous positions. The net squared displacement (NSD), or the squared distance between each position along a track and the first position of the same track, can be used to determine whether an individual is following a correlated random walk (each position is correlated with the previous position), or biased random walk (each position has a directional bias). To test whether each of our individual neonates were moving randomly in their dispersal phase or following a fixed compass direction (biased random walk) as they left the shore, the net squared displacement (NSD) was calculated for each current-corrected position for each neonate (see section "Processing of Acoustic Tracking Data"). The NSD's were found for each of these true neonate tracks by taking the distance between each position along a track and the first location and squaring it. NSD for true tracks was then compared to the NSD expected for a correlated random walk (CRW) and a biased random walk (BRW). The expected CRW and BRW NSD's were calculated using the following equations (Turchin, 1998): Where n is the number of moves from the first location, m 1 is the mean step length, m 2 is the mean squared step length, θ is the mean cosine of absolute move direction, and c is the mean cosine of the turning angles (the angle between sequential moves).

Environmental Drivers of Current Speed Variability
To determine environmental drivers of current speed, and therefore any inter-annual or spatial differences that could impact neonate dispersal, a generalized additive mixed model (GAMM) was fit to the 2016 and 2018 data. This model predicted the speed of the local currents (derived from the raw track speed of the drifters) as a function of local wind speed and the categorical variables of lunar cycle phase (classified as spring tide conditions or full and new moon phases; neap tide conditions or quarter moon phases; and mid-tidal conditions or 3 days before/after full, new, and quarter moon phases), time of day, year, and tidal current direction (classified as ebb, flood, or slack, with slack being within 15 min of high/low water). Tidal data were obtained for the surrounding region of Limón, Costa Rica (10.00 • N, 83.03 • W), which experiences a mixed semi-diurnal tidal cycle (from Mobile Geographics) LLC ©1 . A smoothed interaction function of latitude and longitude was included to account for spatial autocorrelation present in current speed. Exploratory analysis using backward model selection and F-tests to remove non-significant terms and compare terms with and without the non-linear smooth of latitude and longitude indicated the presence of non-linear relationships, indicating a GAMM approach was appropriate. Individual drifter tracks were 1 https://tides.mobilegeographics.com/locations/4273.html included as random effects to account for correlation in current speed within each track. All GAMMs were fit in R using the "mgcv" package (Wood, 2011). The error distribution for the current speed model was Gaussian with an identity link function, and the smoother terms for the interaction of latitude and longitude were fit using penalized regression splines. Backward stepwise model selection was used to select the final model, using Akaike information criterion values (AIC) and only retaining variables with significant p-values (p < 0.05). Model assumptions were checked using a plot of the standardized residuals versus fitted values and auto-correlation function (ACF) and partial auto-correlation function (PACF) plots of the standardized residuals.
To determine the environmental drivers of nearshore current flows measured by the drifters and if they were impacted by large-scale geostrophic current flows, geostrophic current velocity was calculated using extracted mean values for 2016 and 2018 to explore variations in external forcing in each year. For geostrophic current velocity, the delayed time, global altimetry dataset from the AVISO + Live Access Server (LAS 2 ) was used to create plots for the weeks of the 8th August 2016 and the 6th August 2018 for absolute geostrophic vector velocity for the grid area 88 • to 75 • W and 7 • to 19 • N. These dates were chosen based on availability from the LAS dataset and that the second week of August was a shared week of neonate tracking for both 2016 and 2018; as these plots reflect large-scale geostrophic currents, the week chosen for the dataset did not change the velocity patterns observed.

Generalized Additive Mixed Model for Neonate True Speed
We fit a GAMM to examine how neonate true speed was influenced by FPT (derived from the neonate raw tracks), local current speed and direction (derived from the drifters raw track speed and direction), tidal current direction, lunar phase, wind speed, year, latitude and longitude, and time along the track. The GAMM was fit after using a backward model selection exploratory analysis to compare term significance with and without non-linear smooths to determine that nonlinear relationships were present. Current speed, FPT, and the interaction of latitude and longitude each had non-linear smoother terms. Individual neonates were included as random effects to account for correlation in true neonate speed within each track. The error distribution was Gaussian with an identity link function, and the smoother terms for the interaction of latitude and longitude were fit using penalized regression splines. The smoother terms for current speed were estimated by building departures for each individual neonate from a common smooth, using a factor smooth interaction and five basis functions (k = 5) to avoid over-fitting. This allowed for random curves to be built for each neonate, using splines that share the same smoothing parameters and an un-penalized null space. Model selection was performed in the same manner as for the current speed GAMM.
To rule out whether neonates were responding to changes in the environmental conditions they encountered or simply to natural variation in behavior from the time spent swimming, we ran an additional linear mixed effects model predicting the variance in swimming speed for the neonates as a function of time along their track, with individual neonates included as random effects and controlling for the speed of the currents. The variance in swimming speed was extracted as residuals from the final chosen GAMM model fit for the neonates and the time along their track was represented by move number, as each track was broken into regularized 5-minute time steps.

Body Morphometrics and Neonate True Speed Analysis
A principal component analysis (PCA) was used to derive an index of neonate body size, using the body morphometric variables which demonstrated covariance (carapace length, body mass, carapace width). The first principle component was extracted and used as a neonate size index (Burgess et al., 2006). This index was used to compare variability in body size between years and nest treatments (incubator versus hatchery nests), using independent t-tests.
A linear mixed effects model was fit in R using the "nlme" package (Pinheiro et al., 2020) to determine if there was a relationship between neonate body size, represented by the PCAderived size index, and their average true swimming speed.
The year and treatment of origin (incubator versus hatchery nest) were used as random effects in the model after a stepwise model comparison. Neonate swimming speed was standardized for the body length of the neonates by dividing their average true swimming speed by their carapace length. Model assumptions of normality and independence were checked using residual plots.

Acoustic Tracking
A total of 27 neonates tracked from hatchery nests and 15 from incubator-hatched eggs in 2016, and 33 neonates tracked from hatchery nests and 19 from incubator-hatched eggs in 2018. Drifters were deployed for the majority of neonate tracks (n = 76) to determine current speed and direction. Average values for track distance and duration, current speed, raw track and true neonate speed, current direction, and raw track and true heading of all neonates in each year are reported in Table 1, demonstrating any inter-annual differences in current speeds or neonate trajectories. Current speeds in 2018 were higher on average than in 2016, but neonate true speed was similar during the 2 years. The mean current directions and raw track headings of the neonates had a strong south-east direction both years, but neonate true headings were in a more north-east direction and similar between years (Table 1 and Figures 1B,C).

First Passage Time
A radius of 59 m (±57 SD, range 5-228 m) was the mean spatial scale of interest for the neonate tracks (n = 94). Plots of the FPT over time for each track at this scale indicated that the majority of tracks had a higher FPT in the first third of the track (Figure 2A and Supplementary Figure S3), even when a 5-minute acclimation period was accounted for. Lower FPT values coincided with higher values for true speed of the neonates (Supplementary Figure S4). The higher values for FPT observed in the first third of each track reflected some spatial variability in FPT, with areas nearshore having higher FPT values and areas farther from shore having lower FPT values (Supplementary Figure S3).

Random Walk Models
All neonate true tracks during 2016 and 2018 followed a BRW during the first five move numbers, with an NSD following that expected with a BRW (example track shown in Figure 2B), indicating that the neonate current-corrected positions had a directional bias consistent with the movements of an animal with a fixed compass direction. However, plots of NSD for true tracks showed they had a greater NSD than predicted for a BRW with move numbers higher than five (Figure 2B), indicating that they were moving quicker than expected away from their release point with directed movement.

Inter-Annual Variability and Environmental Drivers of Current Speeds
Drifter-derived current speed (41 drifter tracks in 2016 and 93 drifter tracks in 2018) had a south-east orientation and significant inter-annual and spatial variability ( Table 2). Current speeds in the 2018 study period were significantly higher than 2016 (p-value < 0.0001, Supplementary Figure S5). The significant interaction smoother for latitude and longitude in the GAMM showed currents had significant spatial variation in both years (p-value < 0.0001, Table 2), with lower current speeds nearshore (∼300-1,000 m from shore) and higher current speeds farther from shore (∼1,800-3,800 m from shore) ( Figure 3B).
Nearshore drifter-derived current speed and direction was broadly consistent with large-scale geostrophic current flow, suggesting that the inter-annual variations in large-scale current flow are relevant to nearshore current processes in our tracking area. For the geostrophic current velocity, there were stronger currents in 2018 than in 2016 within the spatial area of the obtained data within the south Colombia Basin (Figure 3A), with a higher mean absolute vector velocity of 0.646 m s −1 in 2018 compared to 0.578 m s −1 in 2016.

Environmental Variables and Neonate True Speed
Significant spatial autocorrelation for neonate true speed was seen in plots of the residuals, indicating the need for an autoregressive correlation structure. A model selection process using AIC values revealed that the best model was neonate true speed as a function of year and non-linear, smoothed variables of current speed and FPT (Table 3), with an autoregressive-movingaverage (ARMA) lag 1 correlation structure. Neonates swam at a steady speed until current speeds reached 0.5 m s −1 . A plot of the smoother for neonate true speed versus the current speeds showed three different response strategies to high current speed conditions (>0.5 m s −1 ). Most individuals had speeds that increased rapidly in positive association with the currents (n = 41), while others maintained either a constant swimming speed (n = 14) or had a negative association (n = 21) between their swimming speed and the currents ( Figure 4A).    Current speed was fit with a basis dimension of k = 5 and factor smooth interactions, with separate smoother terms for current speed for each individual track.A correlation autoregressive-moving-average (ARMA) structure was added with p = 1 and q = 1.Overall deviance explained was 35% edf: estimated degrees of freedom.
Neonates that increased their swimming speed in response to high current speed conditions had tracks with directions that deviated in a northerly direction from the currents, while those with a decrease or constant response strategy had tracks that were more closely matched with the direction of the currents ( Figure 4B). This was further evidenced with an increase in current speed resulting in a greater difference in neonate north-south and east-west vector components compared to current vectors ("Mean Difference in East-West/North-South Vectors, " Table 4). The neonates' north-south vector components were stronger when currents were faster, indicating that they were swimming in a more northerly compass direction in faster current speeds (  Figure S7).
Even though individual variability was present in responses to high current speed conditions, neonates at the population level maintained a consistent mean swimming speed between years, with a mean speed of 0.27 m s −1 in 2016 and 0.25 m s −1 in 2018 ( Table 1). The significant interaction smoother for latitude and longitude in the GAMM indicates neonates had spatial structure to their swimming speed ( Table 3). This reflected the spatial structure seen in the current speeds, with neonates swimming faster in areas farther from shore ( Figure 3C). This was further confirmed with the smoothed relationship of their true speed as a function of FPT (Supplementary Figure S6), with neonates having higher swimming speeds with lower FPT values, which correlated with areas farther from shore with faster current speeds (Supplementary Figures S5, S3B).

Body Size and Neonate Swimming Speeds
The PCA results showed that mass and carapace length had the largest correlations and loadings with the first principal component, which explained 56% of total variance and was loaded with carapace length, width, and body mass (Supplementary Figure S8). Scores from the first principal component were extracted for each turtle and used as a size index (Burgess et al., 2006). Body size measurements for neonates appear in Supplementary Table S1. An independent t-test showed a significant difference in the size index for neonates between years, with 2018 having a larger mean size index than 2016 (2018: 0.501; 2016: −0.453, p-value = 7.95e-04) (Figure 5A). Independent t-test results showed that neonates from hatcheries had significantly larger body size indices on average than incubator neonates in FIGURE 4 | (A) Generalized additive mixed model (GAMM) smoothing curve for current speed (m s −1 ) as a function of neonate true swimming speeds for 2016 and 2018, with an estimated random departure from the mean true speed trajectory for each neonate represented by colored lines and k = 5 basis functions. Each color represents the response strategy of each neonate to high current speed conditions (>0.5 m s −1 ): mean increase in true speed (green), mean decrease in true speed (yellow), or no change in mean true speed (purple). The threshold for high current speed conditions is shown with a black dashed line. (B) Rose plots for the frequency of mean neonate true swimming speed plotted by neonate raw track bearing or current direction. Colors represent neonate raw track bearing for each swimming speed response strategy in currents greater than 0.5 m s −1 (decrease/yellow, constant/purple, increase/green), and black represents current direction for current speeds greater than 0.5 m s −1 .
The linear mixed effects model did not show a significant relationship between average neonate swimming speed and neonate body size (Table 5 and Figure 5C).

DISCUSSION
Dispersal of leatherback neonates from Pacuare, Costa Rica is significantly affected by current speed, with spatial and temporal variability in the currents resulting in three different behavioral response strategies in neonate swimming speed and 4 | Mean and standard deviation values (in parentheses) for true neonate swim speed, current speed, current direction, neonate raw track and true heading, mean difference in north-south and east-west vector components for neonates and the current, and mean north-south and east-west vector components for neonates for different current speed categories (low < 0.3 m s −1 , medium 0.3-0.5 m s −1 , high > 0.5 m s −1 ). orientation. This could have implications for individual neonate fitness and survival, as changes in strength and direction of the currents could impact the energy use of neonates (Jones et al., 2007;Shillinger et al., 2012b), the amount of time spent in predator-rich waters (Gyuris, 1994;Pilcher et al., 2000;Bolten, 2003;Stewart and Wyneken, 2004;Santidrián Tomillo et al., 2007), and potentially have carry-over effects into their long-term migratory pathway and ability to find resource patches offshore in older life stages (Gaspar et al., 2012;Shillinger et al., 2012b;Putman and Mansfield, 2015;Briscoe et al., 2016;Cardona and Hays, 2018). Although we identified significant inter-annual differences in current speed and neonate body size, the larger sizes of neonates in 2018 did not significantly impact their swimming speed, showing that variability in average neonate swimming speed was due to environmental, not biological, conditions. These findings are especially important for the future conservation and management of this endangered and declining subpopulation of Caribbean leatherback turtles (Troëng et al., 2007;Stewart et al., 2016;The Northwest Atlantic Leatherback Working Group, 2019), as the faster currents experienced during anomalous years may either enhance neonate dispersal or require them to have higher energy expenditures to maintain their heading away from shore. Our findings substantiate that leatherback neonates swim at a mean constant speed in their frenzy stage to efficiently move offshore and out of high-risk, nearshore zones (Wyneken and Salmon, 1992;Wyneken et al., 1998;Jones et al., 2007;Okuyama et al., 2009), areas where we observed higher FPT's and slower current speeds. Neonates had similar mean true speeds between years and in different current speed conditions (low, medium, high speeds) ( Table 4), despite significantly faster current speeds in 2018 versus 2016. On average, as current speeds increased, neonates increased the difference in their directional components with the currents, indicating they were swimming increasingly in the opposite direction of the currents. An increase in current speed also resulted in most of the neonates swimming more strongly, on average, into a given compass heading opposite of, and in a more northerly direction to, the currents. This reinforces the theory that neonate sea turtles are following a fixed compass direction as they disperse from shore (Lohmann, 1991;Lohmann and Lohmann, 1996;Goff et al., 1998).
Despite relatively constant mean true swimming speeds between years (Table 1) and in different current conditions at the population level (Table 4), there was individual variability in neonate swimming speed in response to faster current conditions, with three main response strategies observed and no indication that these behaviors were simply due to time spent swimming (Supplementary Figure S7 and Supplementary Table S2). When currents surpassed speeds of 0.5 m s −1 , more than half of the individuals (54 percent) increased their true speed in positive association with the current velocity and in an oriented, more northerly direction from the currents (Figure 4). Increased velocity in high current speed conditions (>0.5 m s −1 , see Figure 4A) may have negative implications for the neonates' energy storage, as leatherback neonates have a metabolic scope and morphometrics optimized for continuous, slow swimming and a limited yolk storage that is needed to sustain them in their long journey offshore to find foraging areas (Davenport, 1987;Wyneken et al., 1998;Jones et al., 2007). The other two response strategies observed were: (1) neonates maintaining a constant mean, oriented true swimming speed around 0.25 m s −1 , or (2) decreasing their true speed in high current speed conditions (>0.5 m s −1 ) but having raw track bearings more closely matching the direction of currents ( Figure 4B). Neonates that utilize either a constant or decreased swimming speed strategy in faster currents father from shore could conserve yolk storage and have higher chances of survival as currents transport them to larger current systems offshore in a shorter amount of time (Shillinger et al., 2012b) and decrease the amount of time spent in predator-rich areas (Wyneken and Salmon, 1992;Wyneken et al., 1998;Okuyama et al., 2010;Santidrián Tomillo et al., 2010). This was further suggested by the lower first passage times and raw track speeds for neonates in areas farther from shore (>1,800 m from shore), where the currents were faster (Supplementary Figure S5 and Figure 3B), and their ability to more efficiently displace themselves than predicted by a biased random walk further into their track (Figure 2B).  Year and origin (incubator or hatchery) were included as random effects.
Overall, the variability in individual swimming behavior in high current speeds seen in our results may either boost dispersal and/or incur metabolic costs for this decreasing subpopulation of Caribbean leatherbacks (The Northwest Atlantic Leatherback Working Group, 2019). The optimization of leatherback neonates for continuous, slow swimming (Jones et al., 2007) raises concerns that 53 percent of our neonates appeared to increase their swimming speed under high current speed conditions (Figure 4A), despite a population-level response of relatively constant swimming effort in different current speed categories ( Table 4). These individual-level responses to extreme current speed conditions could possibly prematurely exhaust valuable yolk reserves that would otherwise be used to reach offshore foraging habitats. However, in the absence of these potential metabolic costs, the diversity in neonate swimming responses to high current speeds may increase their survival and aid their dispersal offshore. These fast currents could transport neonates away from dangerous nearshore waters where predation is high (Gyuris, 1994;Pilcher et al., 2000;Wyneken et al., 2000;Wilson et al., 2019) and potentially result in neonates having a wider distribution as they encounter variable current systems offshore, increasing their population resiliency to environmental change (Mansfield et al., 2017).
The combination of statistical models and inter-annual data collected on current speeds with environmental predictors show currents impacting neonate movements in this area are not driven by tidal current direction, time of day, lunar cycle phase, or local wind speed. Nearshore current flow spatial and temporal structure in this region is instead likely driven by changes in largescale geostrophic currents within a cyclonic gyre in the southern Colombia Basin (Centurioni and Niiler, 2003). The inter-annual differences in nearshore local currents measured by our drifters reflected inter-annual differences in geostrophic currents in this region, with faster mean geostrophic current speeds during our 2018 tracking period than in 2016 ( Figure 3A). The impact of large-scale geostrophic current flow in the Colombia Basin on local nearshore current flow and thus neonate dispersal supports that future variability in juvenile and adult dispersal behavior in oceanic gyres may depend on their swimming behavior responses as neonates to these environmental conditions (Shillinger et al., 2012b;O'Connor et al., 2014).
Higher precipitation in 2018 may have affected the incubation duration and body morphometrics of neonates in this study (11.96 mm day −1 for 2016 and 15.57 mm day −1 for 2018, see Supplementary Figure S6). Higher precipitation may result in decreased nest temperatures and longer incubation durations (Lolavar and Wyneken, 2015), which subsequently produce larger body sizes but decreased locomotor abilities and overall fitness (Booth et al., 2004;Fisher et al., 2014;Booth, 2017). This is likely due to a decreased fatigue resistance and inability to power stroke for long periods (Burgess et al., 2006). Although we were unable to retrieve precise data on the incubation duration for individual neonates from hatchery nests, the reported average for all neonates from the hatchery nests in 2018 was 70 days (Williamson, 2018, personal communication), which is longer than the known average of 60 days for leatherbacks (Bilinski et al., 2001). Additionally, neonates in 2018 were significantly larger than in 2016 (Figure 5A), and of these, eggs hatched in the hatchery, where nest temperature was not controlled and vulnerable to precipitation effects, were significantly larger than those hatched in incubators ( Figure 5B). Body mass for neonate leatherback neonates is usually 30-47 g (Mickelson and Downie, 2010), a range which was greatly exceeded by more than 75% of the neonates in 2018 (range: 42-56 g). However, we did not find a significant relationship between neonate body size and average true swimming speed ( Figure 5C). These results show that environmental variables, such as current speed and direction, may have a greater impact on neonate swimming behavior for this Costa Rican rookery than biological factors, such as neonate size, incubation duration, or incubation treatment. However, to fully test the impact of each of these biological variables on their true or current-correct swimming speed, fitness tests would need to be performed on individuals.
Future research efforts that can inform conservation for this population should incorporate longer tracking durations through the deployment of fixed acoustic telemetry arrays in nearshore neonate dispersal corridors highlighted in this study. Use of miniature satellite tags on larger juveniles or other age classes would extend our understanding of the potential carryover effects of environmental conditions on Caribbean leatherback turtle distribution, their ability to successfully find foraging areas, and the risk of interaction with fisheries (Shillinger et al., 2012a;Putman and Mansfield, 2015;Christiansen et al., 2016;Mansfield et al., 2017). The long-term probability of neonate survival and dispersal success leaving the nesting beach during different climatic and oceanographic conditions could be assessed using model simulations with particle tracking techniques, ocean current models, and empirically-informed parameters on neonate movement, growth, and metabolism (Putman et al., 2011;Scott et al., 2012;Shillinger et al., 2012b;Briscoe et al., 2016;Christiansen et al., 2016;Gaspar and Lalire, 2017;Lalire and Gaspar, 2019;DuBois et al., 2020).
The methodological approach used in this study allowed us to successfully track a large number of sea turtles during a vulnerable life stage that has previously had its movement ecology characterized mainly by observational or laboratory studies (Frick, 1976;Lohmann et al., 1990;Wyneken et al., 1990;Lohmann and Lohmann, 1992;Goff et al., 1998;Smith and Salmon, 2009). The methods used allowed us to identify three different neonate behavioral response strategies in high current speed conditions and these results may have long-term implications for neonate energetics, survival, distribution, and overall population success in the face of varied environmental conditions at sea. This methodological approach could be applied to the neonate stage in other marine species, which is understudied in movement ecology due to the difficulty of tag attachments (small body size relative to equipment) and often high mortality rates (Shillinger et al., 2012a). Juvenile movements are generally better understood, such as for sturgeon, where studies have used biotelemetry with acoustic tags to look at home-ranges, seasonal and ontogenetic habitat associations, first passage times, and movements predicted by environmental covariates (Geist et al., 2005;Smith and King, 2005;Thomas et al., 2019). These and other highly migratory marine taxa, such as sharks, billfish, and swordfish, may take advantage of either the affordable miniature acoustic tag and drifter design used in this study and/or the combined analytical approach utilizing statistical models to provide an integrated interpretation of the relationship between an individual's movement and its environment during dispersal. The characterization of these movements of younger life stages can then be used to inform management, conservation efforts, and responses to climate change scenarios, as has been performed for older life stages (e.g., Hazen et al., 2012b).

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

ETHICS STATEMENT
The animal study was reviewed and approved by University of Maryland Center for Environmental Science Institutional Animal Care and Use Committee.