Interspecific Interactions May Influence Reef Fish Management Strategies in the Gulf of Mexico

This study highlights the importance of interspecific interactions among marine organisms and the effect that these trophic interactions have on the development of effective, adaptive management strategies for reef fishes in the Gulf of Mexico. To represent the spatially and temporally constrained, interspecific interactions among reef fishes we employ Atlantis (a spatially explicit, biogeochemical ecosystem model) as our simulation tool. Within Atlantis, we evaluate the performance of a two-point harvest control rule (HCR) that adaptively increases fishing mortality linearly between upper and lower biomass thresholds based on the available biomass of the stocks. This example demonstrated the use of a “blanket” two-point HCR that assessed the available biomass of several reef fish species (often co-caught in fishing gear) both simultaneously and objectively. To estimate the impact of reef fish fishing on species abundance and biodiversity in the ecosystem, we examined four “low” and four “high” fishing mortality (F) scaler scenarios. All model projections are forward looking, representing a 50-year time horizon (2010 to 2060). We evaluated the performance of the two-point HCRs under the eight fishing mortality scenarios using ecosystem metrics that were previously found to robustly track changes in ecosystem function caused by fishing. We found that the lower F scenarios produced an ecologically distinct ecosystem state compared with the higher F scenarios, where relatively higher levels of fishing mortality (particularly on predators such as the deep Serranidae group) resulted in an increase in prey availability in later years of the simulation. This led to an increase in the overall productivity of the ecosystem over time and higher catch and biomass of most other reef fish groups at equilibrium (year 50). Our results suggest that a better understanding of interspecific interactions among targeted reef fishes and their prey is critical to developing ecosystem-based management strategies for the Gulf of Mexico. Subject editor: Kenneth Rose, University of Maryland Center for Environmental Science, Cambridge *Corresponding author: mbonewit@mail.usf.edu Received February 6, 2017; accepted October 26, 2017 This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. Marine and Coastal Fisheries: Dynamics, Management, and Ecosystem Science 10:24–39, 2018

Single-species assessments may not adequately capture uncertainty when strong interspecific interactions exist between targeted species, such as in a reef fish assemblage. These shortcomings may be significant impediments to effective management of depleted and recovering stocks. In the Gulf of Mexico, four reef fish stocks have been considered overfished in recent years: Gag Mycteroperca microlepis, Greater Amberjack Seriola dumerili, Gray Triggerfish Balistes capriscus, and Red Snapper Lutjanus campechanus (GMFMC 2013). All four of these stocks have been under stock-rebuilding plans developed by the Gulf of Mexico Fishery Management Council (GMFMC). Since 1984, there have been multiple changes to the original reef fish fishery management policies (FMPs) in the Gulf of Mexico, including the establishment of individual fishing quotas (IFQs), changes to gear restrictions and size limits, and the implementation of total allowable catch (TAC) limits. However, there has been little advancement in improving our understanding of the role that interspecific interactions have on effectively managing reef fishes in the Gulf of Mexico (Okey et al. 2004;Walters et al. 2008;Chagaris et al. 2015).
Fisheries managers have the potential to enhance reef fish management in the Gulf of Mexico using strategic, whole-system, simulation tools that can improve our understanding of complex ecosystem processes (Link 2010;Gr€ uss et al. 2017). To date, there are several ecosystem simulators that can be used to address a wide range of conceptual, strategic, and tactical hypotheses (e.g., Walters et al. 1997;Fulton et al. 2004b;Kazanci 2007). For example, in the Gulf of Mexico, Ecopath with Ecosim, OSMOSE, and Atlantis models are being used to incorporate multispecies considerations into the management decision process Gray et al. 2013;Gr€ uss et al. 2013Gr€ uss et al. , 2017. The utility of full system models in an ecosystem-based fisheries management (EBFM) context is to represent an extensive suite of ecosystem processes that can affect target species as well as nontargeted (or less valuable) species (Link 2010). Representing these integrated ecosystem processes is essential to EBFM, as these dynamics can strongly influence fisheries productivity and safe harvest rates.
Evaluating the performance of management strategies is particularly useful in a tool like Atlantis that is spatially explicit in three dimensions, has dynamic fleet behavior, can account for the spatially and temporally dynamic, interspecific interactions among modeled species, and produces quantifiable outputs that allow for the explicit evaluation of various harvesting policies (Fulton et al. 2004b). Thus, simulation results provide strategic guidance for managers beyond what is available from a single-species approach. For example, recent studies used Atlantis modeling software (Fulton et al. 2004b) to account for trophic and environmental effects on productivity (Fulton et al. 2007Smith et al. 2011Smith et al. , 2015. In this study, we evaluated the performance of a twopoint harvest control rule (HCR), as two-point HCRs were found to be robust to error in other systems (e.g., Parma 2002;Deroba and Bence 2008;Punt et al. 2008). We examined four "low" fishing mortality (F) and four "high" F scaler scenarios under a "blanket" two-point HCR, where the same harvesting policy is applied to a complex of harvested reef fishes each year. The reef fishes (or model functional groups) evaluated in this study were Gag, Red Grouper Epinephelus morio, deep Serranidae spp., Red Snapper, Vermilion Snapper Rhomboplites aurorubens, and Lutjanidae spp., not elsewhere included) (Table 1). We refer to these six Atlantis model groups throughout the rest of this paper as the "reef fish complex." Since population dynamics within the reef fish complex are closely linked by competition and predation (Masi et al. 2014;Chagaris et al. 2015), we conducted a sensitivity analysis of the initial diet (or availability) matrix on a subset of simulations to bracket the uncertainty in model outputs. The solutions from these scenarios were used to develop equilibrium yield curves and to derive ecological indicators that were found to be robust and important in previous analyses (Masi et al. 2016). Using these policy performance metrics, we quantified the ecosystem-level tradeoffs among the various levels of fishing mortality. The results from this analysis show the potential benefits of a simple, adaptive management policy that can be applied across a range of co-caught species and highlight the importance of accounting for the interspecific interactions among reef fishes in developing adaptive management policies in the Gulf of Mexico.

METHODS
Atlantis ecosystem model of the Gulf of Mexico.-Atlantis (http://atlantis.cmar.csiro.au/) is a three-dimensional, spatially explicit, trophodynamic ecosystem model that integrates biology, physics, chemistry, and human impacts (Fulton 2001;Fulton et al. 2004aFulton et al. , 2004bFulton et al. , 2005Fulton et al. , 2007Fulton et al. , 2011. In this application, we used an Atlantis model to represent the Gulf of Mexico marine ecosystem. The parameterization and calibration of the 2010 Gulf of Mexico Atlantis model used in this study is described in Ainsworth et al. (2015), in which Drexler and Ainsworth (2013) set initial biomass distributions, and Masi et al. (2014) and Tarnecki et al. (2016) developed the diet matrix. Only features of the Atlantis model framework pertinent to this analysis are reviewed here.
A key strength to Atlantis is the use of polygons , which spatially delineates bioregions, management jurisdictions, social or industrial structures, and, of particular value to this study, the spatial and temporal separation of predators from their prey and fishing fleets to their catch. The polygon structure designed for the Gulf of Mexico Atlantis model includes 66 polygons (Figure 1). Each polygon has associated weightings, which represent the prevalence of certain physical habitat types. The prevalence of biogenic habitat types is linked to the biomass of habitat-forming functional groups (e.g., seagrass). Habitat availability in the model affects the distribution of functional groups during the dynamic simulations according to a habitat affinity matrix.
The Atlantis modeling framework also consists of submodels that represent dynamic biochemical processes, exploitation, and the formal assessment of selected groups (Fulton et al. 2004b;Link et al. 2010). In the Gulf of Mexico Atlantis biological submodel, there are 91 model functional groups (Table A.1 in the Appendix). Each group consists of either individual species (e.g., Gag) or aggregated groups of species that share similar diets, habits, or niches. Vertebrate functional groups are tracked by numbers of individuals and mean body weight per individual, while invertebrate groups are treated as a single biomass pool. Changes in biomass for vertebrate or invertebrate consumers are tracked (at 12-h time steps) according to equation (1), where biomass (B) is substituted by abundance per age-class in the case of vertebrate consumers: where G is population growth (in biomass per unit time) (equation 2), M i is mortality due to predator i, n is number of predators, M is natural mortality not captured by trophodynamics (equation 3), I is immigration into model domain, and F is fishing mortality. Population growth is expressed as, where P i is predation by consumer on prey i, ε i is assimilation efficiency on prey i, δ O2 is an oxygen limitation factor, δ space is a space limitation factor, and A is the rate of catabolism. Natural mortality (not represented by predation mortality) is given as, In the Gulf of Mexico Atlantis model, biomass for group i is only constrained by density-independent linear 26 mortality (M lin ), as quadratic mortality (M quad ) and special mortality (e.g., to represent mechanical stress on macroalgae) were set to zero. Interspecific interactions are assigned in the initial conditions file using an availabilities matrix. Predation (P) by vertebrate functional groups on available prey i by predator j (in biomass per unit time) follows equation (4), which is the modified version of the Holling Type II functional response (Holling 1959): where B i is biomass of prey I, B j is the biomass of predator j, C j is the clearance rate (i.e., grazing efficiency) of predator j, G j is the growth rate of predator j, E ij is the growth efficiency of predator j eating prey I, and a ij is the availability of prey i to predator j. Group-specific parameter inputs are presented in Ainsworth et al. (2015). The food web interactions among the six reef fish complex groups and their top five prey items (both adults and juveniles) at model start are shown (Figure 2), though model dynamics allow interspecific interactions to vary among runs (e.g., based on prey availability). Here, we also included the top prey item of the reef fish complex's prey to emphasize the role that dynamic food web interactions can have on ecosystem structure and function ( Figure 2). Feeding rates can also vary dynamically according to specified gape limitations, which directs predation mortality to prey groups and age-classes that fall within accessible size ranges. Gape limitations are represented as a fraction of a predator's body weight. Further explanation of feeding dynamics is available in Fulton et al. (2004b).
In the Gulf of Mexico Atlantis model, recruitment is assumed to follow a Beverton-Holt recruitment relationship for all vertebrate groups except for sea birds, marine mammals, and sea turtles; this assumes a fixed number of offspring per female per year (Ainsworth et al. 2015). Along with initial estimates of growth and availabilities (a ij ), recruitment parameters for all groups that use Beverton-Holt recruitment were adjusted iteratively in tuning until observed biomass trends matched the predicted (Ainsworth et al. 2015). In our implementation of Atlantis, new recruits were assigned to the first of 10 age-classes and entered at a group-specified day after spawning (Ainsworth et al. 2015). The recruitment window varies among functional groups, and reproduction is dynamically dependent on both maternal condition during the spawning window and densitydependent factors (e.g., available habitat). Maternal condition was assessed as the ratio of reserve (soft tissue) nitrogen to structural (hard tissue) nitrogen. Another key assumption in the biology submodel is the density-dependent movement of predator functional groups toward areas with higher prey availability and the seasonal and annual migration of species both among polygons ( Figure 1) and into and out of the model domain. In the Gulf of Mexico Atlantis model, the seasonal movement patterns of each group were set according to the vertebrate and invertebrate concentrations per polygon determined by generalized additive models (GAMs)  or by expert opinion for highly migratory pelagics (E. Orbesen, National Oceanic and Atmospheric Administration [NOAA], Southeast Fisheries Science Center, personal communication).
In the Gulf of Mexico Atlantis exploitation submodel, fishing fleets were assigned based on gear type, targeted species, and selectivity patterns. Although sometimes only as bycatch, the reef fish complex groups are co-caught by several of the Gulf of Mexico Atlantis model's fishing fleets (Table A.2), which holistically represents the way in which reef fishes are co-caught by fleets operating in the Gulf of Mexico (Saul and Die 2016). There are roughly 60 marine protected areas (MPAs) included in the Gulf of Mexico Atlantis model, which spatially and temporally limit the fishing fleets to the selected fish groups during a simulation.
The exploitation submodel supplies the simulated data during a run to the assessment submodel, which houses an integrated assessment routine-the management strategy evaluation (MSE) routine. The MSE routine is designed to simulate a "closed-loop" management decision making process. It relies on an HCR to adjust F each year dynamically, based on the available biomass of an assessed model functional group (or groups). The MSE routine can be set up in Atlantis to assess a variety of harvesting options. However, in this study we were interested in the applicability of establishing a blanket two-point HCR to manage co-caught reef fishes in the Gulf of Mexico. Like a feedback control, a two-point HCR works by linking a control variable (F) to a state variable (e.g., total annual biomass) (Roel and DeOliveira 2007;Froese et al. 2011;Little et al. 2011;Eikeset et al. 2013). Similar to precautionary, single-species, control rules that are used to manage reef fish in the Gulf of Mexico (GMFMC 2017), the two-point HCR applied in Atlantis follows a "brokenstick" shape ( Figure 3), and requires the prescription of an upper biomass threshold, a lower biomass threshold, and a maximum fishing mortality rate. The maximum F allowed by the two-point HCR is defined as F mult . We treated this as a unitless scaler proportionate to the 2010 F rates estimated in Ainsworth et al. (2015) (e.g., F mult × 0.5 refers to one-half of the 2010 F). In our application, the biomass thresholds (lower: B low ; upper: B up ) are based on a proportion of the initial (2010) biomass (Ainsworth et al. 2015).
Each year, the available biomass is passed (internally, within Atlantis) to the MSE routine for each functional group in the reef fish complex. We assumed perfect knowledge in order to characterize the potential benefits of the two-point HCR. If the current biomass (B) is greater than B up , the maximum allowable fishing mortality is applied on that group, F mult (Figure 3). When B is below B low , a fishing mortality rate of zero is applied. When B is between the upper and lower thresholds, the fishing rate for the next year (F Applied ) is determined as in equation (5): In this application, we ran 24 two-point HCR simulations in Atlantis, where we varied F mult and B up . In total, we evaluate three variants on the upper biomass threshold of the two-point HCR (40%, 60%, and 80% of the 2010 FIGURE 3. Typical, two-point (hockey stick) harvest control rule (HCR) relating fishing mortality (F) to biomass (B) thresholds. 28 biomass by group). Within each of those variants, we evaluated eight variants on the F mult , looking at four low F scaler scenarios (F mult × 0.5, 0.8, 1, and 1.5) and four high F scaler scenarios (F mult × 3, 12, 25, and 50). We chose moderately high F scalers in this analysis, but the final realized fishing mortality was a function of this scaler, the spatial and seasonal overlap of each fleet (Table A.2), and the available biomass of the assessed groups throughout the simulation period. The F mult was applied simultaneously (i.e., as a "blanket" policy) to all six reef fish groups in the complex (Table 1). In all 24 Atlantis simulations, the lower biomass threshold was held constant at 20% of the 2010 biomass value for each group. Spatially the fleets operated in the same way among all 24 model simulations.
Policy performance metrics. -Masi et al. (2016) found that reef fish catch, Gag biomass, and biodiversity metrics are good indicators for tracking changes in ecosystem dynamics caused by fishing in the Gulf of Mexico. Thus, the solutions from the two-point HCR scenarios were used to evaluate ecosystem-level tradeoffs in both fishery and ecological performance and to single out the ecological performance of Gag biomass under the two-point HCR for all eight F mult scaler scenarios. To compare the low F scaler scenarios to the high F scaler scenarios, the fishery performance is evaluated using the reef fish complex catch (in tons). Here, reef fish complex catch equals the combined catch of the six reef fish groups per year averaged over the last 10 years of the simulation period (at equilibrium). We used the equilibrium yield curves for each reef fish complex group to estimate F at maximum sustainable yield (MSY), F msy ;, under the two-point HCR solutions by fitting a thirdorder polynomial to both the catch and biomass trends for the two-point HCR solutions. We averaged outputs at equilibrium to summarize and compare the results among simulations and to illustrate what the stable ecosystem state is under a given policy. This is important, as fishing can induce alternative trophic pathways that can lead to unique, stable, ecosystem states (Cury et al. 2000).
To illustrate ecosystem-wide shifts in abundance under different levels of fishing mortality, we further aggregated most of our Gulf of Mexico Atlantis model fish and invertebrate functional groups into six distinct species guilds: (1) the assessed reef groups, (2) all reef fish, (3) forage fish, (4) pelagic fish, (5) demersal fish, and (6) a shrimp, crab, and benthic invertebrates guild. Masi et al. (2016) found that the menhaden functional group dominates the forage fish guild because of the proportionately high biomass of Gulf Menhaden Brevoortia patronus in the Gulf of Mexico Atlantis model. Thus, we omitted menhaden from the forage fish guild to illustrate the impact of these harvesting policies on the remaining forage groups within the forage fish guild. Masi et al. (2016) found biodiversity metrics to be robust indicators for tracking the impact of F on the Gulf of Mexico marine ecosystem. Thus, we derived estimates of Kempton's Q (Kempton and Taylor 1976;Ainsworth and Pitcher 2006) for the low and high F mult scaler scenarios. To evaluate the changes in biodiversity in the ecosystem we compared Kempton's Q over time and at equilibrium between the low F and high F scaler scenarios. This biodiversity metric provides a combined measure of species richness and evenness, as evaluated by Masi et al. (2016). We also quantified biomass and catch (in tons) for the reef fish complex and individually for Gag (at equilibrium).
Bracketing uncertainty.-To bracket the uncertainty that stems from interspecific trophic interactions, we randomly sampled from Dirichlet distributions that were fit to observational diet data in Masi et al. (2014). This diet-randomizing methodology was used for all diet observations obtained from stomach samples in Masi et al. (2014) and in a followon study by Tarnecki et al. (2016). Diets for some Atlantis functional groups were taken from other literature sources (these are also described in Masi et al. 2014). We resampled these assuming a normal error distribution with a wide CV (SD/mean) of 0.4. The diets were randomized in 10 independent random draws and applied to 10 new two-point HCR simulations for both the F mult × 0.5 (a low F scaler) and the F mult × 3 (a high F scaler) two-point HCR scenarios. From the 10, randomized, diet runs we derived the mean and associated 95% confidence limits for the biomass and catch (in tons) under both a low and high F scaler scenario for comparison. We chose to compare F mult × 0.5 and F mult × 3 in this analysis because these two solutions were the most discrete with respect to one another.

Fishery Performance: Tradeoffs in Catch
Although we analyzed our two-point HCR scenarios at different upper biomass thresholds, we found that increasing B up from 40% to 60% (or even 80%) of the reference biomass level in 2010 had little effect on fishery performance of the reef fish complex. Thus, we show only analyses that use a B up of 40% of the 2010 biomass estimates. Although there were no strong differences at the reef fish complex level, altering the upper biomass threshold did affect Gag specifically (these results are discussed below).
Although varying the B up had little impact at the reef fish complex level, varying F mult (the F scaler) did affect the fishery performance of the two-point HCR solutions. Here, we show biomass versus catch for the reef fish complex, but excluded the deep Serranidae and Vermilion Snapper groups in this part of the analysis to highlight the response of the other four groups in the reef fish complex (Gag, Red Grouper, Red Snapper, and Lutjanidae). The results suggest that higher levels of fishing (an F mult scaler ≥ 3) led to a stable ecosystem state that is more productive and more favorable for these four reef fish groups than the stable state under lower levels of fishing (i.e., F mult scalers ≤ 1.5). The two different stable states achieved in the two-point HCR scenarios are referred to as ecosystem state 1 (ES1; less productive) and ecosystem state 2 (ES2; more productive). They have qualitative differences in ecosystem structure.
Under the higher F mult scalers, particularly in F mult × 3 and higher, the biomass of these four reef fish groups reached its peak value at which the reef fish complex biomass is approximately 450,000 tons (Figure 4). In ES2, the model predicted that the biomass of Gag, Red Grouper, and Red Snapper is greater under all high F mult scaler scenarios (e.g., the maximum biomass under F mult × 3 is 33,000, 22,000, and 100,000 tons, respectively, compared with a maximum of 20,000, 40,000, and 60,000 tons, respectively, under F mult × 0.5). Notably, in ES1 the Red Grouper catch has been reduced to zero (at equilibrium), suggesting that all low F mult scalers (and the associated realized Fs applied under these lower scaler values; see the x-axis of Figure 5 for the realized F values) are unsustainable for Red Grouper (e.g., Figure A.1), whereas in every high F mult scaler scenario all reef fish complex groups maintained sustainable catch levels at equilibrium. We showed that under higher levels of F the model predicts a more Pareto-efficient tradeoff frontier (Munro 2007), where we achieve higher levels of reef fish complex catch yet still maintain similar biomass levels (or higher, in the case of the F mult × 3 scenario) compared with the low F mult scaler scenarios (Figure 4; e.g., Figure A.1). Here, Pareto efficiency is defined as the circumstance in which high levels of reef fish complex catch cannot typically be obtained without lowering biomass levels (e.g., the F mult × 3 scenario shows higher biomass levels at higher catch removals). The shift in the Pareto efficiency frontier is based on large-scale changes in ecosystem structure that are caused by interspecific interactions and driven by increased productivity in the reef fish complex (as described in the Ecological Performance section). Differences in species biomass are also described in the Ecological Performance section.
Notably, the realized F per species (Figure 5, x-axis) among simulations ranges from approximately 0 to 0.7, except for deep Serranidae and Vermillion Snapper for which the realized F ranges from~0 to 0.15. These realized Fs are comparable with realistic ranges of fishing mortality applied to these species in recent years (SEDAR 2006(SEDAR , 2013a(SEDAR , 2013b.We plotted the reef fish complex by group to show how catch ( Figure 5, closed circles) and biomass (Figure 5, open circles) changes (at equilibrium) among the assessed reef fish functional groups over the eight F mult scaler scenarios tested. Besides the deep Serranidae group, every group in the reef fish complex reached a peak biomass level under the higher, F mult × 3 scenario ( Figure 5, bluehighlighted circles). These results suggest that moderate levels of fishing mortality (note the realized Fs on Figure 5, x-axis) drive declines in the abundant deep Serranidae group, but this in turn could increase productivity and abundance of other species in the reef fish complex.
As an example of model performance, the two-point HCR solutions predict that F msy for deep Serranidae, Red Snapper, and Vermilion Snapper is approximately 0.02, 0.21, and 0.03/year, respectively. The estimated F msy for Snowy Grouper Epinephelus niveatus, an aggregate species in the deep Serranidae group, is 0.05/year (SEDAR 2013b), for Red Snapper is 0.53/year (SEDAR 2013a), and for Vermilion Snapper is 0.33/year (SEDAR 2006). It is typical for multispecies models to predict lower F msy values than single-species models (Walters et al. 2005;Link et al. 2012), as a plethora of ecosystem components are being considered concurrently within a whole-system model, like the Atlantis model we used.
As an example of adjustments to fishing mortality as a function of stock size (following the two-point HCR) we showed the realized F in the F mult × 0.5 scenario as Red Snapper biomass declines below the upper biomass threshold ( Figure 6). The F current remains low (F = 0.06-0.08/year) in consecutive years, as biomass has not recovered. Thus, the realized F from the model continues to decrease as stock size declines (until biomass reaches the B low ).
Inclusion of uncertainty in HCR model predictions.-In the F mult × 0.5 scenario, the lower limit of biomass for the reef fish complex from the 10 random diet draws is 432,863 tons, and the upper limit is 613,976 tons, with a mean biomass estimate of 555,670 tons. In the F mult × 3 scenario, the mean biomass for the reef complex over all 10 random  draws is comparable (476,138 tons) with the low F mult solution, with an upper limit of 521,152 tons and a lower limit of 437,505 tons. The biomass is slightly higher under the low F mult solutions on average, as the reef fish complex is dominated by the deep Serranidae group. Under the low F scenario, the deep Serranidae group has higher biomass when the realized F is lower (discussed in the Ecological Performance section). On average, we found catch in the high F mult scaler scenario to be much higher in all 10 random draws than in the low F mult solution (28,802 tons compared with 15,567 tons). Thus, the solutions from the 10, randomized, diet runs compare well with our initial model predictions (Figure 4).
Evaluating fishery performance for Gag.-We showed that the Atlantis model still predicts the same distinct separation in ecosystem states between the four low F mult and the four high F mult scaler scenarios (compare Figure 4 with Figure 7) when assessing just the Gag fishery performance. The model also predicts we could achieve a higher MSY -780 tons versus 730 tons-for Gag under a more aggressive, two-point HCR (i.e., a two-point HCR with an upper biomass threshold of 40% compared with 60% or even 80%) (Figure 7). A higher MSY is attainable because the productivity of the reef fish complex has increased under relatively higher levels of fishing mortality (as described in Ecosystem Performance section).

Ecological Performance: Tradeoffs in Biomass and Biodiversity
Averaging the biomass of each species guild across the four low F mult scaler scenarios and the four high F mult FIGURE 5. Biomass (open circles) and catch (closed circles) of each assessed functional group shown (to scale for comparison purposes) over the realized fishing mortality rate, for both the low and high F mult scaler scenarios, where each point on the figure represent a different F mult scenario. The blue-highlighted point represents the biomass under the F mult × 3 scenario. Biomass, catch, and realized F are averaged across the last 10 simulation years (at equilibrium). Solutions are used to produce equilibrium catch curves. The high F mult scaler scenarios predict higher catches at higher biomass levels for all reef complex groups, besides deep Serranidae.  Figure 4, where these two groups were omitted), which is why this guild shows higher biomass under the low F mult scaler scenarios. The deep Serranidae dominates the reef fish complex biomass, so under the low F mult scaler scenarios the deep Serranidae group has a much higher biomass at equilibrium; the biomass of deep Serranidae at equilibrium under the F mult × 3 scenario (in ES2) is 399,805 tons, whereas under F mult × 0.5 (in ES1) it is 1,928,312 tons ( Figure A.1). However, under the low F mult scaler scenarios the biomass of the less abundant stocks (Gag, Red Grouper, and Red Snapper) is much lower than in ES2 ( Figure 5). For example, the biomass of Gag, Red Grouper, and Red Snapper at equilibrium under F mult × 3 (in ES2) is 147,548 tons, whereas in F mult × 0.5 (in ES1) it is only 16,066 tons (Figure A.1). This suggests that varying F introduces strong interspecific interactions between the deep Serranidae group and the lower biomass groups (Gag, Red Grouper, and Red Snapper) in the reef fish complex.
In all cases, the increase in the biomass of the reef fish complex at the end of the scenarios with higher fishing can be attributed to a higher available biomass of their prey ( Figure 2). In both high and low F mult scaler scenarios, the reef complex is fished down in the first few years (Figure A.1), particularly stocks with high catch : biomass ratios (i.e., Gag, Red Grouper, and Red Snapper). During this early period, the pattern of fishing determines the amount of stock productivity later in the simulation. In high F mult scaler scenarios, the groups are fished down immediately, and then the realized F gets scaled back in years 2-4, whereas in the lower F mult scaler scenarios the realized F remains low but constant throughout this same period. With less predator FIGURE 7. Derived equilibrium yield curve for Gag under the twopoint HCR scenarios, where catch and biomass are averaged over simulation years 40-50 (at equilibrium). A higher Gag catch and biomass (tons) was achieved in the higher F mult scaler HCR scenarios than in low F mult scaler scenarios, and a greater maximum sustainable yield (MSY) was achieved using a more aggressive two-point HCR (i.e., an upper threshold of 40% is more aggressive than one set at 60% or 80%).
FIGURE 8. Combined biomass of each species guild across the four low F mult scaler scenarios (ES1) and the four high F mult scaler scenarios (ES2) at equilibrium (averaged over years 40-50). Model functional groups have been aggregated into ecological guilds to illustrate the impact of fishing ecosystem-wide. Note, the forage fish guild includes the Pinfish, small pelagic fish, and medium pelagic fish functional groups. The assessed reef groups includes only the six reef fish complex groups, whereas the reef fish guild includes nontargeted reef fish groups as well. 32 biomass under the higher F mult scaler scenarios (i.e., the large, carnivorous reef fishes in the complex) in those first few years, there is a predation release on the shrimp groups (e.g., predation mortality on the other shrimp group drops from 0.39 in the F mult × 0.5 scenario to 0.13 when F mult × 3) and the biomass of the shrimp groups increases by fiveand sevenfold in the other shrimp and white shrimp groups, respectively) (Figure 2; Figure A.2). Within the model, these shrimp groups are the top prey items for the Lutjanidae group, and the other shrimp group is the second highest prey contributor to the Pinfish Lagodon rhomboides group ( Figure 2). Thus, at around year 5 of the simulation both the Lutjanidae and Pinfish groups increase in biomass. Pinfish and Lutjanidae are stop prey items for the deep Serranidae group (Figure 2). In the low F mult scaler scenarios, deep Serranidae have higher biomass (a maximum of 3 × 10 6 tons) than in the high F mult scaler scenarios (a maximum of 1 × 10 6 tons), thus limiting the increase in Pinfish biomass. However, when the realized F is higher (under higher F scaler values) there is more Pinfish biomass (a maximum of 1.4 × 10 6 tons, compared with a maximum of 8 × 10 5 when F is low) available and less deep Serranidae biomass; thus there is less competition for the Pinfish and Lutjanidae prey resource (Figure 2). This leads to an increase (around year 7) in small pelagic biomass ( Figure A.2), as Pinfish is their top prey item within the model (Figure 2). By year 20 in the high F mult simulation, there is an abundant supply of the reef complex's top prey items-Lutjanidae, shrimp, Pinfish and small pelagic fish groups (2.5 × 10 6 , 3.08 × 10 6 , 1.4 × 10 6 , and 40 × 10 6 tons, respectively)-allowing the reef fish complex to recover at a more productive ecosystem state (ES2) than is achievable under any of the low F mult scaler scenarios (see years 8+; Figures A.1, A.2).
We found that both the low and high F mult tow-point HCR solutions predicted declines in biodiversity, as we removed the large, carnivorous reef fish predators under any level of fishing mortality. As an example, the annual average Kempton's Q value under a high F mult (e.g., F mult × 3) scaler scenario is 9.05, whereas under a low F mult (e.g., F mult × 0.5) scaler scenario the annual average is only 8.47. This decline in Kempton's Q indicates a reduction in the number of high biomass functional groups. At equilibrium, we find the Kempton's Q value from the F mult × 3 solution to be 10.27 and only 8.14 under the low F mult (e.g., F mult × 5) scenario. Overall, our results indicate that under the four higher F scenarios we derived a steady state that is more biodiverse (a Kempton's Q of 10.27 versus 8.14), with more groups at high levels of biomass (Figure 8), than what was achievable under any of the four lower F scaler scenarios.

DISCUSSION
Two-point HCRs are attractive management tools because the upper and lower biomass thresholds can be agreed upon ahead of time by both managers and stakeholders. This allows for explicit planning and the ability to address a variety of needs for the fishery, such as maintaining consistency in quotas, minimizing extinction risk, or maximizing revenue (Parma 2002;Deroba and Bence 2008;Punt et al. 2008). However, as we demonstrated, harvest rules applied across a number of competing species may have unforeseen consequences that can alter the expected outcome of a harvest program. The integrated, ecosystem-modeling methodology used here can help quantify those risks and benefits and predict nonlinear ecosystem responses. Results from our modeling-and potential hypotheses to test via future empirical and modeling analyses-suggest that in the case of the Gulf of Mexico these nonlinear ecosystem dynamics can potentially benefit both stock abundance and fisheries catches for trophodynamically linked reef fishes (Masi et al. 2014;Chagaris et al. 2015;Tarnecki et al. 2016).
We showed that removing more of the large, carnivorous reef fish predators early in simulations resulted in a more Pareto-efficient tradeoff frontier, with higher reef fish biomass and higher catch (simultaneously), at equilibrium. These results suggest that adaptively applying different levels of fishing mortality under a blanket, two-point HCR can shift the state of the ecosystem to one that is more productive overall. Investigation into the predation mortality outputs showed that the improved performance of the twopoint HCR simulations using the higher F scaler is driven by a strong "cultivation effect" (Walters and Kitchell 2001), where reducing the top predator abundance in the short term allows for increased productivity in the reef fish complex groups in the long term. This finding is largely driven by reductions in carnivorous, reef fish predators in the first few years (particularly the competition between the deep Serranidae and Lutjanidae groups and the lower biomass groups in the reef fish complex, i.e., Gag, Red Grouper, and Red Snapper, which are more heavily exploited) under the higher F scaler scenarios. In the Gulf of Mexico Atlantis model, Gag, Red Grouper, and Red Snapper are larger bodied, less productive fish than the generic Lutjanidae functional group. We found our results to be robust to uncertainty in the diets that govern these trophic pathways. A similar analysis looking at the impact of trophodynamics on the catch and biomass of reef fishes in the Gulf of Mexico predicted that removing large predators (e.g., Gag) benefited almost every other reef fish in the ecosystem, particularly Black Sea Bass Centropristis striata (Chagaris et al. 2015). Black Sea Bass are included within our deep Serranidae group in Atlantis because we found them to have strong competitive effects on other less-abundant reef fish complex groups.
As noted, the F mult scaler values we applied in each scenario tested are arbitrary values. What is important are the realized F values, which are obtained from Atlantis simulation outputs. The realized Fs differ from the F scalers because the available biomass of an assessed group is spatially and temporally constrained over time. Our realized Fs reflect the ranges of Fs used in single-species stock assessments for some of our reef fish groups (SEDAR 2006(SEDAR , 2013a(SEDAR , 2013b. Interpretation of our results should focus on the realized Fs, and not the arbitrary F scaler values. In this application, biomass was derived annually for each of the six assessed reef fish groups using perfect knowledge. However, it is unlikely that fisheries managers could perform stock assessments each year and adjust F values accordingly for the subsequent year (i.e., accounting for gear selectivity, implementation error, or other factor). The potential benefits, though, in terms of increased yield, biomass, and biodiversity are quantifiable with our simplified simulation approach. For these reasons, the simulations presented here represent the theoretical maximum benefit offered by this mode of management. Future research will consider two-point HCRs that assume imperfect knowledge by managers, gear selection effects, and environmental variation. Simulations testing less frequent (and therefore more realistic) adjustments of fishing rates are also warranted, though Gr€ uss et al. (2016) noted in a similar study that a two-point HCR was largely insensitive to the frequency of the assessment.
Typically, fisheries managers respond to reductions in the available biomass of a stock by implementing a singlespecies rebuilding plan, similar to how the two-point HCR framework is implemented here. However, our simulations provide a quantitative perspective on the role of interspecific interactions in sustaining the biomass of competing species and the potential impact that different harvesting policies may have across the whole marine ecosystem. We showed that competition within the reef fish complex and shifts in system productivity and the forage base may drive the effectiveness of reef fish management plans for this region.