1 INTRODUCTION

Within the marine environment, the integration of commercial fishing vessel tracking systems (e.g. Vessel Monitoring and Automatic Identification Systems) with animal movement data derived from tracking technologies, is beginning to reveal the extent to which sharks are thought to overlap with commercial fisheries (Queiroz et al., 2016, 2019; White et al., 2019). Marine protected areas (MPAs) are often established to exclude commercial fishing, attempting to reduce the mortality of sharks and other animals, by protecting key aspects of their ecology (e.g. core space use areas, breeding grounds). While the perceived ‘success’ of MPAs for conserving sharks and rays is very much location and goal specific (MacKeracher, Diedrich, & Simpfendorfer, 2018), some MPAs have been shown to benefit tropical shark species conservation (Knip, Heupel, & Simpfendorfer, 2012).

Large, remote and strongly protected MPAs are on the increase (Lubchenco & Grorud-Colvert, 2015). Yet even within high profile MPAs, illegal, unregulated and unreported (IUU) fishing activity still persists (Ferretti, Curnick, Liu, Romanov, & Block, 2018; Graham, Spalding, & Sheppard, 2010), in part because the enforcement of such MPAs faces difficult and long-term financial and logistical constraints (Bradley et al., 2018). Furthermore, those that enforce MPAs rarely have the capacity, or indeed the data, to incorporate baseline ecological information to optimize patrolling decisions, as this too takes time, effort and additional funds. The underlying spatial dynamics of animals within MPAs however, where possible, should be incorporated into the management strategy, particularly in instances where monitoring of target species (e.g. animal tracking data) overlaps with monitoring of illegal activity (i.e. sightings, arrests, seizures).

Established in 2010, the British Indian Ocean Territory (BIOT) in the Central Indian Ocean remains one of the world's largest ‘no take’ MPAs. It was designated not only to protect the near-pristine reefs and associated fish communities of the Chagos Archipelago, but also large expanses of the pelagic environment and numerous oceanic features such as seamounts that attract aggregations of highly migratory species (Koldewey, Curnick, Harding, Harrison, & Gollock, 2010). Both before and since the establishment of the BIOT MPA however, widespread IUU fishing activities have operated there. Acknowledging that some species of shark product fetch more at market than others, reef sharks and large pelagic sharks are primarily targeted, but not necessarily specific species (Moir Clark, Duff, Pearce, & Mees, 2015). Poaching by unlicensed vessels, almost exclusively by vessels of Sri Lankan and Indian origin (Moir Clark et al., 2015), resulted in declines in reef sharks in BIOT of nearly 90% between 1975 and 2006 (Graham et al., 2010; Sheppard et al., 2012). The MPA is currently patrolled by a single vessel, managed by the Marine Resources Assessment Group (MRAG) on behalf of the UK Foreign and Commonwealth Office (FCO). From seizures of IUU catch compiled during these patrols between 2002 and 2015, it is estimated that two of the most abundant predators on the reef, grey reef sharks Carcharhinus amblyrhynchos (hereafter GRS) and silvertip sharks Carcharhinus albimarginatus (hereafter STS), on average (with 95% confidence intervals) are extracted at rates of 1,973 (CI: 1,065–3,710) and 2,124 (CI: 1,253–3,640) individuals per year respectively (Ferretti et al., 2018).

With a wide but patchy distribution across the tropical Indo-Pacific, GRS and STS frequently comprise a large component of the predator assemblage within reef ecosystems (Barnett, Abranteská, Seymour, & Fitzpatrick, 2012; Bond, Tolentino, Mangubhai, & Howey, 2015). While these species are similar in size and appear sympatric in many locations, including BIOT, there is increasing evidence that they partition food resources (Curnick, Carlisle, Gollock, Schallert, & Hussey, 2019; Frisch et al., 2016), segregating spatially in both horizontal and vertical planes (Espinoza, Lédée, Simpfendorfer, Tobin, & Heupel, 2015; Heupel, Lédée, & Simpfendorfer, 2018). However, how the dynamic movement patterns of reef sharks influence their susceptibility to IUU activity remains unknown. A recent reconstruction of historical population changes for these species in the BIOT MPA, indicates that their population trajectories have differed considerably in response to commercial fishing and subsequent protection (Ferretti et al., 2018). Further investigation of the link between reef shark behaviour and vulnerability therefore is warranted. The persistence of a species reflects both its demographic characteristics, in particular its intrinsic rate of population growth (Cortés, 2016), but also its exposure to sources of mortality through patterns of habitat use and access to refuges (Dulvy et al., 2014). As refugia availability changes with anthropogenic influences such as management (positive) or increased exploitation and habitat degradation (negative), it becomes increasingly important to understand how sharks behave over different spatial and temporal scales; this will inform not just how they respond to particular conditions and habitat availability but also allow comparisons both inside and outside of MPAs (Lea et al., 2016; Queiroz et al., 2016).

Here we extend the application of telemetry-derived network analyses to explore the relative movement patterns between these sympatric reef-associated species in BIOT. Specifically, we address the following questions: (1) Do the movement networks of these species differ in relation to dispersal and connectivity? (2) How predictable are movement patterns between monitored locations? (3) Are there hotspots of movement within the MPA where enforcement should be directed to best mitigate shark removal? and (4) How do space use strategies, including both movement and residency behaviour of each species, influence their risk to IUU? To address Q1, we quantify movement using network analyses of acoustic detections for both species monitored across a multi-atoll acoustic array. Then we model the movement patterns of each species using an individual-based model (Q2 & Q3) and measure estimates of their space use, residency and movement connectivity in relation to spatially explicit data on IUU activity derived from eight years of patrolling effort within the MPA (Q4).

2 MATERIALS AND METHODS

To address our research aims, these methods are divided into five sections. First, we describe the tracking of reef sharks inside the BIOT MPA which consists of a number of monitored atolls and reefs separated by large expanses of open water, as well as how the data were manipulated. The subsequent sections then address the four specific research questions outlined above.

2.1 Shark tracking and data manipulation

Acoustic tags were surgically implanted in 101 sharks (47 GRS and 54 STS) during February–April between 2013 and 2016. Sharks were tracked for 3 years across five isolated atoll systems in the north of the BIOT MPA. The study array consisted of a total of 52 Vemco acoustic receivers (Vemco, Nova Scotia, Canada, 46 VR2Ws, plus four VR2W-UWM and two VR4Gs) with additional receivers placed further afield in the archipelago excluded from the analyses due to later deployment dates or lack of data from this cohort of animals (Figure 1a). Sharks were caught at various locations around Peros Banhos and the Salomon Islands using hand-lines with barbless circle hooks and brought alongside or into the boat. On the boat, a seawater hose was inserted into the mouth to irrigate the gills, with a cloth covering the eyes to decrease stress on the animal. Vemco V16 acoustic tags were soaked in betadine, then surgically implanted into the shark by making a small incision in the abdominal region of the animal and inserting the tag into the peritoneal cavity. The incision was then closed with a biodegradable suture. Sharks over 150 cm were brought alongside, restrained with a tail rope, rolled over to place in tonic immobility (Kessel & Hussey, 2015) and then tagged using the same procedure. Animal handling procedures were approved by the Stanford University Administrative Panel on Laboratory Animal Care under permit APLAC-10765, held by Hopkins Marine Station's Block Laboratory and the Zoological Society of London's ethics committee.

Tags acoustically transmit a unique ID code at regular intervals (nominal delay of either 30–90 s or 60–180 s) for the duration of their battery life (~10 years). Tagged animals were detected whenever they came within range (~500 m) of an acoustic receiver. As focal receivers were deployed periodically, inter-annual analyses were performed to test the potential influence of increasing receiver numbers and animal loss on the system, but no significant differences in connectivity year-on-year were found, allowing the data to be considered in its entirety (see Appendix S1; Figure S1). Spurious detections caused by infrequent tag clashes or ambient noise at 69 kHz, were removed and detections were standardized to 1.5 min to achieve comparable detection frequencies (see Appendix S1 for correction procedure). Finally, for each individual, data for the first 24 hr post-tagging were removed to make sure that the detected behaviour was not influenced by the stress of capture (Cliff & Thurman, 1984).

2.2 Quantifying connectivity in movement networks

Movement networks, based on movements between receivers in either direction, were used to define individual or collective space use and movement patterns of the two species. Network ‘nodes’ were the acoustic receiver locations and the ‘edges’ between nodes marked the departure of individual sharks from one location (origin) followed by the arrival at another (destination), the next receiver the shark was detected on. Thus, weighted networks were created where edge weight represented the frequency of transitions between receiver locations (Jacoby, Brooks, Croft, & Sims, 2012). From structural metrics extracted from the movement networks (below), connectivity was defined for both the full network (all tagged sharks of both species) and separate species networks to compare movement patterns between these sympatric species. We compared the number of transitions by each species using a Wilcoxon rank sum test. Structural metrics included:

reciprocity, 〈r〉 for directed networks defined as

where is the edge weight or frequency, at the population level, of directed movements between receivers i and j, is the frequency of these movements that are reciprocal in both directions. For a comparison of how well receivers were connected to one another across the full and species-specific movement networks, standard centrality metrics, mean degree or the unweighted sum of nodes connected to a specific location, can be expressed as

where

n is the total number of nodes and the binary incidence of an edge between receivers i and j; node strength 〈s〉 which is a summed or weighted (w) version of 〈k〉 where

Edge density, D gives the proportion of edges present in the symmetric, undirected network relative to the maximum possible number of edges with a given number of nodes (n), thus

D enables us to explore how receiver array layout, in addition to shark movement strategies, might play a role in constraining connectivity within the system. Finally, a weighted clustering coefficient, C (Barrat, Barthélemy, Pastor-Satorras, & Vespignani, 2004) was calculated

where s_i is the strength of node i as defined above and i, j and h are all network neighbours. C is a measure of how strongly local nodes are connected to one another relative to those in the rest of the network.

2.5 Space use and IUU risk

2.5.1 Activity space estimates

To estimate species-specific activity spaces, we adapt a spatial modelling approach often used to infer animal density from camera trap arrays (Royle, Karanth, Gopalaswamy, & Kumar, 2009). This approach was adopted over other more conventional methods (e.g. kernel utilization distributions) to enable inferences of activity spaces beyond the boundaries of the array. Briefly, this model assumes that each individual has a hypothetical centre of activity, c_x, at the core of their activity space. As the distance from c_x increases, the frequency of occurrence of the individual at any point in space declines.

To capture this process, we modelled the detection frequencies of the individual sharks recorded at the receivers (as a proxy of the frequency of occurrence). We expected that the hypothetical c_x of each shark was nearest to the receiver with the highest overall detection frequency (r_c). As the distance from this receiver increased, we assume the expected detection frequency (r_i) at the other receivers declines, most likely with an exponential decaying function. Detection frequencies for this model were defined as the number of detections per day. This index was calculated by dividing the number of detections recorded at each receiver by the days this receiver could detect pings from the acoustically tagged shark (i.e. the number of days between the first and last detection for an individual and for the duration the receiver was active). Receivers were accumulated over time and therefore were not all active for the same period (see Appendices S1 and S2).

To derive a boundary estimate of activity space we set a threshold of detection frequency beyond which we assumed detections were sporadic incursions outside the usual activity space of the animal. We arbitrarily set this threshold at 0.0054 (i.e. ≤2 detections per year) deciding that if a location had an expected detection rate of fewer than two pings in a year, this was not part of that animal's regular activity space (see Appendix S2 for sensitivity analyses around this decision). So our modelled boundary estimate of activity space was the distance that would predict a detection rate of 0.0054 pings per day.

We assumed that the detection frequencies followed a normal distribution when log-transformed, with mean as a declining function of the distance between r_c and receiver r_i.

To account for variation between individuals within each species, individual was included as a random effect both in the intercept and the slope of the modelled relationship, although the random slope did not account for shark length. A random effect on the intercept would capture the variability expected in the distance between the unobserved c_x and the receiver r_i with the highest detection frequency (i.e. the closer a receiver is to an activity centre, the higher its expected detection frequency). A random effect on the slope would capture the possible intraspecific variability in movements that generate an activity space. For each species, we estimated the decaying function by fitting linear mixed-effects models of the structure:

where d is the detection rate, X is the matrix of predictor variables including log-distance and total shark length, is a vector of fixed effect coefficients, Z is the design matrix of the random effects, is the vector of random effects and is the vector of residuals, which are assumed to follow a normal distribution; models were fit using the lme4 r package (Bates, Mächler, Bolker, & Walker, 2015). Total length (TL) of sharks varied (range = 80–208 cm) and therefore we also tested this covariate as a fixed effect in the activity space analyses. Different model structures were compared according to their akaike information criterion. Finally, sensitivity analyses were conducted to explore the potential impact of different minimum detection thresholds, and of removing receivers, on activity space estimates (see Appendix S2; Figures S2 and S3).

2.5.2 Spatial risk

MRAG are a consultancy organization that manage and enforce the MPA using the patrol vessel, and prior to that the fishery in BIOT for the UK Government. Intercepted illegal fishing events are documented by MRAG with details on the location, boat characteristics, seized catch and other incidental information useful for developing a map of IUU hotspots in the MPA. MRAG score each 0.1° × 0.1° grid cell (~11 km²) throughout the BIOT MPA to reflect the cumulative likelihood of illegal activity. Scores combine measures of risk based on previous known illegal activity (historical arrests), expert knowledge of MPAs, whether an area is deemed a priority based on current UK Foreign and Commonwealth Office (FCO) policy, as well as information on distance to MPA boundary, current transit routes and vessel sightings (Davies, Moir-Clark, Pearce, Mees, & Collins, 2014). They assign each potential risk a likelihood of occurrence and an impact score (both scale 0–3) and values were accumulated to produce a risk per cell between 1 and 50 (low risk to high risk, summarized in Appendix S2; Table S1).

To assess the impact of reef shark movement dynamics on potential risk of IUU fishing, movement networks were generated for each individual for each species and analysed as a collection or stack of undirected, weighted networks. Mean (log-transformed), weighted degree, 〈s〉 for individual movement networks were modelled with a GLM against an individual composite risk score. Using the MRAG spatial risk data, composite risk reflected individual risk associated with residency patterns across all encountered locations, in addition to the potential risk individuals face as part of their wider activity space when not detected by our receivers. For each individual, a polygon with a radius given by our activity space estimate was created, centred on the receiver with the highest number of detections. A mean risk was extracted across all grid cells within this polygon using the packages rasterVis, rgdal and dismo in r. For residential risk, individual residency per location was derived from the number of hours in a day in which an individual was detected at a particular location, relative to the time at liberty (TAL), in hours across the array. This resulted in values between 0 (no residency) to 1 (continuous residency). Residential risk reflected the time spent within the grid cell that the receiver was situated in. Finally, the composite score simply weighted the residential risk and the activity space risk to reflect the proportion of time an individual spent inside and outside the array. To account for individual variation in residency, the response variable in the GLM was weighted by the proportion of TAL spent within the array. Unless stated, all analyses were conducted in the stats and igraph packages in r version 3.2.3 (Csardi & Nepusz, 2006; R Core Team, 2015).

3 RESULTS

3.2 Predictability of movement patterns

In exploring the predictability of movement patterns, the distributions of step lengths for both species were found to follow a lognormal function (Figure 2). The distance-decay IBM revealed a potential predictability of 96% for STS and 99% for GRS (i.e. the proportion of observed, weighted edges that fell within the complete range of the modelled values). The model appeared to capture GRS movement considerably better than STS, under-predicting GRS movements for 13 receiver pairs out of a possible 1,362 (p < 0.05). The IBM for STS underestimated more edges, a total of 35, and overestimated 13 at the p < 0.05 level.

3.4 Space use and IUU risk

Using a minimum detection threshold of ≤2 detections per year and controlling for TL of individuals (TL range: GRS 70–171 cm, STS 83–208 cm), GRS showed a significantly smaller estimated activity space radius of 58 km (34–99) than STS at 175 km (112–275) (Figure 3). Activity space estimates declined steadily for each species with increasing minimum thresholds but importantly the relative difference between species remained fixed across detection thresholds (Figure S2).

Composite risk to STS, a score reflecting risk associated with residency, movement and activity space, was found to increase significantly with increased spatial connectivity, measured as the mean weighted degree of the network nodes 〈s〉 (slope = 1.726, R² = 0.286, p < 0.001). No relationship was found in GRS (slope = 0.691, NS, Figure 4). We reran the model removing the large outlier for STS (notably the individual for which we had the most data), and found this relationship to still hold, albeit with the shallower slope (STS outlier removed: slope = 0.700, p = 0.022).

4 DISCUSSION

Understanding the habitat use and distribution of sharks in remote coral reefs is vital for assessing the protection afforded by the creation of large MPAs as well as for measuring the impact of threats posed by anthropogenic pressures. This study builds on the work of Espinoza, Lédée, et al. (2015) and Lea et al. (2016) who utilize network analyses to explore the residency of reef-associated shark species between managed and unmanaged areas, in addition to Queiroz et al. (2016, 2019) and White et al. (2019) who explore the overlap between shark movements and commercial longline fisheries. Explicitly, we demonstrate that, despite their sympatry, these species connect the MPA in very different ways. Further, the movements of tagged sharks between receiver locations in the north of the MPA at least, are quite efficiently predicted using very simple distance-decay models. This means that with knowledge of an individual's core area, we can estimate how far they are likely to disperse based on short distance movements between locations being much more likely than long-distance ones. We use outliers from the IBM to infer important movement corridors where ecological conditions at origin or destination receivers, or indeed in between, drive unexpectedly high levels of movement, that can be used in the short-term to update the risk map that currently informs routine enforcement patrols. Our detailed analysis of movement networks corroborate previous work suggesting that STS are more wide-ranging and connective of reef ecosystems than GRS (Carlisle et al., 2019; Espinoza, Heupel, Tobin, & Simpfendorfer, 2015; Espinoza, Lédée, et al., 2015). We extend this work however, to show that by occupying broader activity spaces, being less resident and moving more frequently between a greater number of areas, STS are more vulnerable to IUU risk.

So what is driving these predictable patterns of movement and these apparent ecological hotspots? Due to time constraints in territory and the area covered by each of our receivers, our analyses currently do not integrate information about the environmental variables at each monitoring station (i.e. node covariates); further research that incorporates high-resolution satellite remotely sensed environmental data will be important for exploring the drivers of reef shark movements in these remote locations (Williamson, Tebbs, Dawson, & Jacoby, 2019). That said, a recent study exploring the bleaching and recovery of the scleractinian coral reef communities within BIOT during the consecutive, temperature-driven mortality events of 2015 and 2016 (Head et al., 2019), showed that our predicted ecological hotspots had lower percentage change in coral cover in response to warming events (e.g. Northwest Peros Bahnos), compared to other areas. Thus, in order to better tackle the logistical, financial and management challenges associated with the BIOT MPA, we argue for greater integration of concurrent environmental monitoring data (both in situ and remote) with continued tracking of reef-associated predators, all contributing to a frequent dialogue between science and policy (Bradley et al., 2018).

Important movement corridors likely also prove productive grounds for illegal fishing activities. Here we suggest that the IUU fishing affecting the MPA has greater implications for STS that move more frequently and more widely. Even as juveniles, silvertip sharks have a tendency to associate strongly with reef drop-offs, a characteristic that potentially contributed to the loss of >20% of the tagged individuals caught by fisherman around an MPA in the Seychelles (Lea et al., 2016). Indeed, further exploration of our own time-series data have revealed a cessation of detections from 15 acoustic tags over a period of just 10 days in December 2014 in a suspected illegal fishing event, 13 of which were STS and just two were GRS (Tickler et al., 2019). As further information is gathered on illegal vessel sightings and arrests, vulnerability can be more directly measured through the integration of enforcement vessel positional data (e.g. VMS) and shark movement data (Queiroz et al., 2016) that better incorporates the potential temporal variation in IUU threat in BIOT (Davies et al., 2014).

Our receiver array, although covering a geographic area of approximately 1,500 km², unlikely captured the full dynamics of shark movements despite this cohort of animals not being detected on any other receivers (e.g. in the south of the archipelago, ~160 km). Therefore, to estimate species activity space, we used parametric modelling approaches that were better suited than conventional kernel utilization distribution, to extrapolate beyond the geographic restrictions of the receivers analysed. Specifically, we adopted an approach previously used to infer animal density from camera trap arrays (Royle et al., 2009). Our GRS activity space estimate of 58 km, is relatively consistent with White et al. (2017), who found that the median distance moved beyond the boundary of Palmyra Marine Reserve by satellite-tagged GRS in the Pacific Ocean was 66.2 km, helping to validate our model for this species at least. In addition to a significantly larger activity space, STS were less resident and undertook many more movements to and from different places, exemplified by high mean degree and high clustering of the movement network, all behaviours that contributed to their increased exposure to IUU activity. The less expansive movements of GRS here reflect the central place foraging behaviour of this species observed elsewhere (Heupel et al., 2018; Papastamatiou et al., 2018).

The BIOT MPA is currently enforced by a single, multipurpose patrol vessel, that has been previously roughly estimated to detect just 10% of illegal fishing activity (Price, Harris, McGowan, Venkatachalam, & Sheppard, 2010). Patrols to intercept illegal fishing in BIOT are routinely guided by anthropogenic evidence, including previous arrests and radio intelligence, or by covering key geographic features (e.g. seamounts) within the MPA on a semi-regular basis (Davies et al., 2014). Finding an optimal search strategy for the enforcement vessel however, is complicated by the fact that the locations in which seizures are made are not necessarily the same locations in which the sharks are caught. Although the two species analysed in this study are not the sole target for illegal fishers, they are regularly extracted in large numbers (Ferretti et al., 2018; Tickler et al., 2019). Our findings that shark movements are actually quite predictable at least tells us where we might expect to find them in the future, information that can be immediately used by MPA managers. Due to the difficulties of obtaining data on illegal activity in BIOT (e.g. time and location of arrests), a limitation of this study is the temporal aggregation of the risk landscape, although MRAG do update their risk map on a monthly basis.

In conclusion, understanding the movements and predictability of reef predators in remote oceanic MPAs is important for explaining community structure and for predicting the outcomes of conservation measures, recovery plans and enforcement strategies. Through long-term acoustic monitoring of reef sharks and network analyses, we explore how dynamic movements between receivers can highlight important, species-specific ecological corridors in addition to enabling us to estimate species risk to IUU fishing activity. In highly dynamic coral reef environments, which are increasingly subject to substantial and sometimes irreversible changes (e.g. bleaching, destructive fishing techniques, temperature fluctuations), it is crucial to understand how reliant species are to localized sections of an ecosystem and therefore how resilient they might be to putative disturbance and exposure to spatially variable threats.

AUTHORS' CONTRIBUTIONS

A.B.C., T.K.C., D.J.C., J.J.D., R.J.S., F.F. and D.T. collected the field data; B.A.B. led the acquisition of data; D.M.P.J. conceived the study, while D.M.P.J., F.F. and R.F. analysed the data; D.M.P.J. wrote the manuscript and all authors contributed and agreed to the final version.