Culling corallivores improves short-term coral recovery under bleaching scenarios
Our model focused on the trophic interactions among CoTS and two groups of coral within a feedback loop with natural and anthropogenic forcing. Our model draws on accepted features of the published dynamics described by Morello et al.37, Condie et al.28 and Condie et al.17, but is a substantial advance in terms of adding spatial structure and coupling with climate variables. Here we have resolved a fine spatiotemporal model structure, developed a novel recruitment formulation for CoTS, integrated tactical management control dynamics and incorporated the impact of broad-scale drivers upon the population dynamics of corals and CoTS at the local scale. Our model is formally fitted to a subset of the CoTS control program data described by Westcott et al.12. We operationalised our model as a tactical and strategic tool to inform how CoTS management strategies interact with alternative disturbance and ecological realisations at the sub-reef scale, the scale at which management operates.DataWe fitted our model to a subset of four reefs from the dataset described by Westcott et al.12, which were consistently and intensively managed (for a map with reef locations see Fig. 2). We restricted our focus to a subset to avoid parametrisation of reef and management site dynamics. Thus, ~39% of site visits were concentrated over the 13 management sites we considered, with a mean of 20.73 ± 5.5 (mean ± standard deviation) visits across the time series relative to a mean visitation rate of 12.23 ± 4.7 (mean ± standard deviation) for the rest of the sites. Each reef in the subset contained two or more management sites where each site was visited at least 18 times. The subset was used because it contained sufficient data for estimating the 11 model parameters for each management site. Across included sites were a range of CoTS densities, coral abundances and disturbance histories12,72,73. Given the intensity with which these sites were managed, they therefore provided us with a valuable opportunity to formally fit the interactions between management intervention, coral abundance and CoTS dynamics in the presence of regional sequential bleaching events.Model spatial structure and ecological componentsSpatially, we considered a circular 300 km region of the Great Barrier Reef centred between Cairns and Cape Tribulation, and resolved at a daily timescale and a sub-reef spatial scale, matching the scale at which observed data were resolved12,19. Reefs were randomly generated as points to capture possible spatial correlation in disturbance impacts between nearby reefs, as well as to allow variability in reef locations. Coral, CoTS and disturbance dynamics within the management sites of each reef were resolved relative to a 1 ha focal region. That is, each management site was captured as a 1 ha area representative of the whole site. In the Pacific, Acanthaster spp. disproportionally target faster-growing corals, predominantly Acropora, Pocillopora and Montipora22. Coral taxa characterised by slow growth rates and massive morphologies, such as Porites, are generally consumed less than expected based on their abundance22 and are thus non-preferred prey. The two modelled coral groups were the fast-growing favoured prey items of CoTS, and the slower-growing non-preferred prey. Processes resolved in the model included reproduction, density dependence, the effect of bleaching and cyclonic disturbances on corals and the impact of manual control (culling) upon CoTS and coral dynamics.CoTS population structureWe used an age-structured approach to model CoTS population dynamics. We defined our age classes to encapsulate plausible size-at-age variation due to plastic growth. This was achieved through linking catch size classes of the management control program19 to age classes through size-age relationships developed from observations spanning multiple environmental realisations, manipulated scenarios and methodologies55,70,74,75. Delayed growth in juvenile CoTS due to deferral of their switch to coral prey or composition of their pre-coral diet, may induce variability in the size-at-age of juveniles52,53. However, the population-level consequences of prolonged juvenile phases are not easily observed nor understood. For example, juveniles are subject to high mortality rates in situ, delayed growth may reduce lifetime fitness and there have been no observations of juveniles during spawning periods that would indicate protracted juvenile phases55,56,57. Consequently, suggests size-at-age is—due to an early life history mortality bottleneck or otherwise—predominantly concordant with growth curves of the literature55,70,74,75 and the size classes we have used here. Age classes comprised annual 0, 1, 2 and 3+ groups, with 3+ being an absorbing class – once there, they stay there. Age-0 ( ; 32.5)). This induced a slope change in the relationship between maximum wind velocity and its radius at a wind velocity of 32.5 m.s−1 (≥ category 3 intensities). However, whilst maximum wind velocity was modelled to determine ({d}_{{{{{{rm{m}}}}}}}), the overall size of the cyclone was uncorrelated with its intensity. The overall size was uniformly sampled from 130 to 460 km diameter which allowed for the potential of complete focal area coverage and for a range of intensity-size relationships to be captured. Given a cyclone footprint of radius ({d}_{0}) (km), wind velocity, (V) (m.s−1), at a distance, (d) (km), was interpolated104 through:$$Vleft(dright)=left{begin{array}{c}{V}_{0}+left({V}_{{{{{{rm{m}}}}}}}-{V}_{0}right){left(frac{sqrt{{d}_{0}}-sqrt{d}}{sqrt{{d}_{0}}-sqrt{{d}_{{{{{{rm{m}}}}}}}}}right)}^{alpha },,dge {d}_{{{{{{rm{m}}}}}}}\ {V}_{{{{{{rm{m}}}}}}}, , d ; < ; {d}_{{{{{{rm{m}}}}}}}end{array}right.$$ (32) The distance from the cyclone centre to the reef perimeter, (D) (km), is calculated through:$$D=sqrt{{left({x}_{{{{{{rm{rf}}}}}}}-{r}_{1}-{x}_{{{{{{rm{cyc}}}}}}}right)}^{2}+{left({x}_{{{rm{rf}}}}-{r}_{1}-{y}_{{{{{{rm{cyc}}}}}}}right)}^{2}}$$ (33) Thus, given a reef strike occurs ((sqrt{{d}_{0}}-sqrt{d}ge 0) required from non-integer (alpha)), the wind velocity experienced at said reef due to the tropical cyclone was calculated as (Vleft(Dright)). Wind velocity was subsequently categorised and damage to reef zone corals calculated as per Supplementary Table 4.We resolved stochasticity in cyclone dynamics in projection scenarios. In projected scenarios cyclone arrivals, locations and intensities were probabilistically sampled and their inflicted damage upon coral communities sampled from damage ranges. Cyclone locations, their footprints, intensity ranges and corresponding damage ranges were sampled from uniform distributions. Cyclone arrivals were sampled from a Poisson distribution and considered in scenarios from 2018 to 2029. Projections were averaged over 80 simulations to capture mean dynamics and bound trajectory uncertainty due to said stochasticity.Our cyclone model was calibrated to parameters sourced from the literature (Supplementary Tables 4-5). This was necessary since our data time series did not encompass a cyclone event and/or impacts upon a reef and cyclone-induced mortality is typically a key coral mortality source30. Consequently, we were unable to validate the impacts of cyclones through formal estimation in our model. However, our endeavours to source parameters from empirical and modelling studies in conjunction with our formulation allowed us to plausibly capture the cumulative outcomes of a cyclone event at discrete locations. Our cyclone model offers a limited complexity approach that is empirically grounded to simply resolve cyclone impacts in local-scale models without the need to be coupled to a regional-scale model.Cyclones, induced thermal stress and tactical managementThe occurrence of cyclone events was modelled to directly interact with both management interventions and thermal stress events. Cyclones were assumed to realistically preclude co-occurring co-located management interventions. This was such that a management site control visit was abandoned if a cyclone preceded or was forecast within five days of a control voyage. The later interaction of cyclones with thermal stress events operated through an induced thermal cooling of sea surface temperatures (SST) at impacted locations.In the case of the overlapping cyclone and thermally induced bleaching events, we first accounted for cyclone impacts. This was because, in addition to physical damage to corals, cyclones have the potential for regional-scale cooling of SST which can reduce coral bleaching43,107. To capture this interaction, we resolved the duration108,109 and amplitude107 of tropical cyclone-induced cooling. We captured this interaction through Degree Heating Weeks (DHW) which is a useful metric for the accumulated thermal stress experienced by corals94.The duration of tropical cyclone-induced cooling was modelled through a temporal-SST response curve consistent with the work of Lloyd and Vecchi108 and Vincent et al.109. Cooling rapidly occurs once a tropical cyclone arrives at a location and decays in an asymptotic manner over a period of ~40–60 days108,109. Temperatures however do not return to pre-cyclone levels and plateau at ~1/4 of the cooling signal amplitude below pre-cyclone levels108,109. We expressed this cooling response curve as it related to bleaching-induced coral mortality through DHWs.We based the average expected DHW cooling signal on the work of Carrigan and Puotinen107. This was achieved through scaling the difference in amplitude of overlapping thermal stress-tropical cyclone events and thermal stress only events—a cooling signal amplitude of ({{{{{{rm{DHW}}}}}}}_{{{{{{rm{Amp}}}}}}} sim 1.5) DHW. Consistent with the model of Carrigan and Puotinen107, we then resolved cooling within the radius of gale-force winds (category 1, 17 m.s−1) to model tropical cyclone-induced cooling. Depending on the size of the tropical cyclone, this meant that an individual cyclone would not necessarily cool all reefs within the model region. However, the culmination of multiple cyclones may have limited bleaching exposure for corals across the region107.We did not treat the cooling consequences of multiple cyclones additively nor the complex interplay of oceanic feedbacks upon cyclone intensity and cooling. Such processes were beyond the scope of our study and model. If multiple cyclones occurred within our model, then the cooling signal timeline was re-initialised at impacted reefs for the last tropical cyclone at said location. Non-impacted reefs maintained the timeline for the decay of the cooling signal originating from their previous tropical cyclone interaction.Once a tropical cyclone impacted a reef, the duration of the induced cooling signal was modelled. Price et al.110 found that cooling decays exponentially which is reflective of the recovery of SST following tropical cyclones as demonstrated by Lloyd and Vecchi108 and Vincent et al.109. We operationalised the exponential functional form in conjunction with the decay timelines of Lloyd and Vecchi108 and Vincent et al.109 and the DHW amplitude of Carrigan and Puotinen107. We modelled the level of cooling ({{{{{{rm{DHW}}}}}}}_{{{{{{rm{cool}}}}}}}) after ({d}_{{{{{{rm{postTC}}}}}}}) days post-cyclone event by:$${{{{{{rm{DHW}}}}}}}_{{{{{{rm{cool}}}}}}}left({d}_{{{{{{rm{postTC}}}}}}}right)=frac{1}{4}{{{{{{rm{DHW}}}}}}}_{{{{{{rm{Amp}}}}}}}+frac{frac{3}{4}{{{{{{rm{DHW}}}}}}}_{{{{{{rm{Amp}}}}}}}}{{e}^{{d}_{{{{{{rm{postTC}}}}}}}/10}}$$ (34) This ensured that once a reef experienced a tropical cyclone event, the cooling signal initialised at ({{{{{{rm{DHW}}}}}}}_{{{{{{rm{Amp}}}}}}}) and decayed to (sim frac{1}{4}{{{{{{rm{DHW}}}}}}}_{{{{{{rm{Amp}}}}}}}) after 40–60 days108,109. The rate of decay was given by the e-folding time (days required for the cooling signal to be reduced by a factor of (e)) which we took to be 10. This is consistent with the results of Price et al.110, Lloyd and Vecchi108 and Vincent et al.109 who found e-folding times ranging from 5 through to 20 days. Thermally induced bleaching mortality of corals was computed after cyclone physical damage and cooling had been accounted for.Formal model fittingWe formally fitted our coral-CoTS model simultaneously to coral cover data, catch-per-unit-effort data and catch numbers obtained from the management control program with dive effort (minutes) treated as an input (visits summarised in Supplementary Table 7)12. Simultaneously fitting CoTS and coral dynamics at concurrent locations was useful here as it allowed for coral cover trajectories to help inform local CoTS abundance (sensu CoTS feeding vs. coral trajectories63,79 and local site fidelity24). Our model also used Long Term Monitoring Program (LTMP) data (based on manta tows and provided by the Australian Institute of Marine Science) which provides an independent index of relative abundance of CoTS. This was such that our model here was developed and parametrised based on an earlier version37,111 which did not use CPUE information but was fitted to the LTMP data on CoTS relative abundance, as well as the corresponding coral cover, to estimate a number of CoTS-coral interaction parameters used in the present model (Supplementary Table 3).Fitting and estimation of our model were achieved through Maximum Likelihood Estimation (MLE). Our objective function was the outcome of combining the negative log-likelihood contributions arising from fitting the model to multiple sets of location-specific data, across a range of environmental and ecological realisations, in conjunction with penalty terms. Specifically, we fitted coral cover (data series ({x}^{{{{{{rm{Coral}}}}}}})) and CoTS CPUEs (data series ({x}^{{{{{{rm{CoTS}}}}}}})) at each management site which contained ({n}_{{{{{{rm{Coral}}}}}}}) and ({n}_{{{{{{rm{CoTS}}}}}}}) data points respectively. This involved fitting parameters that were specific to management sites (e.g. thermal stress - DHW), reefs (e.g. recruitment variability) as well as those that were common amongst reefs (e.g. CoTS consumption rates). A parametrisation that optimised one contribution was unlikely to optimise all contributions and hence we obtained a parametrisation across all reefs and sub-regions. For a modelled catch of (N) (sum of catches across age classes), a catchability coefficient (a constant of proportionality) of ({q}_{{{{{{rm{LL}}}}}}}^{{{{{{rm{prop}}}}}}}), and data standard deviation of ({sigma }_{{{{{{rm{LL}}}}}}}) our likelihood contribution arising from a management site CPUEs was given by:$$-{{log }}{{{{{rm{L}}}}}}left({q}_{{{{{{rm{LL}}}}}}}^{{{{{{rm{prop}}}}}}}N,{{sigma }_{{{{{{rm{LL}}}}}}}}^{2}{{{{{rm{|}}}}}}{x}_{i}^{{{{{{rm{CoTS}}}}}}}right) = {n}_{{{{{{rm{CoTS}}}}}}},{{{{{rm{ln}}}}}}left({sigma }_{{{{{{rm{LL}}}}}}}right)+{sum }_{i=1}^{{n}_{{{{{{rm{CoTS}}}}}}}}frac{{left({{{{{rm{ln}}}}}}left({x}_{i}^{{{{{{rm{CoTS}}}}}}}right)-{{{{{rm{ln}}}}}}left({q}_{{{{{{rm{LL}}}}}}}^{{{{{{rm{prop}}}}}}}{N}_{i}right)right)}^{2}}{2{{sigma }_{{{{{{rm{LL}}}}}}}}^{2}}$$ (35) From which the data series variance and catchability coefficient were computed for the maximum likelihood estimate. The derived variance and the catchability were respectively computed as per:$${sigma }_{{{{{{rm{LL}}}}}}}=sqrt{frac{1}{{n}_{{{{{{rm{CoTS}}}}}}}}{sum }_{i=1}^{{n}_{{{{{{rm{CoTS}}}}}}}}{left({{{{{rm{ln}}}}}}left({x}_{i}^{{{{{{rm{CoTS}}}}}}}right)-{{{{{rm{ln}}}}}}left({q}_{{{{{{rm{LL}}}}}}}^{{{{{{rm{prop}}}}}}}right)right)}^{2}}$$ (36) and$${q}_{{{{{{rm{LL}}}}}}}^{{{{{{rm{prop}}}}}}}=frac{1}{{n}_{{{{{{rm{CoTS}}}}}}}}{sum }_{i=1}^{{n}_{{{{{{rm{CoTS}}}}}}}}left({{{{{rm{ln}}}}}}left({x}_{i}^{{{{{{rm{CoTS}}}}}}}right)-{{{{{rm{ln}}}}}}left({N}_{i}right)right)$$ (37) Similarly, the likelihood contribution arising from fitting to a management site coral cover with standard deviation ({sigma }_{{Coral}}) was described by:$$-{{log }}{{{{{rm{L}}}}}}left(frac{{C}_{y,d}^{{{{{{rm{f}}}}}}}+{C}_{y,d}^{{{{{{rm{s}}}}}}}}{{K}^{{{{{{rm{coral}}}}}}}},{{sigma }_{{{{{{rm{Coral}}}}}}}}^{2}{{{{{rm{|}}}}}}{x}_{i}^{{{{{{rm{Coral}}}}}}}right) = {n}_{{{{{{rm{Coral}}}}}}},{{{{{rm{ln}}}}}}left({sigma }_{{{{{{rm{Coral}}}}}}}right)+{sum }_{i=1}^{{n}_{{{{{{rm{Coral}}}}}}}}frac{{left({ln}left({x}_{i}^{{{{{{rm{Coral}}}}}}}right)-left(frac{{C}_{y,d}^{{{{{{rm{f}}}}}},i}+{C}_{y,d}^{{{{{{rm{s}}}}}},i}}{{K}^{{{{{{rm{coral}}}}}}}}right)right)}^{2}}{2{{sigma }_{{{{{{rm{Coral}}}}}}}}^{2}}$$ (38) Where the standard deviation was given by:$${sigma }_{{{{{{rm{Coral}}}}}}}=sqrt{frac{1}{{n}_{{{{{{rm{Coral}}}}}}}}{sum }_{i=1}^{{n}_{{{{{{rm{Coral}}}}}}}}{left({{{{{rm{ln}}}}}}left({x}_{i}^{{{{{{rm{Coral}}}}}}}right)-{{{{{rm{ln}}}}}}left(frac{{C}_{y,d}^{{{{{{rm{f}}}}}},i}+{C}_{y,d}^{{{{{{rm{s}}}}}},i}}{{K}^{{{{{{rm{coral}}}}}}}}right)right)}^{2}}$$ (39) We computed the negative log-likelihood objective function by summing the contributions from all management sites across considered reefs.Fitting was conducted through the modelling language Automatic Differentiation Model Builder (ADMB) which implements a Quasi-Newton optimisation algorithm for estimation of parameters and provides Hessian based estimation of standard errors112. Penalty terms were added to our likelihood function to integrate a prior understanding of system dynamics and to reduce model variability. Penalty terms encompassed recruitment variability and the magnitude of catches observed in the data.Recruitment was expressed through recruitment deviations, ({r}_{y}), given a standard deviation of ({sigma }_{{{{{{rm{R}}}}}}}) about underlying modelled recruitment (sum of self-recruitment and immigration sources described previously). The recruitment variability negative log-likelihood penalty contribution was given by:$$-{{log }}{{{{{rm{L}}}}}}left(0,{sigma }_{{{{{{rm{R}}}}}}}^{2}{{{{{rm{|}}}}}}{r}^{{{{{{rm{rec}}}}}}}right)={sum }_{y=1}^{{{{{{rm{#Years}}}}}}}{sum }_{{{{{{rm{reef}}}}}}=1}^{{{{{{rm{#Reefs}}}}}}}{r}_{y,{{{{{rm{reef}}}}}}}^{{rec}}/2{sigma }_{{{{{{rm{R}}}}}}}^{2}$$ (40) An additional penalty term for model deviations from the magnitude of observed catches was encompassed. This was such that a constant of proportionality relating modelled catches to observed catches tended to one. For an allowed standard deviation of ({sigma }_{{{{{{rm{CM}}}}}}}), the likelihood function was penalised for deviations from unity proportionality, ({r}^{{{{{{rm{CM}}}}}}}), through:$$-{{log }}{{{{{rm{L}}}}}}left(0,{sigma }_{{{{{{rm{CM}}}}}}}^{2}{{{{{rm{|}}}}}}{r}^{{{{{{rm{CM}}}}}}}right)={sum }_{{{{{{rm{zone}}}}}}=1}^{{{{{{rm{#Zones}}}}}}}{r}_{{{{{{rm{zone}}}}}}}^{{{{{{rm{CM}}}}}}}/2{sigma }_{{{{{{rm{CM}}}}}}}^{2}$$ (41) Model simulations were conducted in ADMB with output analysis and visualisation conducted in MATLAB.Sensitivity to CoTS controlTo test whether our projected scenarios were consistent with the period over which data were collected, we conducted a model-based before and after comparison to the impact of control. Specifically, we used the fitted trajectory for sites, including both the coral data and CoTS control data (voyages and time spent), and compared this to the model-suggested coral trajectories if CoTS control had not taken place. These were modelled over the fitted period (2013–2018) and, unlike the projected scenarios (2019–2029), were variable in terms of the timing of control (amount of time between visits was variable), the amount of time spent at sites (not a consistent number of dive minutes per visit), CoTS dynamics (recruitment was fitted and hence different annually and between reefs), and in the level of thermal stress they experienced (different sites experienced different effective levels and some sites experience back-to-back events).Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More
