Cross-jurisdictional larval supply essential for eastern Australian spanner crabs (Ranina ranina)

Keywords: connectivity, East Australian Current, Lagrangian, larval subsidy, modelling, particle tracking, pelagic larval duration, temperature dependent growth.

Variability in the productivity of marine fisheries is primarily driven by dynamics of the exploited stock (Anderson et al. 2008). It is increasingly being recognised that environmental conditions have a significant influence on larval survival, dispersal and recruitment to fisheries (Dixon et al. 1999; Szuwalski et al. 2015). However, the mechanistic relationships between environmental and population dynamics are often poorly understood. This limits our ability to model variation in productivity, which has implications for the sustainable management of exploited stocks (Yatsu et al. 2005; Kurota et al. 2020).

For many exploited marine species, dispersal is driven by ocean currents during a pelagic larval stage. For such species, connectivity patterns are defined during the larval stage and can cross the jurisdictional boundaries within which exploited species are typically managed (Everett et al. 2017). It is difficult to convincingly measure larval connectivity in oceanic systems (White et al. 2019). Genetic techniques are often employed to examine genetic transfer between regions (Lowe and Allendorf 2010; Saenz-Agudelo et al. 2011). Unfortunately, genetic techniques provide little insight into the dispersal mechanisms of the larvae, which are often driven by oceanography, and genetic connectivity can be driven by the movement of few individuals (Cowen et al. 2007). Provided appropriate models are available, particle tracking simulations can show the influence of oceanographic processes, including dispersal mechanisms and pathways.

The coupling of hydrodynamic models with particle tracking simulations, which can incorporate species-specific biological rates (e.g. growth, mortality), allows the re-creation of past dispersal events (Gallego et al. 2007; Hinrichsen et al. 2011). Examining these past dispersal events alongside oceanographic conditions can show both patterns and drivers of connectivity and recruitment. Particle tracking techniques have been applied successfully in many investigations at scales from small embayments (Norrie et al. 2020) to global oceans (Doblin and van Sebille 2016). Further, these techniques have been applied to a variety of taxa including microbes (Doblin and van Sebille 2016), kelp (Coleman et al. 2013), invertebrates (Everett et al. 2017; Cetina-Heredia et al. 2019) and fish (Schilling et al. 2020). By modelling larval dispersal of fished species, it is possible to investigate the influence of oceanographic conditions on interannual fluctuations in recruitment and the impacts on variability in fisheries productivity.

Spanner crabs (Ranina ranina) are large brachyurans found throughout the Indo-Pacific region. Like many marine decapods, spanner crabs have a two-phase life cycle with a pelagic larval phase and a benthic post-settlement, non-migratory juvenile–adult phase (Brown et al. 2008; Kapur et al. 2019). The species is commercially harvested throughout much of the Indo-Pacific (Kapur et al. 2019), with the largest fishery occurring in the waters off eastern Australia (Scandol and Kennelly 2002; O’Neill et al. 2010). The stock in Australia is managed independently by the states of New South Wales (NSW) and Queensland (Qld), with the latter accounting for most of the harvest (e.g. ∼80% in 2017).

The combined commercial harvest of spanner crabs from Qld and NSW declined from >4000 to <1000 tonnes between 1994 and 2019 calendar years. The reduction in the overall catch is partially attributed to large reductions in effort owing to management changes in both jurisdictions, reducing the overall fleet size. The Qld fishery is recognised to have declined substantially since the mid-1990s, with reduced harvests and catch rates in both the fishery and fishery-independent surveys, and overall spanner crab biomass is currently estimated to be below target reference levels (Queensland Government 2021). Within the NSW fishery, which accounts for ∼20% of the overall harvest in this stock (Roelofs et al. 2021), the catch has also declined since the 1990s (from ∼465 to ∼110 tonnes in 2018). However, the main performance indicators (standardised catch rates) for the NSW fishery were tracking above long-term averages from 2015 to 2019, although declines are evident in recent years (Roelofs et al. 2021; Supplementary Fig. S1). The total allowable catch in both NSW (169 to 140 tonnes) and Qld (1631 to 847 tonnes) has been reduced in recent years to restrict fishing mortality, increase protection of the spawning stock, and halt the decline in spanner crab stocks (NSW Total Allowable Fishing Committee 2020; Queensland Government 2021). Within NSW, in 2015, a peak in sublegal-size crabs was observed in the independent survey, which corresponded to a peak in commercial catch-per-unit-effort (CPUE) the following year (Johnson 2021). The peak in both subadult and subsequent legal-size crabs was likely to be driven by favourable environmental conditions increasing the catchability of crabs (Spencer et al. 2017), or a strong recruitment event. However, source–sink dynamics and larval dispersal patterns, which may provide insight into recruitment, are unclear for the species.

In general, spanner crabs spawn throughout the fished areas of the species distribution (Kennelly and Watkins 1994; Brown et al. 2008; Zones 2–7 in Fig. 1). The spawning biomass is protected by temporal closures, in both NSW and Qld, which align with periods of peak reproductive activity (November–December–January), as well as minimum-size limits (NSW: 9.3-cm orbital carapace length, CL; Qld: 10-cm rostral CL; these are approximately equal-sized crabs but measured differently) that aim to allow mature individuals to spawn at least once before recruiting to the fishery, while helping keep a higher proportion of females in the population. In eastern Australia, regionally variable environmental conditions, including temperature, chlorophyll-a, bottom current velocity, and wind-driven upwelling, have large influences on catch rates (Spencer et al. 2017, 2019a, 2019b, 2019c). However, little is known about larval connectivity, and resulting recruitment. This study, therefore, aimed to use particle tracking simulations to (1) understand the patterns of larval dispersal and connectivity in the eastern Australian spanner crab stock; and, (2) investigate correlations between predicted settlement rates and NSW commercial CPUE for the species.

Fig. 1.  Map of eastern Australia, showing (a) the six main fishing regions for spanner crabs overlaid on blended depth colours and relief shading of the ocean bottom derived from CleanTOPO2 data. The dashed line represents the 200-m isobath. Zones 2–6 are managed by Queensland (Qld), whereas Zone 7 is managed by New South Wales (NSW), and (b) an example (1 January 1993) of the oceanographic conditions common to this region as modelled in the BRAN2020 ocean model. The vectors show the direction and speed of the currents, and the colour represents the sea-surface temperature. For our simulations, particles were released on the benthos at depths of <200 m in all six zones shown in (a).

The eastern Australian spanner crab stock is primarily located between 23 and 29°S, with a very small amount of fishing occurring north and south of the main area (Fig. 1). Fishing occurs within six main zones, with Zones 2–6 being managed by Queensland Department of Agriculture and Fisheries (QDAF) and Zone 7 managed by New South Wales Department of Primary Industries – Fisheries (DPIF; Fig. 1). Crabs are harvested using flat dillies, designed to entangle rather than trap, deployed in depths between 30 and 80 m (Dichmont and Brown 2010). In NSW, following the introduction of daily reporting (2009–10), fishers have been required to report the number of net-lifts per fishing day, allowing CPUE (kg fisher day⁻¹) estimates to be standardised to the mean number of net-lifts per fisher-day.

The occurrence of ovigerous females in a large-scale, stratified, randomised survey of NSW identified their reproductive period to span from November to January (austral summer), peaking in December (Kennelly and Watkins 1994). In this survey, the estimated fecundity for individual crabs varied from ∼78 000 to 256 000 eggs per batch and was positively correlated with body size (Kennelly and Watkins 1994). Spanner crabs show limited movement following settlement and are estimated to enter the fishery at Age 4 (Brown 1986; Brown et al. 1999, 2008); however, age and growth remain uncertain for the species.

To simulate the dispersal of larval spanner crabs, velocity and temperature fields from the most recent (2020) version of the Bluelink ReANalysis (BRAN2020; Chamberlain et al. 2021) were used. BRAN2020 is a configuration of the Ocean Forecasting Australian Model, ver. 3 (OFAM3) and is an updated version of the previous BRAN models (Oke et al. 2013a, 2013b). A comprehensive description of OFAM3 is presented by Oke et al. (2013b). Briefly, OFAM3 is a near-global, eddy-resolving, z* configuration of the Modular Ocean Model, developed principally for the purpose of reanalyses and forecasting upper ocean conditions worldwide. BRAN2020 takes observations from satellite sea-surface temperature (SST), satellite sea-level anomaly (SLA) and in situ temperature and salinity and uses these to match the model output to real world conditions (full details in Chamberlain et al. 2021).

The analysis used the daily averaged outputs of zonal (east–west, u) and meridional (north–south, v) velocity (m s⁻¹) and temperature (°C) within a subset of the model (10–40°S, 140–165°E). This domain encompasses the entire eastern Australian region from the southern tip of Papua New Guinea (10°S) to the Bass Strait (40°S), stretching east to New Caledonia (165°E; Fig. 1). Within this region the South Equatorial Current flows westward into the Coral Sea and Great Barrier Reef where it splits into equatorward (north) and poleward (south) components (Oke et al. 2019). The poleward component, the East Australian Current (EAC), is a western boundary current and is the dominant oceanographic feature along eastern Australia (Oke et al. 2019). The EAC forms a jet of fast-flowing water along the coastline, which generates upwelling and eddies along its path (Suthers et al. 2011; Cetina-Heredia et al. 2014; Schaeffer and Roughan 2015).

The model has a 0.1° (∼10-km) resolution grid, with a vertical resolution of 5 m in the upper ocean (top 40 m) and 10-m vertical resolution down to 200-m depth (Chamberlain et al. 2021). BRAN2020 is a data-assimilating model that continually adjusts the model to match observations, thereby recreating an accurate representation of past conditions, including mesoscale eddies with high accuracy (Chamberlain et al. 2021). The re-analysis model runs for 26 years (January 1993–December 2019).

Although BRAN2020 is fully validated and performs well at reconstructing mesoscale features (Chamberlain et al. 2021), an important but often overlooked part of using large-scale or global models in specific regions is the validation of the model against local, fine-scale observations. To assess the appropriateness and determine our confidence in the BRAN2020 model, we compared velocities and temperature at 10–30-m depth with observations at 13 ocean-observing moorings in our study region (Supplementary Table S1, Fig. S2). These moorings are independent and not assimilated into BRAN2020. We compared the modelled and observed velocities using Taylor diagrams to infer the relative strengths and weaknesses of BRAN2020 for our applications (Taylor 2001; North et al. 2009).

Lagrangian particle-tracking simulations were conducted using PARCELS (‘Probably A Really Computationally Efficient Lagrangian Simulator’, ver. 2.2.0, see https://oceanparcels.org/; Lange and van Sebille 2017; Delandmeter and van Sebille 2019), run using Python (ver. 3.8.3, see https://www.python.org/downloads/release/python-383/; Van Rossum and Drake 2009). Lagrangian simulations were run using the horizontal velocities from all depth levels, with all vertical movement forced manually to simulate the biology of spanner crab larvae. During the spawning season, we released 1000 particles per day from each of the ocean zones (<200-m bathymetry). The location within each zone was randomised for each day. Because spawning is known to occur on the benthos and megalopa are found at the surface (Brown et al. 2008), all particles were released in the bottom ocean layer and then rose gradually to the surface over the next 10 days. As the swimming ability of spanner crab larvae is unknown, we have modelled all particles with passive horizontal dispersal (no swimming ability). Following a similar procedure to Schilling et al. (2020), the paths of each particle were interpolated using 5-min time steps, by using a fourth-order Runge–Kutta scheme, and recorded once daily. Each particle included a small random diffusion function to add natural variation to the movement of each particle. Following Okubo (1971), the scale of the horizontal diffusion constant (8.8 m² s⁻¹) was chosen to capture the effect of flows occurring at spatial scales smaller than that explicitly solved by the hydrodynamic model (e.g. sea breeze, turbulence). This parametrisation was successfully used with a previous version of the Bluelink Reanalysis Model at the same resolution (Cetina-Heredia et al. 2015). Particles were, therefore, advected by both the BRAN2020 velocity fields and then in a random direction by the velocity scaled by the diffusion coefficient (van Sebille et al. 2018).

To account for the temperature-dependent nature of larval growth rates and pelagic larval durations, we modelled the ‘settlement’ of our larvae using degree-days (DD; Neuheimer and Taggart 2007; O’Connor et al. 2007). This approach calculates the cumulative daily temperatures experienced by larvae, and settlement is deemed to occur when the DD threshold is passed, resulting in faster settlement in warmer water. Experiments have shown that spanner crabs metamorphose to megalopa after 41.3 days at 25°C (Minagawa 1990). This is equivalent to 1032.5 DD, which was used as our threshold for settlement. As a sensitivity analysis, we tested multiple DD thresholds (800, 1000, 1200 and 1400 DD) and visually compared the predicted settlement patterns in our region of interest (NSW; Zone 7). The sensitivity analysis showed that DD thresholds of 200 (∼7–8 days) apart result in very similar settlement patterns, with differences more obvious with larger DD thresholds (Fig. S3). Regardless of the threshold tested, the broad patterns were consistent with all scenarios, showing putative peaks in settlement between 1997–2000 and 2012–2013. A longer threshold reduced the overall number of particles settling but because these are a relative measure, conclusions are unaffected.

We explicitly modelled natural mortality as part of our simulations by releasing many particles and applying a daily mortality rate to each cohort (particles released at the same time and place). Because natural mortality is assumed to be a daily constant for each cohort and each cohort started with the same number of particles, it can be assumed that the number of particles dying each day is equal among cohorts until particles begin to settle and the number of particles within each cohort may begin to vary (Schilling et al. 2020). Following this logic, to save computing time and reduce the required number of starting particles, we applied natural mortality only from Day 34 (prior to any particle reaching 1032.5 DD and settlement occurring). The actual mortality rate for spanner crab larvae in this region is unknown and we, therefore, applied the mortality rate from a laboratory study of spanner crab larvae, similar to Everett et al. (2017), which therefore does not include natural mortality from predation and is a conservative estimate of mortality. Regardless, by modelling a daily mortality rate (even if it is not entirely accurate) with a temperature-dependent growth rate, we are effectively able to model the effect of cumulative mortality on larval cohorts that have a longer PLD and are more vulnerable to predation before settlement. A simple sensitivity analysis showed that the predicted settlement under 10× mortality was highly correlated with the original predicted settlement (r = 0.96, t₂₂ = 17.31, P < 0.001), demonstrating the outcomes in terms of relative settlement do not change in response to the specified mortality (because mortality is a constant). Larvae that did not reach the settlement threshold while on the continental shelf (<200-m bathymetry) within a recognised fishing zone were considered mortalities in terms of the analysis (dispersal mortality; Everett et al. 2017; Schilling et al. 2020). Owing to the course grid cells (∼10 km) of BRAN2020, we assumed that settlement in <200-m bathymetry was representative of settlement in the 30–80-m bathymetric range.

Assuming equal spawning output from each of the zones, we released 1000 particles per day in each zone during the spawning season (1 November−31 January) for 26 spawning seasons (1993–94 to 2018–19). Spawning seasons were identified as the year in which spawning began. On each day of the spawning season, a random location in each zone was selected and all particles were released there. This stratified random-release strategy ensures that over the spawning season there is a spread of spawning from shallower and deeper areas. All particles were then tracked for 60 days to ensure that they had sufficient time to reach the temperature-dependent settlement threshold.

Connectivity among the six fishing zones was assessed by two connectivity matrices, which can be viewed as a method to determine source–sink relationships within a population (Mitarai et al. 2009; North et al. 2009). For each spawning season, we calculated the percentage of larvae settling in each of the fishing zones (i.e. sinks) from each of the spawning regions (i.e. sources) as

where PS_ifh is the percentage settlement in Zone i, for particles released in Zone f in Spawning season h, N_ifh is the number of particles settling in Zone i, released in Zone f in Season h, T_fh is the total number of particles released in Zone f in Spawning season h.

We also calculated a second metric regarding the origins of particles landing within each zone as

where G_fh is the total number of particles landing in Zone f in Spawning season h. Averaging the percentage particles settled in each source–sink combination over the 26 spawning periods gives relative metrics (mean settlement) of connectivity among the different fishing areas. We also calculated the standard deviation of this mean, which provides information regarding the variability in the connectivity among the zones. For each settlement zone, the first matrix shows the ‘origin percentages’ of particles settling in each spawning zone. The second matrix shows the ‘destination percentages’ of all spawned particles in each of the zones.

To test the hypothesis that strong recruitment events may be due to favourable oceanography within the NSW fishery, we quantified the number of particles predicted to settle in the NSW fishing zone (Zone 7) and assumed that this was a proxy for larval recruitment in each spawning season. This assumes equal spawning output each year across the stock, which is not realistic but without detailed information regarding the distribution and abundance of spawning crabs, this is a necessary assumption to examine our specific hypothesis. No detailed information exists regarding the distribution or interannual variation in the spawning stock biomass. We therefore expect that our metric at least partially captures overall recruitment, acknowledging that we are not accounting for variations in spawning stock biomass. Because the fishery comprises multiple year classes, we calculated a 3-year rolling average to our settlement predictions (to capture a main age class and the fast or slow growers from the following or preceding age classes). This was then lagged by 4 years, matching the broad assumption of the age at which male spanner crabs are estimated to recruit to the fishery (Brown 1986; Brown et al. 1999, 2008; Chen and Kennelly 1999). Because adult spanner crabs are non-migratory (Brown et al. 2008; Kapur et al. 2019), we consider this metric a good proxy for recruitment to the fishery, and suitable to assess how recruitment variability can contribute to variations in fisheries productivity. We, hereafter, refer to this metric as the ‘predicted recruitment index’.

To represent variation within the fishery, two CPUE metrics were calculated from NSW logbook data, both reported as the annual (July–June) average. The first source is a metric of kilograms of harvest per fisher-day, which extends from 1998–99 to 2019–20 (hereafter ‘long-term CPUE’). The second source is a metric of kilograms per net-lift and ranges from 2009–10 to 2019–20 (hereafter ‘short-term CPUE’) and which acknowledges that fishing effort (number of net lifts) may vary among days. The two CPUE metrics are highly correlated (r = 0.82, t = 4.24, P = 0.002). Within the NSW context, previous investigation showed no difference between using standardised catch rates or raw catch rates, so we used raw catch rates to maximise relevance to the fishery (Johnson 2022; detailed below). We used Pearson correlation tests to quantify how well our predicted recruitment index matched the observed fluctuations in each CPUE metric.

Standardised NSW catch rates predicted from generalised linear models (GLM) are presented in the NSW stock assessment for spanner crab (Johnson 2022). Explanatory model terms considered different catch rates among fishing years, seasons, individual fisher operations, their transformed fishing effort (the number of net-lifts, which was a function of the number of ground-lines used, nets per ground-line and ground-line lifts per day; log or cube-root scale), the spatial locations of catches based on 6 × 6′ latitude and longitude grids. However, there is minimal difference between NSW raw and standardised catch rates (D. D. Johnson, unpubl. data), likely owing to stable management arrangements that were consistent, with no changes in permitted effort (capped at 30 dillies) or vessel capacity from 1997–98 to 2017–18 in NSW (https://legislation.nsw.gov.au/view/html/inforce/current/sl-2006-0738#sch).

In addition to the fishery-dependent CPUE, we also used a fishery-independent survey dataset that controls for vagaries of fisher behaviour, and has been identified as the best data source to investigate environmental influences on recruitment (Filar et al. 2021). Briefly, the survey is conducted annually each year (since 2005 in NSW), stratified by Fishing regions 2–7 (for full details of this survey, see Department of Primary Industries and Fisheries 2005; McGilvray et al. 2006). To focus again on NSW, we used data only from Region 7, which consisted of data from 2005 to 2018 inclusive. Data were extracted from Filar et al. (2021) by using the R package ‘digitize’ (ver. 0.0.4, see https://CRAN.R-project.org/package=digitize; Poisot 2011), and reported as the number of crabs per ground-line. The data were reported as sublegal crabs (<100-mm rostral carapace length), legal crabs (≥100-mm rostral carapace length) and all crabs combined. The legal-size crab data were correlated with both fishery-dependent metrics (kg fisher-day⁻¹: r = 0.84, P < 0.001; kg net-lift⁻¹: r = 0.64, P = 0.046).

Spawning intensity is unlikely to be evenly distributed throughout the spawning season, and there is a clear peak during December (Kennelly and Watkins 1994). To investigate whether spawning in particular months was driving the relationship between our predicted recruitment index and CPUE, or whether spawning times may have changed during the past 25 years, we repeated the above analysis monthly from November to January. We hypothesised that the December relationship would be the strongest on the basis of the previously observed peak in spawning.

In our simulations, we assume that spawning is equal across space; however, this is unlikely to be the case. To assess the effect of violating this assumption, we conducted a sensitivity analysis where we sequentially removed spawning (i.e. released no particles) from within one region at a time. For each of these simulations, predicted (simulated) recruitment to the fishery in NSW was recalculated and compared with output from our full simulation (i.e. spawning in all regions) and other reduced simulations.

To investigate whether specific oceanographic features were contributing to years of high and low predicted settlement in terms of recruitment, we plotted the paths of particles that successfully settled and particles that did not. By visually comparing the two groups of paths, we were able to make some broad conclusions of factors contributing to successful and unsuccessful settlement.

Overall and meridional velocity from BRAN2020 was moderately correlated with observations at most of the moorings used in the comparison (Fig. S4). BRAN2020 performed well for the EAC and Coffs Harbour mooring, with a correlation generally between 0.4 and 0.5, but poorly for the shallower waters of the Great Barrier Reef moorings, with a correlation of <0.3, although the magnitudes of the velocities were quite similar (Fig. S4). The zonal velocity in BRAN2020 was a much poorer match to the observations, with mostly low correlations (<0.2), except for the Coffs Harbour moorings, which were moderate (∼0.4). Temperature was well resolved by BRAN2020 with all correlations >0.7 and most >0.9. Overall, this suggests poor performance in the shallow areas with complex topography near the Great Barrier Reef to the north, and more reliable results in the EAC-influenced areas. The large-scale predominately along-shore patterns are well resolved (as already demonstrated by Chamberlain et al. 2021), whereas the small-scale patterns and cross-shelf transport are poorly resolved. BRAN2020 is therefore suitable for investigating broad-scale larval dispersal and connectivity (among zones), but is likely to not be suitable for small-scale movements such as within zones.

Larvae mostly settled south of their spawning location, although high levels of self-recruitment were evident within the northern zones. Zones 2 and 3 show relatively high amounts of self-recruitment and over 98% of larvae that settled in them originated within those two zones (Fig. 2a). By contrast, particles that settled further south originated from a larger area. Larvae settling in Zone 4 were primarily spawned in Zones 2, 3 and 4, whereas particles settling in Zones 5 through 7 were primarily spawned in Zones 3, 4 or 5 (Fig. 2a).

Fig. 2.  Connectivity matrices based on the particle tracking simulations. (a) The origins of particles (y-axis; source) that settle in each zone (x-axis; sink). (b) The settlement locations (x-axis) of particles spawned in each zone (y-axis). Numbers show the average percentage over the model duration, with standard deviation shown in parentheses. The referred to zones can be seen in Fig. 1. Combinations of zones where connectivity was never observed are blank. Note that the different-magnitude colour scales differ among subplots.

Taking a broader view and analysing the distribution of larvae as a proportion of total amount spawned in each zone, there is overall low survival (typical of natural mortality and dispersal), but a much higher settlement of particles spawned in the northern-most zones (2 and 3) than those spawned further south (Fig. 2b). Larvae spawned in Zones 4–7 were typically exported south from the inshore regions by the strong EAC, resulting in very low proportions of larvae from each of these zones settling because of a high dispersal mortality (Fig. 2b).

The Qld region (Zones 2–6) received ∼20× more larvae each year in our simulations than did NSW (Zone 7; Fig. 3). Within both regions, there were interannual variations; however, in terms of relative magnitude of the predicted settlement, the NSW region had larger fluctuations. In the NSW region, there was up to a 3× difference in predicted settlement among years. There were distinct periods of both high (1998–2000, 2008, 2012, 2013, 2015) and low (2001–2005, 2009–2011, 2014, 2018) settlement (Fig. 3).

Fig. 3.  Annual relative predicted settlement of spanner crab larvae in (a) Qld management region (Zones 2–6; Fig. 1) and, (b) the NSW management region (Zone 7; Fig. 1).

There was evidence of a positive correlation between the predicted recruitment index (based on the 4-year lagged, rolling mean of settlement) and both CPUE time series as well as the fishery-independent survey data (Fig. 4, S5). There was evidence of a moderate correlation between the long-term CPUE and our predicted recruitment index (r = 0.45, t₂₀ = 2.262, P = 0.035). Similarly, the short-term CPUE was positively correlated with the predicted recruitment index, but the evidence for the trend was weaker (r = 0.48, t₉ = 1.623, P = 0.139), which may be due to the timeseries being shorter (i.e. fewer observations and lower power). The strongest relationships were identified with the legal-size crab data from the fishery-independent survey (r = 0.73, t₁₂ = 3.793, P = 0.003). In the long-term series, there was an apparent breakdown in the correlations in the early 2000s, where predicted recruitment was high whereas CPUE was not.

Fig. 4.  Relationship between the NSW CPUE metrics and the predicted recruitment to the spanner crab fishery based on the larval settlement estimates. The black line shows the predicted recruitment to the fishery based on a 3-year rolling mean of predicted larval settlement and a 4-year lag. The red line shows the CPUE based on kilograms of harvest per fisher-day. The blue line shows the CPUE based on kilograms per net-lift. The green line shows the legal-size crab CPUE from the fishery-independent survey. Note that the y-axis was rescaled to assist interpretation.

The proportion of total harvest (Qld and NSW combined) taken in NSW showed strong evidence of being driven by recruitment. There was a moderate correlation with the predicted recruitment index (r = 0.51, t₂₁ = 2.262, P = 0.013). When the predicted recruitment to NSW was higher, NSW made up a higher proportion of the total harvest. Although the relationships appear to correlate better in early years than in later years, the relationships generally correlated well with the exceptions of 2012 and 2013 (Fig. 5, S6).

Fig. 5.  Relationship between the predicted recruitment index (black) and the proportion of total harvest taken in NSW (red; Zone 7, Fig. 1). The predicted recruitment index is based on a 3-year rolling mean of predicted larval settlement with a 4-year lag.

Overall patterns in settlement differed among months (Fig. 6). Settlement of November-spawned larvae was the most consistent, with high settlement predicted only in 1 high year (1993), contrasting with two peaks for December-spawned larvae (1998–2000, and 2012–2013). January-spawned larvae showed variation, but was overall more stable than was December-spawned larvae, with only one notable peak in predicted settlement in 2008.

Fig. 6.  Predicted settlement in NSW (Zone 7, Fig. 1) for larvae spawned in each month of our simulated spawning season: (a) November; (b) December; and (c) January. The y-axis is comparable among subplots, but is a relative number.

The pattern observed in the December-spawned larvae was noted to qualitatively match the observed patterns seen in both CPUE time series and fishery-independent survey data. After generating a predicted recruitment index to the fishery timeseries only on the basis of the December-spawned larvae, there was evidence of positive correlations with the CPUE time series (Fig. 7, S7). The long-term CPUE time series (kg fisher-day⁻¹) showed strong evidence of a moderate correlation with the predicted recruitment index (r = 0.49, t₂₀ = 2.481, P = 0.022). The shorter kilograms per net-lift time series showed a stronger positive correlation strength, but the evidence for the trend was only moderate because of the shorter time series (r = 0.59, t₉ = 2.168, P = 0.059). The fishery-independent survey data again showed the strongest relationships, with both all and legal-size crabs showing strong positive correlations (all crabs: r = 0.78, t₁₂ = 4.3378, P < 0.001; legal-size crabs: r = 0.84, t₁₂ = 5.447, P < 0.001). There was again a breakdown in the relationship between predicted recruitment index and CPUE in the early 2000s, seen in the longer time series.

Fig. 7.  Correlation between the two NSW CPUE metrics and the predicted recruitment to the spanner crab fishery based on December-spawned larvae settlement estimates. The black line shows the predicted recruitment index based on a 3-year rolling mean of predicted larval settlement and a 4-year lag. The red line shows the CPUE based on kilograms of harvest per fisher-day. The blue line shows the CPUE based on kilograms per net lift. The green line shows the relative CPUE for legal-size crabs from the fishery-independent survey. Note that the y-axis was rescaled to assist interpretation.

The sensitivity analysis testing the removal of spawning from each region showed that removing spawning entirely from any of the regions did not significantly alter the predicted recruitment to the NSW fishery. All variations were highly correlated with the equal spawning assumption (r ≥ 0.92) and showed the same peaks and troughs (Fig. S8), suggesting that the reported patterns are robust to spatial variations in spawning.

When comparing the paths of particles that successfully settled and those that did not, it was evident that excess southward transport was detrimental, as no particles returned north and successfully settled if they went south of 30°S (Fig. 8). At ∼25°S, Fraser Island (K’gari) extends into the Pacific Ocean, forming a ‘choke point’ where many larvae were retained to the north-west of the island, driving the high settlement rates in this zone (Zone 3), with very few of the particles being retained by the island, then leaving (Fig. 8a). In the northern fishing zones, there was a large persistent eddy that enabled some particles in 23–24°S to be transported off the shelf and still be returned (Fig. 8b). Further south, there were decreasing numbers of tracks that successfully returned to the shelf after being transported into deeper waters (Fig. 8b).

Fig. 8.  A random subset of tracks from larvae that (a) failed to settle in an appropriate area and, (b) successfully settled in an appropriate area. Successful settlement was defined as reaching the degree days threshold while on the continental shelf and within the fishing zones (2–7, Fig. 1).

This paper has demonstrated that oceanographic drivers are likely to influence dispersal and connectivity across the Australian spanner crab stock. The model suggests that the eastern Australian spanner crab stock is highly likely to be supported by larvae spawned within the northern half of the available fishing area, within the Qld jurisdiction. Modelling also indicated that larvae spawned within the southern portions of the Qld fishery and NSW are likely to be transported away from the suitable habitat areas (<200-m bathymetry within the fishing zones) by the East Australian Current, limiting the potential contribution of the southern zones to recruitment. Variation in CPUE in the southern fishing zone (both fishery-dependent and fishery-independent) correlated with predicted recruitment to the fishery solely on the basis of simulated larval dispersal. This indicated that fluctuations in the NSW spanner crab fishery are likely to be recruitment driven and that there may be potential to forecast abundance. Further, predicted recruitment can be used as a performance indicator for the stock and may support the development of decision rules within a harvest strategy.

There were distinct patterns in connectivity among spawning and settlement locations for spanner crabs. The two northern zones (Zones 2 and 3) showed high levels of self-recruitment, suggesting that these zones may be self-sustaining in terms of long-term larval supply. This contrasted with zones further south, with the probability of self-recruitment decreasing with increasing latitude. This is likely to arise because of the prevailing oceanography of the region, particularly the increasing influence of the EAC because its higher proximity to the southern spawning areas increases the likelihood of transporting the larvae south, away from possible settlement zones. The general observation for the southern zones is that almost all the larval supply comes from the north, transported by the EAC, an observation consistent with previous larval dispersal studies south of the current study (Roughan et al. 2011; Everett et al. 2017; Schilling et al. 2020).

A similar latitudinal gradient was seen in the differential ‘survival’ from the different zones. Successful settlement was an order of magnitude higher in the northern zones than in the southern zones, with successful settlement dropping from >5% of all spawned larvae in Zones 2 and 3 to <4% in Zone 4 and <1% in Zones 6 and 7. This is again likely to be due to the influence of the EAC, in addition to the fact that the area available for settlement in the southern zones is substantially smaller than in the northern zones, which have a wider continental shelf. The extension of Fraser Island out into the Pacific Ocean was also a significant feature that increased settlement by capturing larvae inshore, out of the EAC. The eastern side of Fraser Island is recognised as where the EAC is most coherent (Kerry et al. 2018), and larvae in this region will be transported southward rapidly. It is interesting to note that much of the stock is now located in the regions around Fraser Island, potentially contributing to the high recruitment in NSW via the dominant southward current.

We found moderate to strong correlations between the predicted recruitment and the observed variation in catch rates, depending on the metric used. The strongest relationship was found using December-simulated spawning and the legal-size crabs in the fishery-independent survey, aligning with the current belief that November–December is the peak spawning period (Skinner and Hill 1987; Kennelly and Watkins 1994). The predicted recruitment index from the larval dispersal simulations closely matched NSW CPUE in recent years, but the relationship broke down in the early 2000s. We believe this could be the effect of varying levels of spawning stock biomass resulting in an incorrect number of particles released in those years (e.g. perhaps the overall numbers spawned should be higher or lower in some years), although it is recognised that it is likely that abundances within the Qld fishery declined form the early 2000s to 2017 (Roelofs et al. 2021). A more likely alternative may be a shift in the dominant spawning biomass from Region 3 to Region 4, aligning with shifts in both the harvest and fishery-independent survey, which occurred around this time (Filar et al. 2021). This dynamic would not have been captured in our assumption of equal spawning in time and space; however, our sensitivity analysis showed that the overall patterns in the results would not have been altered.

When the settlement and survival patterns are viewed together, it strongly suggests that both the Qld (Zones 2–6) and NSW (Zone 7) zones are supplied by larvae from the Qld zones. As such, changes in larval supply from Qld (arising through a reduction in egg and larvae production or changes in oceanographic features that affect larval transport) may pose a risk to both NSW and Qld fisheries. The limited survival of larvae spawned within the NSW reporting zone suggests that the NSW fishery is heavily subsidised by the Qld spawning stock.

The dependency of the entire stock on larval supply from the Qld jurisdiction highlights the need to maintain a healthy spawning stock in this region. Any decline in spawning stock biomass and subsequent larval supply from Qld region, in particular from Zones 3–5, could have implications across the whole fishery. Going forward, the strong connectivity results here suggest that it may be wise to combine Qld and NSW management and assessment procedures.

A detailed study regarding the reproductive biology, including fecundity over the latitudinal range, would be advantageous to advance these particle-tracking modelling scenarios further beyond the current study. Recent research has also shown 4-year lagged effects of both chlorophyll-a (positive effect) and sea-level anomaly (negative effect) on catch rates in Region 3 (Filar et al. 2021); this suggests that there could be important influences beyond larval supply, such as the influence of food availability for larvae (Leggett and Deblois 1994), or mesoscale oceanographic effects such as eddies (Bakun 2006). Management of the Qld fishery is currently undertaken with a range of input and output controls, with the catch being limited by a total allowable catch (TAC). A target spawning biomass ratio of 60% has been imposed as part of a harvest strategy based on a maximum economic yield (Queensland Government 2019). If achieved, this 60% target should provide some security for both Qld and NSW fisheries in terms of continued larval supply.

An increasingly important factor in eastern Australia is climate change. Both northern and southern eastern Australia are being rapidly affected by warming waters and shifting currents. The Great Barrier Reef to the north is experiencing sustained heatwaves, resulting in increased bleaching and impaired coral recruitment (Hughes et al. 2017, 2019), occurring concurrently with a strengthening EAC (Sun et al. 2012; Li et al. 2021), resulting in increased water temperatures at higher latitudes, creating a global warming hotspot (Hobday and Pecl 2014). Located between these dramatic changes, the spanner crab habitat is warming, which may lead to a southward redistribution through time if the north of the range is unfavourably warm (Pecl et al. 2017; Gervais et al. 2021), although it is unknown whether the habitat further south would be suitable. It is likely that a climate-driven shift may have already occurred, with the catch in Region 3 declining, resulting in Region 4 becoming the most productive (Filar et al. 2021). If such a southward shift continues, larval supply could be detrimentally affected, with an increased southerly transport.

The interannual variations in predicted settlement in both Qld and NSW have implications for future CPUE-based stock assessments, and potential limitations on predicted biomass that is available for harvest. For example, models that operate on the assumption of average recruitment across all stock sizes, including Surplus Production Models (SPMs) previously applied to the Hawaiian spanner crab fishery (Kapur et al. 2019), are sensitive to potential changes and the lagged effects of recruitment (Punt and Szuwalski 2012). On the basis of predicted recruitment, estimates of maximum sustainable yield (MSY) and associated fisheries reference points derived from fitting SPMs to data from jurisdictions separately without adjustments based on predicted recruitment dynamics, may be biased.

Although our simulations showed strong agreement to patterns of catch, there are several factors potentially contributing to the disparity between predicted recruitment and fisheries harvest in some years. Some of these factors include (1) the 4-year time-lag between recruitment and harvest, (2) our estimation of larval mortality, (3) the disparity between the hydrodynamic model and reality, (4) the exclusion of swimming for our larvae, and (5) the uncertainty surrounding the biology of spanner crabs. Our predicted recruitment to the fishery was based on estimates that male spanner crabs enter the fishery at 4 years old (Brown et al. 2008). This is an assumption that carries considerable uncertainty, with previous minimum legal-size estimates ranging from 1.75 to 8.83 years for females and from 1.08 to 3.58 years for males (Kirkwood et al. 2005). It would be extremely valuable for further studies to better define and understand growth rates as well as the settlement of post-larvae (e.g. where 1–3-year-old juveniles primarily reside). In these 4 years, it is probable that other factors, such as environmental fluctuations, influenced the number of juveniles reaching maturity. Environmental factors such as temperature and current speed have been shown to influence the catchability of spanner crabs in this region (Spencer et al. 2017, 2019a, 2019b, 2019c), and it is possible that such factors may also influence juvenile survival. Environmental variables such as onshore winds have also shown correlations with larval fish supply in this region (Schilling et al. 2022). Swimming ability has also previously been demonstrated to significantly improve dispersal modelling of juvenile and larval fish, and turtles (Putman and Mansfield 2015; Leis 2021). Unfortunately, there is scant information regarding the swimming ability (speed or direction) of spanner crab larvae and it was, therefore, not possible to include it in our models. As the adult habitat of spanner crab larvae is the continental shelf of Australia, we suggest that swimming would generally occur in a westward direction and the overall connectivity patterns would remain similar because of the predominant flow of the currents along the coast, although overall recruitment might be higher.

There is considerable uncertainty around the reproductive biology of spanner crabs because of most of the reproductive information being over 20 years old. Predicted recruitment was simulated from different months as a control to test whether there have been possible shifts in spawning season. The collection of new data to describe the size–fecundity relationship at smaller spatial scales than previously investigated, and ongoing monitoring of the timing and duration of the spawning period in the Qld fishery, may provide a more accurate representation of the actual egg production to refine future simulations. Replicating existing annual standardised fishery-independent surveys during the identified spawning period would provide direct estimates of the abundance, size structure and the proportion of female crabs spawning. This would also allow for the collection of data to monitor spatial and temporal variations in the distribution and abundance of the spawning population, possibly even with reference to optimal sea-surface temperatures for spawning. This is particularly important, with rapid ocean warming occurring within the EAC region (Wu et al. 2012), and temperature clearly affecting fecundity in other decapods (Johnston and Yeoh 2021; Nolan et al. 2022). However, the cost–benefit of increased fishery-independent monitoring and the relative value of data on spawning populations (e.g. cost–benefit ratio on TAC) within an integrated fishery model for the East Coast biological stock would need to be determined to support this recommendation (Dennis et al. 2015).

Our assessment of the BRAN2020 model identified that it performed well in terms of temperature reconstruction and has previously been demonstrated to accurately represent large-scale dynamics (Chamberlain et al. 2021). However, BRAN2020 did not resolve small-scale dynamics, particularly cross-shelf velocities, adequately (Fig. S4). This may have resulted in the fine-scale movement dynamics of our larvae being poorly represented. To overcome this, it would be necessary to use a higher-resolution ocean model that better represents fine-scale dynamics on the continental shelf. Such models are available, but currently do not fully cover the area of interest with the main areas of interest in the models, either focusing north on the Great Barrier Reef (Steven et al. 2019; ∼9–28.5°S) or south on the Separation zone of the EAC (Kerry et al. 2018; ∼26–40°S). These models would need to have an extended domain or be nested within each other; both options are not currently available. Regardless, although the lack of fine-scale representation may mean that within-zone movement is poorly represented, we are confident that the large-scale connectivity findings are robust, with transport dominated by southerly currents. The strong correlations between lagged predicted recruitment and observed fluctuations in CPUE provide us with confidence in the general ability of our model and suggest that the dispersal is most influenced by larger mesoscale ocean dynamics.

The evidence presented here suggests that recruitment within the spanner crab stock in eastern Australia is likely to be heavily supported by spawning in the northern half of the fishing area. Larvae spawned from the northern and middle zones showed higher rates of successful settlement in all fishing areas than did larvae spawned in the more southern zones where the EAC has strong influence. In the NSW managed zone (Zone 7), we showed significant correlations between predicted recruitment to the fishery (using commercial and independent CPUE metrics) and modelled larval settlement. Combined with the connectivity results, this suggests that the NSW fishery is supported by recruits originating from the Qld fishery (Zones 2–6) and that joint management and assessments between the jurisdictions may be beneficial. By monitoring the arrival of larval or juvenile spanner crabs in NSW and using improved ocean models, it may be possible to reliably forecast future years of anomalously high and low CPUE in the NSW fishery.

Supplementary material is available online.

All modelling code is available in the github repository (see https://github.com/DEHewitt/portunid_particle_tracking). The NSW commercial data are available only upon request from NSW Department of Primary Industries Fisheries.

The authors declare that they have no conflicts of interest.

H. T. Schilling was supported by a NSW Government Research Attraction and Acceleration Program grant awarded to SIMS. D. E. Hewitt was supported by an Australian Government Research Training Program and a NSW DPIE top-up scholarship.