A detailed understanding of how protected species use their habitats can guide management interventions in areas of high human use. For marine turtles, different food availability and physical habitat characteristics can underpin turtle presence at anthropogenically modified compared to unmodified sites. We develop telemetry-based habitat models with boosted regression trees to identify the environmental characteristics underpinning foraging habitat suitability for green turtles in the Great Barrier Reef region. We fit models to green turtle Fastloc GPS tracks from both modified and unmodified inshore foraging sites and using  pseudo-absences (simulated correlated random walks). We assess model performance by the ability to predict known foraging areas, true skill statistic, explanatory power (percent deviance explained) and predictive skill (AUC) of the models. We then predict potentially suitable foraging areas for green turtles in the Great Barrier Reef region using the model for unmodified habitats. Between 2010 and 2022, the total area of suitable foraging habitat declined by 41.2%, and nearshore habitat suitability retracted. These areas are likely affected by floods, development and increased turbidity. In 2022, 50% of predicted suitable habitat fell within habitat protection zones, and 19.4% in Marine National Park Zones of the Great Barrier Reef Marine Park. A detailed foraging distribution of the species has not previously been compiled at this regional scale. Identifying biophysical drivers of habitat suitability can inform identification of possible foraging habitat in less data rich regions in Australia and overseas. Evaluating changes over time in habitat distribution provides insights into the degree to which broad-scale environmental changes and anthropogenic activities influence the condition and function of habitats, even within protected area boundaries.

Code supporting this dataset is available at https://github.com/egwebster/SSM-SDM-public

Turtles were tracked during foraging at Port Curtis, Shoalwater Bay and post-nesting from Raine Island. Turtles were tracked with Fastloc-GPS (mean error 40m) and/or ARGOS tags (mean error ranging <1km to >10km) opportunistically between 2010 and 2021, resulting in a mixed sample of sex and maturity cohorts across the three study sites. The tracking dataset included 85 turtles at Port Curtis (mean±SD CCL= 84.78±21.66 cm), four at Shoalwater Bay (mean±SD CCL= 92.93±2.33 cm) and five who migrated to foraging grounds from Raine Island (mean±SD CCL= 102.89±2.72 cm, Appendix 2). The Port Curtis turtles comprised of 35 females, 36 males and 14 unidentified sex. Eighteen were juveniles, 24 were subadult turtles and the remainder were mature (49% of turtles in modified habitats were immature). At Shoalwater Bay  one was female andthree were males, only one of the four was a juvenile. The Raine Island turtles were all post-nesting adult females (11% of turtles in unmodified habitats were immature). Capture of turtles, transmitter attachment and data collection for this study were approved by either the JCU Animal Ethics Committee or the Department of Agriculture and Fisheries Ethics Committee and conducted within an approved Queensland Government project.

We discarded the first 24 hours of the tracks from Port Curtis and Shoalwater Bay. We removed duplicates and empty rows, spurious locations determined from unrealistic travel speeds and turning angles (max speed 9.9km/h and threshold inner angle 90 degrees), and locations falling above the high tide line using the SDLfilter R package (25). We treated tracks as independent when a single track contained gaps of more than 72 hours between successive locations, and only retained independent tracks consisting of more than 20 locations. We used both Fast-loc GPS and ARGOS locations where available, but only retained ARGOS points with location classes 1, 2 and 3. We removed nesting and post-nesting tracks of Raine Island turtles by visually identifying the end of directed travel from Raine Island. We also removed possible breeding migrations of Port Curtis turtles identified as a long-distance (greater than 150 km) departure from their release site to known breeding sites for the species during summer months. The filtered tracking data contained 21,578 ARGOS and 56,409 Fast-loc GPS locations.

To identify the environmental characteristics that make foraging habitat suitable for green turtles we developed telemetry-based habitat models using the tracking data. These models combine species presences from telemetry data, and pseudo-absences generated by models detailed below, with environmental data to first identify environmental predictors underpinning species presence and then predict distributions in space and time. We adopted the approach of Hazen et al. 2021 (26) whereby telemetry-based habitat models are constructed using interpolated locations of tracked animals as presences, and pseudo-absences are generated using correlated random walks. Though similar to species distribution modelling, telemetry-based habitat modelling does not use abundance or true presence data to simulate a realized distribution of a species. Instead, presences derived from tracking are spatially and temporally autocorrelated and therefore spatially biased, and pseudo-absences are simulated to control for these biases. For example, pseudo-absences can be generated from simulations that mimic the movement process so that both presences and pseudo-absences have the same autocorrelation structure. The model output reflects the habitat conditions that best explain turtle presence.

a)   Presences and pseudo-absences

We derived presences from the tracking data by interpolating locations at regular intervals, accounting for telemetry error, with time-varying move persistence models in the aniMotum R package (27). We have attempted to reduce autocorrelation and imperfect detection biases with state space models (movement persistence models), which interpolated the irregularly sampled data points to produce an estimated location at at least 12-hour intervals (27,28). We assumed the data had symmetric spatial error, using aniMotum’s default error margins for ARGOS locations, and estimated error from (21) for Fast-loc GPS. For each independent track, we examined state-space models with diagnostic visualisations of 1) a time-series, b) qq-plots and c) autocorrelation functions. By visual examination of candidate state-space model outputs we determined the most appropriate timestep, 12 hours where possible, or lower resolution models for tracks with sparser data (24 or 48 hours). Sparse data resulted in non-convergence or unrealistic patterns (straight lines or perfect loops), which were discarded.

We generated pseudo-absences via correlated random walks (CRW, i.e., where the animal could have gone based on track simulations defined by step lengths and turning angles from the interpolated tracking data; each simulated point is overlayed with the ‘worldHires’ map from the mapdata R package (29)and if the point is on land it is not used and a new point is generated,  Appendix 2) as per (26). CRW models were used to infer drivers of turtle presences on local scales. Our approach assumes presences and pseudo-absences are independent, despite tracking data being inherently autocorrelated in space and time. Turtles tracked for longer are therefore more represented in the presence/pseudo-absences, as are turtles tracked in Port Curtis compared to other sites (Appendix 2).

b)   Environmental predictors

To investigate the relationship between probability of turtle presence and environmental conditions, we developed a set of environmental predictor layers from which to ascertain the conditions at each tracked location in space and time. We catalogued freely available environmental datasets with spatial coverage of the coastal Great Barrier Reef region, and selected variables with known relevance to green turtle foraging ecology based on current literature, with at least 4 km spatial resolution, and appropriate temporal coverage of the tracking period (static datasets developed since 2010 or time-varying spatial data covering late 2010 to late 2019). Selected datasets are summarized in Table 1.

We sampled values of each environmental predictor for each presence and pseudo-absence point corresponding to their timestamp and location. For static datasets (a single, non-time-varying layer) (Table 1) we used terra::extract() for grid data, and sp::over() for shapefile (categorical) data. We extracted values of eReefs variables via a custom R script to query the AIMS THREDDS server. We downloaded annual mangrove maps (30) from the Digital Earth Australia data cube using Dask and used these to generate annual maps of Euclidean distance to pixels of at least 20% canopy cover [i.e., minimum class where mangroves were present in the mangrove dataset (30)] in QGIS. Additionally, we generated 30m (to match the high resolution of bathymetry data) grids of distance to recreational boating facilities, distance to coast, distance to reefs, and distance to rivers in QGIS (details in Table 1)]. Prior to modelling, we evaluated the temporal coverage of the variables derived from the eReefs biogeochemical model with data hosted on the THREDDS server only being available from late 2010 until April 2019. We therefore removed presences and pseudo-absences corresponding with timestamps after April 2019, corresponding to entire tracks of 22 turtles that were tracked in early 2010 or after April 2019. The final collection of turtle tracks used for model development is summarized in Appendix 1.

c)   Modelling approach

We used boosted regression trees to develop the telemetry-based habitat models for 1) modified (Port Curtis) and 2) unmodified (Shoalwater Bay and foraging sites of Raine Island turtles) sites separately, with the gbm.fixed() function of the dismo R package (67). These models use presences and pseudo-absences as the response, and the environmental values as predictors. This approach assumes that presences and pseudo-absences are not autocorrelated (26,68). For each model formulation we used five-fold cross validation, whereby the dataset was split into five parts, with each part acting as a test set for one model, and the remaining parts as the training set. We used a tree complexity of five, learning rate of 0.005, a maximum of 10000 trees, and a bag fraction of 0.75. We plotted the relationship of each predictor in the input data with the response, to inform our specification of the var.monotone argument of gbm.fixed(). The initial folds were generated including all of the input predictors and a large (>10,000) number of trees. We then re-ran the cross-validation, removing those predictors that were relatively unimportant (in the lowest 10% in every fold) to simplify the model and pruned using the out of bag method to adjust the n.trees parameter to encompass the best number of iterations of all the folds. Variable importance is assessed as the total contribution by each variable to a reduction in the loss function.

We evaluated model performance with averages of the performance metrics from the five folds: predictive skill with area under the curve (AUC), true skill statistic (TSS), true positive rate (TPR) and explanatory power with percent deviance explained. The final model was trained on 100% of the input data with the refined hyper-parameters obtained during the cross-validation process.

To investigate environmental predictors in relation to foraging habitat for green turtles, we obtained a ranking of the most influential environmental predictors. We examined the shape of the relationship with the response for each predictor with variable importance greater than 100/number of predictors in the model, with partial deviance plots. Influential two-way interactions were identified with gbm.interactions().