Juvenile hawksbill residency and habitat use within a Caribbean marine protected area

Understanding the spatial ecology of highly mobile marine vertebrates is necessary for informing conservation and management strategies aimed at protecting such species. Buck Island Reef National Monument (BIRNM), off the coast of St. Croix, US Virgin Islands, harbors critical foraging habitat for Critically Endangered juvenile hawksbills Eretmochelys imbricata that exhibit high site fidelity until sexual maturation. Using an array of fixed passive acoustic receivers that covered over 20.2 km2 at its largest configuration and in-water biannual sampling, we analyzed residency patterns and habitat use of 29 hawksbills. High recapture rates allowed for longterm data collection for some individuals, with 1 individual being detected within the array 1952 d (mean ± SD: 411 ± 444 d). We used detection data to construct a resource selection function based on a generalized linear mixed model in order to determine relative habitat selection, or the use of different habitat types within an area proportional to the ‘true’ selection. Our covariates in the model were benthic structure, bathymetry, time of day, and year. Results showed selection by tagged individuals for shallow (<20 m) depths in areas with high rugosity characterized by a high density of reef or rock. However, individuals also selected areas comprised primarily of sand interspersed with seagrass pastures. We also used the best supported model to predict relative selection across BIRNM and found that the total area of high relative selection decreased significantly at night. The information provided will help guide both future in-water surveys for cryptic hawksbills and management decisions about public area-use at BIRNM.


INTRODUCTION
Coastal marine ecosystems across the globe are threatened by direct and indirect anthropogenic influences (Halpern et al. 2015). Similar to terrestrial reserves, marine protected areas (MPAs) have become a popular conservation tool to mitigate detrimental influences within designated areas (Keller et al. 2009). Designation and management of MPAs often takes into consideration large, charismatic mar-ine vertebrates despite the inherent difficulties of protecting such species using place-based management methods (Browman et al. 2005, Keller et al. 2009). One argument for this is the recognition that marine megafauna can have important influences on coastal ecosystem communities, typically through top-down processes (Heithaus 2013, Hays et al. 2016. Understanding how and the extent to which such species use a protected area (including preference for areas outside of an MPA) can help better inform management decisions as well as provide insight into the state of other protected resources.
Distributed globally throughout the tropics, hawksbill sea turtles Eretmochelys imbricata are one such species of marine megafauna shown to directly impact coral reef ecosystems through their consumption of macroalgae and sponges (Meylan 1988, León & Bjorndal 2002, Gaos et al. 2012, Bell 2013. However, habitat degradation, direct take, and incidental bycatch over the past century have greatly reduced global hawksbill populations, causing the IUCN to list the species as Critically Endangered (Mortimer & Donnelly 2008). Similar to other marine turtle species, hawksbills exhibit ontogenetic shifts over long lifespans wherein they use both pelagic and coastal environments (Musick & Limpus 1997, Godley et al. 2008. Post-hatchling hawksbills spend a variable number of years foraging at or near the surface in drifting oce anic weed lines during the 'lost years' before actively recruiting to neritic habitats and switching to a demersal foraging behavior (Musick & Limpus 1997). Juveniles remain resident in these developmental foraging areas until sexual maturation (between 10 and 20 yr of age), after which they typically exhibit more migratory behavior between mating, nesting (for females), and foraging grounds. The transitory nature across relatively large and often remote areas makes place-based protection for posthatchling and adult hawksbills unrealistic. In contrast, the juvenile develop mental stage offers a practical opportunity to protect hawksbills using MPAs.
Designing and effectively managing MPAs for hawks bills requires knowledge about their spatial ecology. A variety of techniques have been used to investigate the spatial ecology of juvenile hawksbills within the Caribbean region. Several studies using in-water capture-mark-recapture (CMR) have provided important information on the population demographics, residency time, and horizontal displacement between capture events of juvenile hawks bills (van Dam & Diez 1998, Blumenthal et al. 2009a, Hart et al. 2013, Strindberg et al. 2016. Active tracking, whereby a tagged individual is followed from a boat using a directional receiver, produced preliminary estimates of home range, whereas time− depth re corders provided data on diurnal dive patterns (van Dam & Diez 1997, 1998, Blumenthal et al. 2009b, Witt et al. 2010, Scales et al. 2011, Berube et al. 2012, Wood et al. 2017. Results from these studies indicate long-term residency patterns of juvenile hawksbills within coastal Caribbean habitats. However, most previous studies had small sample sizes (<10 individuals), short tracking durations (< 25 d), limited spatial data (< 20 unique locations), or some combination of these. Passive acoustic telemetry is increasingly being used to track larger numbers of individuals while mitigating some of the problems outlined above inherent to these other techniques.
Fixed passive acoustic telemetry provides continuous, fine-scale location data for numerous individuals over extended periods of time (years) with less potential for influencing animal behavior compared to active tracking (Heupel et al. 2006, Kessel et al. 2014, Hussey et al. 2015, Zeh et al. 2015. Fixed passive acoustic technology combines the portable hydrophone and receiver used for active acoustic tracking into a single submersible unit referred to simply as a receiver. Receivers anchored to the substrate at predetermined locations record the date, time, and tag ID (collectively termed a detection) of any acoustic tag within the receiver's range. The technology is used extensively to understand the movement of numerous fish species, but only re cently has it been used to examine marine turtle spatial ecology (Hazel et al. 2013, Chevis et al. 2017. The array at Buck Island Reef National Monument (BIRNM), US Virgin Islands, is unique in both its coverage and longevity and was made possible through the continued collaboration of academic and government (both local and federal) researchers. These factors provided us the opportunity to quantify previously suggested behavior patterns of juvenile hawksbills and establish a baseline for future comparisons. Higher-resolution spatial and temporal habitat-use data provide more power to detect anomalous changes in the environment through shifts in behavior. This is especially important in the highly connected and dynamic marine environment, where adverse anthropogenic threats far from an MPA boundary, such as an oil spill, or regionally imposed threats, such as climate change, can have detrimental effects on those areas protected by the MPA over time.
Using this fixed passive acoustic array of receivers, we investigated the residency and habitat-use trends of juvenile hawksbills captured within the protected waters of BINRM. Specifically, our objectives were to (1) examine residency through temporal detection patterns of tagged hawksbills within the array, and (2) investigate space use within BIRNM using a resource selection function (RSF) framework to understand relative selection based on 4 covariates: benthic substrate, bathymetry, time of day, and year. Results from our study will give resource managers insight into areas of potential importance for juvenile hawksbills at BIRNM.

Study site
Designated in 1961, BIRNM originally encompassed 3.6 km 2 , including the uninhabited Buck Island, until the boundary was expanded in 2001 to 73.4 km 2 . Managed by the US National Park Service, BIRNM is one of only a few no-take marine protected areas in the Caribbean and harbors a diverse assortment of habitats (Guarderas et al. 2008). An emergent reef crest surrounds Buck Island, creating a shallow 50− 150 m wide lagoon starting on the southern side and continuing counterclockwise to the northwest corner, where it breaks apart into a series of isolated patch reefs. South/southwest of Buck Island are sea grass beds comprised of Syringodium sp., Thalassia sp., and Halophila sp. interspersed with low-rugosity sand. North of Buck Island are densely clustered remnant stands of dead elkhorn coral Acropora palmata that rise to the surface from a depth of 9−15 m, called 'haystacks' (Bythell et al. 1989). The receiver locations were centralized around Buck Island and dispersed throughout the various habitats (Fig. 1). The passive acoustic array at BIRNM began as 6 Vemco VR2W receivers (Vemco Amirix Systems) deployed in September 2011 and through subsequent collaborations reached 141 receivers at its pinnacle. Receiver downloads and station maintenance were performed twice a year.

Field sampling
We began capturing and tagging juvenile hawksbills in the spring of 2012. Subsequent tagging trips lasted between 8 and 10 d and were conducted twice annually (spring and fall) until 2016. To capture individuals, 2 teams of 3 snorkelers swam designated transects and pursued any observed hawksbills from the surface until they settled on the sea floor. One team member would then free-dive to the bottom and grasp the turtle firmly by the carapace at the nuchal and anterior marginal scutes before slowly Relative location of BIRNM within the greater Caribbean region ascending with the turtle's head pointed towards the surface. We recorded morphometric measurements be fore tagging individuals with either Vemco V16-4L (16 × 88 mm, 24 g in air, 69 kHz, 152 dB, with 30−90 s delay interval) or V13-1L (13 × 36 mm, 11 g in air, 60− 84 kHz, 147 dB, with 30−90 s delay interval) acoustic tags depending on the individual's size; no tag exceeded 0.8% of an individual's body weight in air. We attached acoustic tags to the flattest section of posterior marginal scutes using rubber-coated wire and Devcon Marine Plastic Putty with the transmitting end left open so as not to interfere with signal transmission. All turtles were released within 1 h at the point of capture.

Data filtering
False positive detections, more commonly referred to simply as false detections, occur when tag transmissions are distorted by other nearby acoustic signals with similar frequencies (Simpfendorfer et al. 2015). Some false detections appear as new unknown tag IDs, which are easily filtered and removed. More problematic are type B false detections, which record as a known deployed transmitter that was not in actual proximity of the receiver (Simpfendorfer et al. 2015). We flagged potential type B false detections using 2 filters created in R version 3.4.3 (R Core Team 2017). The first basic filter flagged any detections that occurred before a tag was deployed or after the estimated battery life provided by the tag manufacturer (2808 d for V16s and 623 d for V13s). The second filter flagged any detections that indicated an unrealistic swim speed for juvenile hawksbills (>10 m s −1 ) between receivers more than 1 km apart (Revelles et al. 2007, MacDonald et al. 2013. Distances between receivers were calculated using a least cost pathways analysis in R with raster data of BIRNM's major habitat types and bathymetry. Land (Buck Island) and shallow (<1 m) reef were scored as highly resistant so that the resulting distance matrix be tween receivers reflected the shortest distance be tween receivers that a hawksbill could potentially swim (Zeller et al. 2012).

Residency
To investigate individual detection history within the array, we created a residency plot visually depicting when an individual was detected in the array over time. We also calculated a residency index (I R ) for each individual by dividing the total number of days detected in the array since initial tag deployment by 365 d (Abecasis & Erzini 2008, O'Toole et al. 2011. We chose 365 d from initial tag deployment to calculate I R because it gave us a year of detection history and was the maximum number of days between the last acoustic tag deployment and the latest array download. An I R = 1 indicates that an individual was detected within the array at least once every day for the 365 d period, whereas an I R = 0 implies that an individual was never detected within the array.

RSF model
RSF models are a type of point process model that compare spatial locations with randomly assigned background points within a predetermined available area for tagged individuals (Boyce & McDonald 1999). By definition, RSFs are proportional to the true probability of selection; therefore, selection strengths are considered relative to each other and not in absolute terms (Manly 2002). Thus, throughout this paper, we will refer to the quantities estimated by our RSF as 'relative selection.' The resulting binary data structure can be analyzed using a generalized linear mixed model (GLMM) with a binomial distribution where individuals' locations are coded as '1' and the random available points as '0' (Hebblewhite & Merrill 2008). Covariates associated with those locations are evaluated to determine relative selection for each habitat, with individual included as a random effect to account for individual variability and unequal sampling.
After filtering out type B false detections using the methods described above, we calculated centers of activity (COAs) for tagged individuals based on the mean position algorithm (Simpfendorfer et al. 2002). The mean position algorithm calculates a COA for an individual by taking the average of the receiver locations where detections occurred within a certain time interval, t (Worton 1987, Simpfendorfer et al. 2002. We chose a Dt of 60 min to reduce temporal autocorrelation between locations, but also to maintain as fine-scale temporal resolution as possible (Scales et al. 2011, Chevis et al. 2017. Individuals with £100 COAs were removed from the data, leaving 24 individuals for analysis. Critical to implementing an RSF model, especially with highly mobile species, is defining a meaningful scale for available habitat (Hooten et al. 2017). This is further complicated in the case of passive acoustic technology, where location information is inherently inferred from static receivers with dynamic ranges (Selby et al. 2016). We must therefore restrict our inference to areas covered by the array. Based on previous range-testing of the array at BIRNM, we chose a 400 m buffer around each receiver to represent our 'available' habitat and removed any COAs outside of the buffer region (Selby et al. 2016). We chose bathymetry, benthic structure, time of day (day 07:00−18:00 h, night 19:00−06:00 h), and year as our covariates. Benthic structure was categorized based on rugosity and structure type as either high-rugosity relief, low-rugosity sand, or low-rugosity hardbottom. High-rugosity relief was characterized by dense, highly rugose physical benthic structure such as reef and rock. Low-rugosity sand comprised homogenous sand interspersed with sea grass pastures, and lowrugosity hardbottom contained low-lying rock interspersed with occasional relief and soft corals. To ac count for location error associated with COAs, we re sampled the benthic structure and bathymetry rasters to a coarser 200 × 200 m resolution. We recalculated the benthic structure raster using the mode of adjacent cell values and used the nearest neighbor technique to recalculate the bathymetry. Year was in cluded as a covariate to account for array expansion. We randomly generated locations restricted to the available habitat where the number of random points equaled the number of COA locations for each individual. This established a null probability of selection of 0.5 for detecting active selection of a habitat type in comparison to random selection.
We fit a GLMM with a logit link function and a binomial error distribution where our response variable was observed (individual location data = 1) or background (randomized point within available area = 0) using the package 'lme4' in RStudio with R version 3.4.3 (Bates et al. 2015, R Core Team 2017. Individual turtle identity was included as the random effect to account for variability among individuals and unequal sampling. We scaled and centered our continuous variables bathymetry and year to help achieve model convergence. Model selection was done using Akaike's information criterion corrected for finite samples sizes (AICc; Burnham & Anderson 2002). We visualized hawksbill habitat selection throughout BIRNM by using coefficients from the best-supported model to predict relative habitat selection at the population level using the combined raster layers and holding the other covariates constant.

Residency
Between March 2012 and March 2017, we captured and tagged 29 hawksbills with an acoustic transmitter (

RSF
We had sufficient detection data (>100 COAs) for 24 tagged individuals for RSF analysis. We compared 4 models using AICc, including a null model with only the random effect of individuals included ( Table 2). The best-supported model had time of day and year as interactive covariates with quadratic 6 Fig. 2. Residency and capture history for acoustically tagged hawksbill turtles. (d) individual was detected within the array on that day. (d) date on which the initial acoustic tag was attached. (s) date on which an individual was captured prior to initial tagging. (×) dates on which an individual was re-captured and a previously attached tag was no longer present, and no new tag was applied. (m) individual previously tagged was re-captured and the tag was still present and assumed to be functioning. (h ×) previously tagged individual was recaptured and the existing tag was removed and replaced with another to continue collecting data. (j) a previously tagged individual whose tag fell off prior to that capture date and was replaced bathymetry and benthic structure. This model had 17 parameters and a model weight of 1. The DAICc was 6129.8 between our top model and a slightly reduced model which used only time of day as an interactive covariate, reducing the parameters to 13. Finally, our most reduced model, which discarded interactive covariates, had only 9 parameters and a DAICc of 8076.9.
Using the estimated coefficients from our top model, we examined relative selection in each of the benthic habitats during the day and night as it varied by depth while keeping year constant at 2015. Relative selection was greater for high-rugosity relief and low-rugosity sand during the day, but selection for those habitats decreased at depths ³20 m (Fig. 3A,B). At night this trend accelerated as depth increased, but there was greater relative selection for lowrugosity sand habitat during the night as opposed to the day (Fig. 3). We found an avoidance of lowrugosity hardbottom structure during both the day and night at any depth (Fig. 3C). Holding bathymetry constant at 5 m and year at 2015 highlighted the relative selection described for each benthic structure during the day and night (Fig. 4).
The best-supported model was used to predict juvenile hawksbill relative selection across the bathy metry and benthic structure rasters available within BIRNM. Holding year constant at 2015, we predicted relative selection across the monument during the day and night (Fig. 5A,B). The areas adjacent to Buck Island to the south and southwest (which include the shallow south forereef and sand/seagrass flats) had the highest relative selection by tagged individuals. Relative selection decreased in the deeper depths farther from Buck Island, especially to the north where low-rugosity hardbottom habitat is predominant. This pattern was more pronounced at night, with tagged individuals selecting those shallow high-rugosity relief and low-rugosity sand habitats close to Buck Island (Fig. 5B). The overall area of selected habitat also appears to decrease noticeably at night.

DISCUSSION
For many highly mobile marine species such as hawksbills, protection from broad scale anthropogenic threats via MPAs can be limited, especially in smaller or less well-regulated MPAs. In contrast, the observed high site fidelity during the post-oceanic juvenile phase presents an opportunity for effective protection using place-based conservation methods and the ability to track changes in protected environments by achieving a quantitative baseline of key species' behavior. Our results help quantify the longterm habitat use by juvenile hawksbills within a coastal Caribbean protected area. Understanding how these individuals use BIRNM is critical to evaluating the efficacy of current management practices, but also provides information for comparison with future cohorts to understand potential shifts in be havior as a result of anthropogenic stressors on the environment there. The RSF analysis provides in sight into those habitats/areas of the monument currently being selected for by resident juvenile hawksbills. We expect that this result would generalize well to other coastal areas throughout the Caribbean where juvenile hawksbills are found, and we predict that increased protective efforts (i.e. habitat restoration and law enforcement) in those shallow water, high-rugosity reef and sandy bottom habitats could improve the survival rate of juvenile hawksbill cohorts in the region. Our results improve our knowledge of thirdorder (i.e. patch-level; Johnson 1980) habitat selection by juvenile hawksbills. More information is still needed at the scale of the specific re sources (i.e. fourth-order habitat selection) within each habitat class to better frame the conservation goals within an area where juvenile hawksbills are found. Finally, as the use of passive acoustic technology continues to grow, especially in studies of marine turtles, our results provide a substantive base for controlled com- Previous research between 1994 and 1998 at BIRNM using CMR methods determined that residency ranged from 59 to 1396 d (mean ± SD = 621 ± 402 d) for 30 hawksbills (including adults) with sufficient recapture data (Hart et al. 2013). These results were consistent with a study by Blumenthal et al. (2009), another relatively long-term CMR investigation conducted in the Cayman Islands. Using this approach, residency is assumed between capture events. Moreover, the in-water sampling at BIRNM was performed within a limited area primarily to the east/ northeast fringe reef surrounding Buck Island. Passive acoustic technology, in conjunction with inwater sampling, provided continuous long-term data on the residency of tagged individuals within BIRNM between in-water sampling events. Residency profiles from tagged individuals reveal frequent use of areas within the array by at least 10 individuals that, despite intense biannual in-water sampling, were never subsequently recaptured (Fig. 2). Compared to Chevis et al. (2017), who similarly used passive acoustic technology to track 22 juvenile hawksbills at Lighthouse Reef Atoll in Belize, we found a wider range of days detected within the array (1−1952 d [our study] vs. 10−1414 d) and a higher mean number of days between first and last detection (693 ± 559 d [our study] vs. 570 ± 484 d). Our results help to further quantify the notion of site fidelity as it pertains to juvenile hawksbills, showing consistent and longterm use of areas within the array over multiple years by some individuals, but also highlighting the plasticity in habitat use among individuals.
The RSF analysis revealed interesting insights into relative habitat selection by juvenile sea turtles at BIRNM. Previous research has shown that hawksbills primarily associate with coral reefs and other complex hardbottom structures (Cuevas et al. 2007, Blumen thal et al. 2009a, Gaos et al. 2012, Hart et al. 2013. These highly rugose habitats provide an abundance of prey items (sponges and macro algae) as well as protection from large predators (Musick & Limpus 1997, Blumenthal et al. 2009b. Results from our RSF analysis supported these findings, showing greater relative selection for high-rugosity benthic structure during both the day and night at depths < 20 m. These habitats are clearly important for juvenile hawksbills. Additionally, our results showed selection for low-rugosity sand benthic structure during the day with even greater selection at night at depths <10 m. This result should be interpreted with caution, as nighttime resting microhabitat is most likely at a finer scale than the 200 × 200 m cell size used to distinguish habitat classes. This scale for habitat classification was chosen to match the scale of the data collected within each receiver's effective detection range, and matching the scale of the location data to the scale of the covariates is critical in estimating resource selection (Boyce 2006). In some areas, such as along the south forereef, where receivers fringe both high rugosity reef and low-lying sea grass pastures/sand, the raster cell size may have been too coarse to capture the microhabitat detail. Juvenile hawksbills have been observed using fringe habitats, similar in composition to those encapsulated within our low-rugosity sand and high-rugosity reef categories, to forage for food and rest in the nearby safety provided by the rugose reef structure (Bjorndal & Bolten 2010, Gorham et al. 2014). Regardless of the particular microhabitat, hawksbills in this study did select these fringe areas at night (Fig. 5B). Increased array coverage in these fringe habitats at BIRNM could help further elucidate microhabitat selection by juvenile hawksbills.
Results from the RSF analysis also showed that the predicted overall area of high relative selection decreased substantially at night. Several studies have suggested that hawksbills use shallow rugose hardbottom structure for resting at night by wedging themselves beneath structure to assist in combating buoyancy and increase dive time (Blumenthal et al. 2009b, Chevis et al. 2017, Wood et al. 2017). Furthermore, Wood et. al. (2017) observed sub-adult individuals using the same 'familiar refuges' at night and a 'nearly 50% reduction in the area occupied' based on kernel density from GPS locations (Wood et al. 2017). Our results suggest a similar behavior by juveniles at BIRNM, and this finding could have potential management implications. First, this resting habitat is presumably crucial to juvenile hawksbills, and without it, the suitability of an area as an MPA for these turtles would be diminished. Second, illegal direct take (i.e. poaching) is a substantial threat even in MPAs, and this highly predictable microhabitat site fidelity could be a behavior that poachers could exploit. Third, identifying these diel patterns informs the study design for ongoing and future surveys of this cryptic species, particularly when it comes to timing surveys when hawksbills are most detectable.
Passive acoustic telemetry, like any technology, has its limitations which are further emphasized when deploying receivers in an open system where entry and exit from the array are not limited by geographic barriers (i.e. rivers, inlets, or bays). Ideally for the RSF analysis, location data and available habitat would be associated with a single receiver or a location determined by triangulation using a VPS system (Freitas et al. 2016). Despite being a rather coarse and simplified method for determining individual location data, we feel that the use of COAs and associated location errors are not drastically different from other more classic telemetry methods, especially at the resolution used to define our covariates. Beyond estimating individuals' locations, inherent to any analysis using passive acoustic technology is an effective array design. The array at BIRNM was a collaboration of numerous researchers investigating many different species, from conch to sharks. Thus, the array grew with areas of extensive and sometimes overlapping coverage and other areas with little to no coverage, depending on study-specific objectives. The RSF framework used with passive acoustic telemetry highlights the notion that re ceivers with no detection data can be as informative as receivers with tens of thousands of detections. Nevertheless, the array coverage for this study was sufficient to provide important insight into the ecology of juvenile hawksbills at BIRNM.
The collaborative effort at BIRNM regarding the passive acoustic array highlights one of the more distinct advantages and disadvantages of the technology. Typically, a study on a single species would not have the capability to finance and maintain a passive array of this magnitude for so long. Even if money and time were no issue, GPS tags would provide location data for hawksbills free from fixed locations at a fraction of the cost. However, we consider the necessary collaborative nature of the array at BIRNM a distinct advantage, as it has not only helped to foster species-specific research for multiple key species found in BIRNM, but it has also created the opportunity to look at important interspecific relationshipsresearch which is currently on-going. Collaborative arrays of this nature are arising around the world, as passive acoustic research shifts from isolated studies to regional networks of arrays (Hussey et al. 2015, Ellis et al. 2019. MPA managers must inevitably determine how to optimally allocate limited funding and effort to protect entire marine ecosystems. Quantifying key species' use of MPA resources provides MPA managers with information in the short-term for analyzing the efficacy of current management practices but also provides a baseline for detecting future changes within the area of interest. Even without high-resolution habitat information, acoustic arrays can yield important information about movement patterns, residency times, and broad scale spatial trends. Collaborative arrays such as the one at BIRNM can provide timely, critical data (absent from other sampling methods) on such species to help inform those decision-making processes.