Habitat models and assessment of habitat partitioning for Kemp’s ridley and loggerhead marine turtles foraging in Chesapeake Bay (USA)

: Understanding the spatial distribution of a species is required to enact effective conservation. Complications to effective conservation can arise when the distributions of multiple target species are non-overlapping. Conservation efforts meant to protect one species may shift threats into the distribution of another species. Two species of marine turtle, loggerhead Caretta caretta and Kemp’s ridley Lepidochelys kempii , are common seasonally in Chesapeake Bay, a large estuary on the US east coast. Both species are protected under the US Endangered Species Act and face spatially complex threats in the region. We created habitat suitability models for these 2 species to inform conservation efforts in the region and explore the extent of overlap be tween their distributions. Argos satellite tags were deployed on 24 Kemp’s ridley and 10 loggerhead turtles to record animal locations within the Bay. Boosted regression tree models were created for each species using presence-only animal locations, predicting suitable habitat within the Bay. Habitat for Kemp’s ridley turtles was predicted in shallow, coastal areas of the southern Bay as well as in brackish areas of rivers. Loggerhead turtle habitat was predicted to extend farther north than Kemp’s ridley habitat and was generally found in deeper areas of the middle Bay. There is some evidence that these 2 species are partitioning habitat. Any conservation measures adopted to conserve marine turtles in the Chesapeake Bay should consider the habitat of both species holistically to avoid shifting impacts from one species to another.


INTRODUCTION
Conservation measures for widely ranging marine fauna often focus on protecting critical life stages or activities such as breeding, foraging, or migrating (Corrigan et al. 2014, Lascelles et al. 2014, Calambokidis et al. 2015. The effectiveness of these conservation measures is dependent on understanding the distribution of a species during these activities, thus al lowing for targeted intervention at appropriate spatial and temporal scales (Norse 2010, Davies et al. 2012, Guisan et al. 2013. Implementing conservation measures for a species can be complicated by the presence of other protected species in the same region with non-congruent distributions. Implementation of conservation measures for one species can shift deleterious effects from one species to another (Salas & Gaertner 2004, Barrows et al. 2005, Abbott & Haynie 2012. For example, by imposing fishing closures in one area to avoid bycatch of a protected species, that fishing effort may shift to where another protected species resides.
Chesapeake Bay (hereafter 'the Bay') is a large, salt-wedge estuary on the mid-Atlantic coast of the USA that is home to economically and strategically important commercial and naval ports, commercial fisheries (Zohar et al. 2008, Richards & Rago 1999, Rick et al. 2016, substantial recreational boating and fishing fleets, and protected marine species (Barco et al. 2018a, Kahn et al. 2019, Aschettino et al. 2020. Overlapping use of the Bay by commercial and naval vessel traffic, fishing gear, private vessels of all sizes, and protected species is complex and not fully understood in cases where information on the habitat and distribution of protected species is limited. More knowledge of the in-water distribution of protected species in the Bay is required to manage this complex situation effectively. The Bay is a seasonally important foraging ground for loggerhead Caretta caretta and Kemp's ridley Lepi do chelys kempii turtles, as was previously documented from the mid-1980s to the early 2000s (Lutcavage & Musick 1985, Keinath et al. 1987, Seney & Musick 2005, Mansfield 2006, Seney & Musick 2007, Mansfield et al. 2009). Loggerheads found in the Bay are part of the Northwest Atlantic Distinct Population Segment and are listed as threatened under the US Endangered Species Act (National Marine Fisheries Service 2011). Kemp's ridleys are a single population and are listed as endangered under the US Endangered Species Act. Both of these protected species face numerous threats in the Bay, including, but not limited to, vessel strikes from commercial, naval, and recreational vessels; dredging activities for shipping channels and beach re-nourishment; commercial and recreational fishery bycatch; climate change; and naval training and testing activities (Barco et al. 2015), but the extent of their distribution is not fully understood.
More recent satellite-telemetry studies found that both species spend most of their time engaged in area-restricted search behavior, generally interpreted to be foraging (Barco et al. 2017(Barco et al. , 2018b, a critical life function. Home range analyses (Barco et al. 2018b) identified presumed foraging habitat for Kemp's ridley turtles in the southwestern corner of the Bay, the James and York Rivers, and several other nearshore locations, although there was extensive variation in home range size and location among individuals (Di-Matteo et al. 2021). In contrast, presumed loggerhead turtle foraging habitat was found primarily in the center of the Bay, into the waters of southern Maryland, with some additional habitat closer to shore (Barco et al. 2017). Like Kemp's ridleys, there was extensive variation in the size and location of individual loggerhead home ranges. These animals are predominantly juveniles, but adults have been documented as well.
A sensitivity analysis using a Kemp's ridley satellite telemetry dataset and tag simulation found that too few tags had been deployed to identify all suitable habitat likely to be present in the Bay (DiMatteo et al. 2021). Habitat modeling was identified by that study as an approach to identify additional potential habitat for both Kemp's ridleys and loggerheads in the Bay.
Habitat models for marine species, in this case habitat suitability models using presence/absence data, can: (1) describe complex relationships among species and environmental covariates, (2) be extra polated in space and time (with caution), and (3) provide insights into species' distributions (Robinson et al. 2017). Many frameworks for habitat suitability models exist, including, but not limited to, generalized additive models (Hastie & Tibshirani 1987, Guisan et al. 2002, maximum entropy (Elith et al. 2011), and boosted regression trees (BRTs, Elith et al. 2008). Here we used the machine-learning technique of BRTs to create habitat suitability models for both species.
BRTs allow for the fitting of complex environmental relationships, explicitly explore covariate interactions, can include factors as covariates, and are robust to outliers in the dataset as well as missing covariate values (Elith et al. 2008). BRTs are a combination of regression trees, a type of decision tree model, and the boosting technique, which produces many simple tree models and then combines them to maximize predictive utility (Elith et al. 2008). These features make BRTs a good option for this tagging dataset and location, as sea turtles may be responding to complex environmental cues that create ephemeral conditions which aggregate prey or create conditions advantageous to their physiology (Schofield et al. 2009, Howell et al. 2015. The Bay is a complex estuarine environment with high variability in environmental conditions (Preston 2004, Werdell et al. 2009), increasing the potential for animals to react to complex environmental cues.
We also explored the possibility that these 2 sea turtle species were partitioning their habitat (e.g. using different resources) by comparing the outputs of the 2 habitat suitability models. Habitat partitioning may complicate conservation efforts if animals are found in different areas of the Bay. Habitat partitioning allows species to avoid or reduce competition by occupying different ecological niches (Whittaker 1967, May & MacArthur 1972. While both loggerhead and Kemp's ridley turtles feed on benthic inver-tebrates (Byles 1988, Keinath et al. 1987, Mansfield 2006, as a population, Kemp's ridley turtles from the region have shown a preference for blue crabs Callinectes sapidus, which are found predominantly in shallow areas of the Bay (van Engel 1958, Millikin & Williams 1984, Seney & Musick 2005, 2007, Barco et al. 2015. Horseshoe crabs Limulus polyphemus and whelks (Busycon spp.) are the predominant components of the diet of loggerheads in the Bay and can be found in deeper waters (Seney & Musick 2007, Barco et al. 2015.
The habitat suitability models presented below can provide managers with information to conserve these species in a complex estuarine environment with multiple stakeholders and uses. The extent of overlap between the habitat suitability of the species will yield insights into whether habitat partitioning is occurring, which will impact the complexity of conservation measures needed to protect both species. Additionally, the habitat models can be used as an input for designating critical foraging habitat for these species, which is not currently designated.

Turtle-tagging and satellite-telemetry data
Data were analyzed from Argos satellite tags deployed on 24 Kemp's ridley and 10 loggerhead turtles (Table 1). Details on capture and tagging methods can be found in Barco et al. (2015Barco et al. ( , 2017Barco et al. ( , 2018b brief, tagged animals were a mix of wild-caught and rehabilitated animals. All animals, regardless of origin, underwent a health assessment prior to tagging and release. Animals were tagged in accordance with United States Fish and Wildlife Service requirements for minimizing drag and proportion of tag weight to animal weight (Eckert et al. 1999, Jones et al. 2011, USFWS 2016, with a variety of Argos satellite-linked platform terminal transmitters. Multiple tag models were used based on budgeting constraints in different years of the program and the need to have a variety of Argos tag sizes to attach to animals with different carapace lengths and weights. Some existing deployments were not included in the analysis: animals whose deployments were entirely outside the Bay; animals that spent fewer than 5 d in the Bay; animals deployed prior to 2014, which was the earliest year for which the selected environmental covariates were available; and 1 Kemp's ridley turtle whose entire deployment was in a small inlet not covered by the environmental covariates. Information on excluded animals was not included in Table 1. Argos satellite deployments for the remaining turtles were run through the Douglas filter (Douglas et al. 2012) to remove unrealistic locations by using settings recommended for hardshell turtles by the Turtle Expert Working Group (2009). Argos location errors can be up to 5 km, depending on the quality of the satellite fix, and Argos locations are not true animal locations. The Turtle Expert Working Group Douglas filter settings included parameters for the maximum redundant distance filter and the distance angle rate hybrid filter algorithms, which were used to account for unrealistic animal speeds and turning angles. Additional filtering included removing locations on land and removing locations 24 h post-deployment to account for behavioral changes associated with animal handling and tagging.

State−space modeling
Locations retained after filtering were processed using a hierarchical state−space model (hSSM), segregated by species, to create estimated locations at regular time intervals from the irregular locations reported by the Argos tags using the R package 'bsam' in R version 4.0.2 (Jonsen et al. 2005, Jonsen 2016, R Core Team 2019). These estimated locations reduced spatial autocorrelation from reported locations caused by animal behavior and satellite coverage. An hSSM was selected because movement parameters are estimated jointly for all individuals, along with individual effect parameters, allowing shorter deployments to benefit from the information in longer deployments. This assumes that animals' behavior is broadly similar, which we feel is reasonable given that animals were the same species, similar age class based on carapace length (juvenile), engaged in similar behavior (foraging), and were in the same region.
We did not include a behavioral component in the hSSM to distinguish between traveling and area-restricted search (ARS) behavior. Previous state−space modeling found that both species spent the majority of their time in the Bay engaged in ARS movements, assumed to be foraging (Barco et al. 2017(Barco et al. , 2018b. This is consistent with the known ecology of the species in this area (Lutcavage & Musick 1985, Keinath et al. 1987, Seney & Musick 2005, Mansfield 2006, Mansfield et al. 2009). Given that relatively few locations from previous studies in the Bay represented traveling behavior, we assume areas identified by the habitat model primarily represent foraging habitat.
A 6 h time step was chosen as the interval between locations in the hSSM, as this was the finest-scale time step that could be fit to the Kemp's ridley data. The loggerheads had more locations reported per day than the Kemp's ridleys, which allowed for hSSMs as fine-scaled as 3 h in exploratory models, but we opted to keep the intervals equal so that each deployment would have the same number of locations reported per day, regardless of species. hSSM diagnostics were examined to ensure that Monte Carlo Markov chains were mixing, that parameters estimates were converging, and that autocorrelation be tween chains was acceptably low. After state−space modeling, estimated locations from the hSSM were reviewed visually to ensure that the resulting deployments were reasonable. Points on land were removed using the global self-consistent, hierarchical, high-resolution geography database full-resolution shoreline (Wessel & Smith 1996). Fig. 1 displays the locations retained for the analysis and the study area.

Environmental covariates and absence data
Candidate covariates for the habitat models included a mix of static and dynamic physical covariates, habitat maps, and temporal factors (Table 2). Biological covariates such as chlorophyll a were not available at appropriate geographic and temporal scales and were not considered. Biological covariates were only available as remotely sensed or ocean model products and were at too coarse a spatial reso-lution, generally 5 km 2 or greater, to reflect fine-scale habitat use in the topographically complex Bay.
Because the Bay is a highly dynamic environment, covariates were chosen to make predictions on the finest temporal and spatial resolution possible, to reflect potential fine-scale features. The Chesapeake Bay Operational Forecast System (CBOFS) provides daily readings of temperature (°C) and salinity (parts per thousand, ppt) from sampling stations within the Bay and some of the major adjoining river systems. Both temperature and salinity are potentially important environmental conditions for turtles and their prey.
Temperature and salinity rasters were derived from CBOFS in situ sensors (n = 157) by taking the noon reading from each sensor and using the diffusion interpolation tool in ArcGIS 10.7 (Esri 2019) to interpolate values to match the extent of high-resolution bathymetry data for the region. The resulting daily rasters had a resolution of 500 m.
CBOFS data were not available prior to July 2014, so all ARGOS locations prior to 2014 were removed from the analysis. Rasters were processed for the months of May−November (with the exception of May and June 2014), the months when telemetered turtles deployed for this study were present in the Bay.
Two temporal covariates, day of year and year, were derived from the times of estimated locations, and were included to account for intra-and interannual variability in habitat not accounted for by other covariates. Latitude was also included as a candidate covariate to account for north−south variability in habitat not captured by other covariates.
Static maps of benthic habitat (National Marine Fisheries Service 2020) and submerged aquatic vegetation (Lefcheck et al. 2018) were also assessed as potential predictors based on the assumption that they may reflect important habitat for prey species. Ultimately, these habitat maps were dropped from consideration as there were many missing values compared to the extent of other available covariates and the submerged aquatic vegetation data also had missing years (the benthic habitat map did not vary temporally). Even though BRTs are generally robust to missing covariate values, here the missing data were extensive both spatially and temporally.

BRTs
Interpolated animal locations from the hSSMs for loggerheads (n = 2247) and Kemp's ridleys (n = 2738) were used as presence data for their respective BRT models. The BRTs also require absence data to fit the model. Because satellite-tag-derived locations were used for presence data, confirmed absences did not exist for this study. Pseudo-absences were used instead and were generated in 1 of 2 ways: (1) Random absences equal to the number of presence samples were created by selecting a random location on a random day from within the study area and time period. Using the same number of presence and pseudo-absence locations can have the best predictive performance in machine-learning techniques such as BRTs (Barbet-Massin et al. 2012). Random locations were drawn from the centroids of raster cells of the bathymetry covariate. Using randomly generated pseudo-absences assumes that animals are not distributed randomly in the study area and that the selected covariates will be able to distinguish between the random absences and the true habitat of the target species (Hirzel et al. 2002). The same method was used to generate 2 sets of random absences, one for Kemp's ridley and one for loggerhead datasets (2) Target-group absences, where confirmed presences of a related but distinct species are used as absences (Phillips et al. 2009). Target-group absences can outperform randomly generated pseudo-absences in some cases (Cerasoli et al. 2017). Loggerhead turtle locations were used for absences in the Kemp's ridley model, and vice versa. Target-group absences were unweighted, as the number of locations of Kemp's ridley and loggerhead turtles were similar. The targetgroup absence method assumed that the 2 species were utilizing distinct habitats within the study area. This may be the case in the Bay (Barco et al. 2018b).
The modeling process was started by fitting exploratory models with random absences that included all possible covariates and examining covariate interactions and contributions to the model. For covariate pairs that were highly correlated, the covariate with a lower contribution to the model was removed to simplify the models. In this exploratory phase, tree complexity, learning rate, and bag fraction, the primary parameters for adjusting how BRTs are fit, were adjusted manually to improve model performance.
Model performance was assessed by examining re sidual deviance, coefficient of variation (CV), and the area under the receiver operating characteristic curve (AUC). Respectively, these assess the explana-tory power of the model, the associated uncertainty, and true positive rate compared to the false positive rate at various thresholds.
Sample predictions for random days, the finest temporal scale of the available covariates, were made throughout the exploratory modeling phase to assess whether predictions contained obviously spurious artifacts. If a covariate was causing clearly unrealistic predictions, it was dropped from subsequent models.
After the exploratory phase, the retained covariates were used to create 4 final models for each combination of species and absence generation method -Kemp's ridley with random absences, Kemp's ridley with target-group absences, loggerhead with random absences, and loggerhead with target-group absences -in a more systematic fashion. Final models were selected from a set of candidate models where tree complexity, learning rate, and bag fraction were systematically changed be tween commonly used values for each (Elith et al. 2008): 1, 2, and 3 for tree complexity; 0.1, 0.05, 0.01, 0.005, 0.001, and 0.0005 for learning rate; and bag fraction values of 0.1−0.9 in increments of 0.1. This resulted in a total of 162 candidate models for each final model, based on all possible combinations of parameter values.
The final models were selected from the candidate models by comparing residual deviance, CV, and the AUC among candidate models. If one candidate model did not score the best on all 3 metrics, professional judgment was used to assess the relative quality of top-scoring models and make the final selection. The selected final models of the 2 absence-generation methods were compared and the better of the 2 was carried forward. This selection was made by examining both model diagnostics and a qualitative assessment of the resulting predictions.

Habitat prediction and analysis
The study area (e.g. the area covered by predictions) was the Bay, including Virginia and Maryland waters, and connected riverine systems as defined by the extent of the available bathymetry, salinity, and temperature covariates. This included areas farther north and farther upriver than individuals of either species occurred based on the available satellite tele metry locations. These areas were included to see if suitable foraging habitat existed that telemetered turtles had not utilized, or if environmental conditions were different from areas frequented by turtles in our dataset.
Habitat suitability values from the BRT models range from 0 (poor habitat) to 1 (excellent habitat).
These habitat suitability scores were subsequently classified into habitat versus non-habitat by selecting a cutoff value that maximizes the ratio of true positives to true negatives predicted by the model (Bradley 1997).
To assess habitat quality over the entire period of the study, which covered 5 seasons from 2014 to 2018, averaging habitat suitability scores over many days did not seem reasonable. This would have the effect of 'washing out' suitable habitat as conditions changed day to day.
Instead, the number of 'habitat days' present at locations in the Bay were assessed at various temporal scales. Habitat days were defined as the number of days suitable habitat was found in each predicted location within the Bay. For each final model, habitat suitability was predicted for each day within the study period. Daily habitat suitability surfaces were reclassified into habitat and non-habitat using a cutoff value.
The cutoff value was determined by randomly splitting the modeling data sets (presences and random absences) into training (70%) and testing (30%) datasets. The model was then fit with the training data and used to predict the testing data, which allowed the cutoff value (Bradley 1997) to be determined. This process was repeated 10 times. The cutoff value from the 10 replicates was averaged and subsequently used to reclassify the daily prediction rasters. This process was performed independently for each of the 4 final models.
This process yielded a set of rasters representing daily foraging habitat within the Bay, classified as 1 = habitat, or 0 = non-habitat. These rasters were summed by month and across all days in the study to assess the number of days suitable foraging habitat could be found in a location. Higher numbers of 'habitat days' were an indicator of higher-quality foraging habitat, with suitable conditions found there more often than elsewhere.
The results of summing the rasters across all days of the study were reclassified into 'important habitat' for the selected final model of each species. Important habitat was defined as the quartile of locations (raster cells) with the highest numbers of habitat days. This provided a map of the most suitable habitat areas within the Bay, for each species, over the entire study period. These maps were then used to assess habitat partitioning between the 2 species.
Habitat partitioning was examined using 3 different metrics: Schoener's D (Schoener 1968), Hellinger'sbased I (Van der Vaart 1998, Warren et al. 2008), and Syrjala's test (Syrjala 1996). Schoener's D calculates the range of a species based on probability distributions of abundance over a set of locations and calculates niche overlap based upon species abundance in those locations. Hellinger's I is based on probability distributions without the assumptions of Schoener's D (Warren et al. 2008, Hosseinian Yousefkhani et al. 2016). Syrjala's test assesses whether 2 distributions are equivalent, invariant of abundance.
We assumed habitat days corresponded proportionally to species occupancy and converted the overall habitat-day rasters to proportional occupancy by summing all cells (total habitat days across all cells) and then dividing the values of individual cells by that total sum. It is unlikely that habitat days correspond exactly to occupancy, but we could not test this assumption without extensively surveying the Bay.

State−space modeling
hSSMs for both species performed acceptably, converging with 80 000 posterior samples for the Kemp's ridley model and 50 000 samples for the loggerhead model. Both models used 10 000 samples as an adaptation phase, and a span parameter of 0.2. hSSM diagnostics for both models indicated that Monte Carlo Markov chains were mixing, that parameter estimates converged based on the Gelman-Ruben shrink factor (Gelman & Rubin 1992, Brooks & Gelman 1998, and that autocorrelation between chains was acceptably low. After additional manual filtering, 2738 Kemp's rid ley and 2247 loggerhead turtle locations remained as presence locations for the BRT models. Modeled locations closely followed Argos deployments, with outlier locations and points on land removed. Locations for both species generally matched previous descriptions of the species' distributions, though Kemp's ridley turtles were found farther upstream in several rivers than previously reported. The northward distribution of the locations extended past the Patuxent River for loggerhead turtles and to the Nanticoke River for Kemp's ridley turtles. This matched the previous northward extent of marine turtles recorded by aerial surveys in the Bay (Barco et al. 2018a).

BRT models
Exploratory BRT models indicated that day of year was highly correlated with temperature and that lat-itude was highly correlated with salinity. Additionally, models that included latitude only predicted suitable habitat at latitudes where animals were located. Because of these factors, the day of year and latitude covariates were dropped from consideration in the final models. While BRTs can include correlated covariates, we preferred to limit the complexity of the models and to preferentially include covariates directly related to habitat conditions where possible.
Changes in predictions that included year as a covariate were more closely related to the number of animal locations in that year than underlying changes in environmental conditions. Because there were not similar numbers of animals tracked or locations in each year, or a way to reasonably standardize effort between years, year was removed as a candidate covariate.
This left temperature, salinity, and bathymetry as available covariates. The best model for all 4 combinations of species and pseudo-absence generation method included all 3 remaining covariates. Table 3 shows the values for BRT and model assessment parameters for each selected model. The selected final models all had at least 1000 trees, as recommended by Elith et al. (2008).
For Kemp's ridleys, the model with randomly generated pseudo-absences was selected as the best final model. Even though the target-group absence model performed better when examining the model assessment values, on qualitative review of the predictions, the target-group absence model predicted most of the suitable habitat north of the Choptank River in Maryland, much farther north than Kemp's ridley turtles have been detected previously. Little suitable habitat was predicted in the southern Bay. This did not accurately reflect either the known ecology of the species, the distribution based on stranding data, or detections from the acoustic array present in the Bay (Barco et al. 2018b).
Based on the functional plots of covariates, Kemp's ridley turtles showed a preference for depths shallower than 15 m, temperatures ranging from 17 to 28°C, and salinities of 15−28 ppt (Fig. 2a). For reference, the average depth of the Bay is 6.4 m. Each of the 3 covariates had relatively equal importance in the model: bathymetry = 37.3%, temperature = 31.7%, and salinity = 31%.
For loggerheads, the model with target-group pseudo-absences was selected as the final model. The target-group model performed better than the random model when comparing model assessment values (Table 3). In the qualitative assessment, both models predicted similar extents of suitable habitat, but the target-group model predictions were more compact and consistent across time periods.
Based on the functional plots of covariates, loggerheads appeared to avoid depths less than 10 m and salinities less than 20 ppt. Loggerheads appeared to prefer temperatures warmer than 25°C, and salinities greater than 15 ppt (Fig. 2b). This contrasts with Kemp's ridley turtles, which occurred in shallower waters with a more defined temperature range. For loggerheads, bathymetry was the most important covariate (43.6%), followed by salinity (29.5%) and temperature (26.9%).

Predictions and habitat partitioning
For the Kemp's ridley model, a cutoff value of 0.48 was used to partition daily habitat suitability predictions into habitat versus non-habitat. Values ≤ 0.48 were classified as non-habitat, and values > 0.48 were classified as habitat, which maximized the ratio of true positive to true negative predictions in tests where Kemp's ridley locations and randomly generated absences were randomly split into training and testing data sets.  Table 3. Best models for each combination of species and absence generation method, including boosted regression tree (BRT) parameters, residual deviance, coefficient of variation (CV), and the area under the receiver operating characteristic curve (AUC). Values with an asterisk indicate that they were the best from within each set of models Aggregating the daily habitat predictions by month, we saw the most suitable habitat for Kemp's ridleys, both total number of days and total area, in June and October and the least in November (Fig. 3). Suitable habitat occurred first in the southern Bay, early in the season, extending northward into Maryland as the summer progressed, and retreating into the southern Bay in November. These latitudinal shifts in prediction were driven largely by temperature, which peaks in mid-to late summer. The dip in predicted habitat suitability in July and August was likely caused by the model predicting Kemp's ridley turtles do not prefer habitat warmer than 28°C.
Aggregating all daily habitat predictions into a single surface, the areas with the highest suitability had over 700 suitable habitat days over the study period (5 yr), or >120 d yr −1 . The areas with the highest number of habitat days were in nearshore Virginia waters characterized by shallow depths, moderately high temperature, and higher salinity Percentages are the relative importance of the covariate to the model than the northern reaches of the Bay. Some suitable habitat was identified in Maryland waters, but it was, in general, predicted to be less suitable than the lower and middle portions of the Bay in Virginia. Almost no important habitat (top quartile of habitat days by grid cell), was identified in Maryland waters. The total predicted important habitat covered 2808 km 2 (see Fig. 5a) and was primarily located in nearshore areas of the southern Bay and the James and York Rivers. For the loggerhead model, a cutoff value of 0.46 was used to partition daily habitat suitability predictions into habitat versus non-habitat. Aggregating the daily habitat predictions by month, we saw the most suitable habitat, both total number of days and total area, in August and the least in June (Fig. 4). In general, suitable habitat was in the southern Bay in spring and early summer, extending northward into Maryland as the summer progressed, and retreating into deeper areas in November. In July−September, habitat was predicted farther north than loggerheads are generally seen. This is likely driven by the importance of depth to the model and the occurrence of loggerhead turtles in more saline conditions compared to Kemp's ridleys (Fig. 2). In the southern Bay, shipping channels, which are dredged regularly, were clearly highlighted 100 Fig. 3. Number of suitable habitat days by month for Kemp's ridley turtles using the random pseudo-absence model as habitat in all months ex cept July and August (Fig. 4). The shipping channels are not explicitly marked in Fig. 4 but can be seen as bright red and orange linear features in the southern Bay in most months.
Aggregating all loggerhead daily habitat predictions into a single surface, the areas with the highest suitability had over 1008 suitable habitat days over the study period, or > 200 d yr −1 . The areas with the highest number of habitat days were in the southern, central Bay with the deepest depths, moderately high temperature, and relatively higher salinity. Some suitable habitat was identified in Maryland waters, but those areas were predicted to be less suit-able than lower portions of the Bay. The small amount of important habitat found in Maryland waters was limited to deep areas in the central Bay and the central reaches of a few larger rivers. The total predicted important habitat covered 2775 km 2 (Fig. 5b) and was predominantly in the southern, central Bay in contrast to 2802 km 2 in nearshore southern areas for Kemp's ridley's (Fig. 5a).
Comparing the selected Kemp's ridley and loggerhead turtle models, the Schoener's D-score was 0.52 and Hellinger's distance (I) was 0.82. These scores indicate there is some evidence that these 2 species are inhabiting different areas. Syrjala's test p-values were < 0.0001 for both the Cramer-von Mises and Kolmogorov-Smirnov tests (1000 permutations). Rejection of the null hypothesis indicates that the spatial distribution of the 2 species is significantly different. On the assumption that our habitat model accurately represents the species' spatial distribution in the region, this is strong evidence that the species are partitioning habitat and using different resources. Overlaying animal locations on top of predicted habitat (not shown) indicates the assumption that our models reflect species' distributions appears to be reasonable.

DISCUSSION
This study provides the first habitat suitability models for Kemp's ridley and loggerhead turtles in the Chesapeake Bay. These models are an important resource for managers, as previous tracking and home range studies indicated that it was not feasible to track enough of these animals to identify all potential foraging habitat in the Bay (DiMatteo et al. 2021). Kemp's ridley habitat was identified in rivers, river mouths, and shallow areas of the Bay, whereas loggerhead habitat was predicted in deeper areas of the central and southern Bay, showing evidence of habitat partitioning between these 2 species.
Examining reasons for why the Kemp's ridley target-group absence model predictions performed poorly, we saw that, like the loggerhead targetgroup absence model, bathymetry was by far the most important variable (51.4%). This meant that shallow areas in the north of the Bay were classified as suitable habitat even though animals were not sighted there. This could have been ameliorated by limiting predictions to the southern Bay, but one of the intents of this study was to identify possible habitat areas where animals were not tracked.
It could be argued that the areas identified by the target-group absence model in the north of the Bay meet this objective. However, the identified habitat was extensive and was the artifact of a single covariate relationship, rather than what we would consider to be a realistic depiction of unutilized habitat. This River, which divides the Maryland and Virginia portions of the Bay, and no Kemp's ridley turtles stranded north of latitude 38.8°N (Barco et al. 2015). The functional relationships of the Kemp's ridley randomly generated absence model were more defined, and the 3 variables had similar influence. This may more accurately reflect the preferences of Kemp's ridley turtles, as the randomly generated absences sampled a broader range of the environmental covariates than the target-group absences. Preferred habitat values from the randomly generated absence model were consistent with previous studies of Kemp's ridleys within the Bay (Lutcavage & Musick 1985, Keinath et al. 1987, Seney & Musick 2005, 2007, Mansfield 2006, Mansfield et al. 2009).
The randomly generated absence model still predicted some suitable habitat north of areas where Kemp's ridleys have been generally sighted. It may be that traveling to northern areas of the Bay is not energetically cost effective for animals during the foraging season or that other factors exist that make these areas less suitable which were not elucidated by the covariates used in our model.
Model performance may be limited by the scale and availability of covariates (Scales et al. 2017, Cox et al. 2018. Solar insolation, bottom temperature, density of submerged aquatic vegetation, and other covariates related to sea turtle physiology could refine our predictions (Abecassis et al. 2013) but are not currently available in the Bay at appropriate temporal or geographic resolutions. Likewise, prey availability in the form of crustacean and mollusk distribution and abundance, which is likely a driving factor of sea turtle distribution in the Bay, is difficult to quantify on small temporal and spatial scales (King et al. 2005). These data would be difficult to acquire because of both the logistics involved with surveying and the highly dynamic nature of prey availability.
There are important insights that can be drawn from the covariates that were available. The functional relationships from the Kemp's ridley randomly generated absence model predict that Kemp's ridleys avoid very warm surface temperatures (> 27°C). This may be true, as juvenile Kemp's ridleys are small, may not be able to thermoregulate as well as larger turtles (Still et al. 2005), and may not have cold-water refugia in the shallower depths that they appear to prefer. Bottom temperatures would be a better re flection of the actual temperatures to which turtles were exposed while foraging, but high-resolution bottom temperature data were not available for the spatial and temporal scale of this study. This predicted relationship, whether true or not, is likely the driving factor for lower suitably being predicted for Kemp's ridley turtles in the months of August and September, which are the warmest months in the Bay.
Of concern for loggerhead turtles are the shipping channels being highlighted as suitable/important habitat, based on the predicted avoidance of shallower waters. These channels are dredged regularly, and loggerheads are the primary species identified to be affected by dredging operations in the Bay (Mansfield & Musick 2003, NMFS 2018. If loggerhead turtles do use these areas regularly, they are at in creased risk from dredging operations and ship strike from the shipping and naval vessel traffic using these deepwater channels. The importance of depth to the loggerhead turtle habitat model could be the reason that these areas are highlighted (e.g. an artifact of the model) but the relationship merits further investigation given the potential conservation implications.
Loggerhead turtles may be recorded as the primary species impacted by dredges because they are larger and more likely to be reported. Kemp's ridley turtles using river habitat may also be susceptible to dredging in rivers. The US Department of Defense regularly dredges the York River to ensure access to the naval facilities there. In the future, collecting dive data for both species to characterize 3-dimenstional habitat use within the Bay could yield better insights into depth selection and whether animals are utilizing the entire water column. Currently no published dive data for these species exist within the Bay.
Comparing the selected Kemp's ridley and loggerhead models, we saw evidence that these species may be partitioning their habitat, given the assumption that our habitat suitability models are equivalent to occupancy (which would have to be derived from extensive surveys). We posit that partitioning is driven by the distribution of the preferred prey of each species, i.e. blue crabs for Kemp's ridley turtles, compared with horseshoe crabs and other benthic invertebrates for loggerhead turtles (Barco et al. 2015). In the Bay, blue crabs are found primarily in shallow, vegetated areas. A study of the habitat of these 2 species in the Gulf of Mexico also showed evidence of habitat differentiation (Fujisaki et al. 2020).
This habitat partitioning complicates conservation efforts, as area closures or restrictions targeting one species may shift risks to the other. We recommend that any mitigation or conservation measures for turtles be applied to the entire southern Bay and southern rivers, as these will be the most effective for protection efforts of both species. Restrictions in the southern Bay may shift impacts to the northern areas of the Bay, where turtles do not appear to be as abun-dant and where less suitable habitat exists, though this could impact non-turtle species, which are discussed below. The timing of any restrictions targeting mitigation of impacts to turtles should cover the core months of the foraging season for the 2 most common species of marine turtles (June−September).
Given the economic importance of the region, it is unlikely that conservation efforts that include outright closures will be reasonable. The shipping channels that naval and commercial traffic depend on have minimum depth requirements and are unlikely to be moved and changed. Additionally, commercially and recreationally valuable species such as blue crabs may only be available in certain habitats, resisting efforts to shift impacts elsewhere.
Other protected marine species utilize the Bay, including 2 additional species of marine turtles, and several species of marine mammals and fish. Green turtles Chelonia mydas and leatherback turtles Dermochelys coriacea are less common in the Bay than Kemp's ridleys and loggerheads based on historic surveys (Keinath et al. 1987) and recent stranding data (Barco et al. 2015, Swingle et al. 2018. Little tracking has been performed on green turtles in the region and none on leatherbacks, making assessing these species in the context of loggerhead and Kemp's ridley habitat use impossible at the current time. If adequate tracking of green and leatherback turtles occurs, overlap of these species with Kemp's ridley and loggerhead turtles can be explored, although the diets and likely ecological niches of green and leatherback turtles are significantly different from Kemp's ridley and loggerhead turtles (reviewed in Bjorndal 1997, Jones & Seminoff 2013 and may preclude significant niche overlap.
Atlantic sturgeon Acipenser oxyrinchus utilize multiple river systems adjacent to the Bay for spawning and pass into and out of the mouth of the Bay while migrating (Welsh et al. 2002, Hager et al. 2014. Humpback whales Megaptera novaeangliae, gray seals Halichoerus grypus, and harbor seals Phoca vitulina use the mouth of the Bay and nearshore waters outside the Bay seasonally (Aschettino et al. 2020, Ampela et al. 2021. Bottlenose dolphins Tursiops truncatus can be found in the Bay year-round, their range extending farther north than that of loggerhead turtles (Barco et al. 1999, Richlen et al. 2018. Given these additional conservation targets in the Bay, comprehensive and dynamic marine spatial planning is required to adequately balance the complex spatial and temporal overlaps of protected species, threats, and economic and recreation resources. Dynamic ocean management can lessen the time and space needed for closures (Maxwell et al. 2015) and can incorporate the needs of many different stakeholders and priorities. This is particularly important given the economic importance of fisheries and vessel traffic in the Bay, where dynamic management can accomplish the management goals of minimizing time and space closures while adapting to complex bio logical and ecological processes (Dunn et al. 2016). Cooperation will be required from multiple state agencies, the federal government, conservation organizations, and commercial and recreational interests. An adaptive management framework is in place for the Bay (Chesapeake Bay Program 2020), bringing together many of these parties across all states that border the Bay. Incorporating and informing more dynamic management processes into this framework can improve the overall management of the Bay. These habitat models can inform marine spatial planning and adaptive/dynamic management processes.
Lastly, critical habitat has been designated for the North Atlantic loggerhead turtle distinct population segment (NMFS 2014). However, no foraging areas were designated as part of this determination, due to lack of prey availability data fine scaled enough to identify areas as foraging habitat (without identifying the entire Bay). Critical habitat has not been designated for Kemp's ridley turtles. The models presented here could be used to inform critical habitat designations for these species which take account of the critical life function of foraging.