Exploring effects of hypoxia on fish and fisheries in the northern Gulf of Mexico using a dynamic spatially explicit ecosystem model

The formation of an extensive hypoxic area off the Louisiana coast has been well publicized. However, determining the effects of this hypoxic zone on fish and fisheries has proven to be more difficult. The dual effect of nutrient loading on secondary production (positive effects of bottom-up fueling, and negative effects of reduced oxygen levels) impedes the quantification of hypoxia effects on fish and fisheries. The objective of this study was to develop an ecosystem model that is able to separate the two effects, and to evaluate net effects of hypoxia on fish biomass and fisheries landings. An Ecospace model was developed using Ecopath with Ecosim software with an added plug-in to include spatially and temporally dynamic Chlorophyll a (Chl a) and dissolved oxygen (DO) values derived from a coupled physical–biological hypoxia model. Effects of hypoxia were determined by simulating scenarios with DO and Chl a included separately and combined, and a scenario without fish response to Chl a or DO. Fishing fleets were included in the model as well; fleets move to cells with highest revenue following a gravitational model. Results of this model suggest that the increases in total fish biomass and fisheries landings as a result of an increase in primary production outweigh the decreases as a result of hypoxic conditions. However, the results also demonstrated that responses were species-specific, and some species such as red snapper (Lutjanus campechanus) did suffer a net loss in biomass. Scenario-analyses with this model could be used to determine the optimal nutrient load reduction from a fisheries perspective. ublis © 2015 The Authors. P


Introduction
Nutrient rich waters flowing from the Mississippi River into the Gulf of Mexico result in high primary productivity in this coastal area . Bacterial decomposition of this organic matter in combination with summer stratification has led to the occurrence of an extensive area of low bottom oxygen since at least the early 1970s . While often referred to as the 'dead zone', the effect on living marine resources of this annually reoccurring area of hypoxic bottom waters off the coast of Louisiana is not necessarily lethal.
Hypoxia refers to oxygen levels of 2 mg/l or lower, which can lead to decreased feeding and growth rates, changes in activity level, avoidance behavior, and death in fish and shellfish (Bell and Eggleston, 2005;Robert et al., 2011;Goodman and Campbell, 2007). The exact level of dissolved oxygen that results in effects on physiology or behavior is species-specific, which can results in community structure shifts and changes in species interactions (Essington and Paulsen, 2010). Indirect effects occur through predator-prey relationships; fish could be affected not by hypoxia, but by the response of their prey or predators to hypoxia, and the effects could be either positive or negative (Altieri, 2008;Pierson et al., 2009;Eby et al., 2005). Effects on fisheries may be even more complicated, as catch per unit effort (CPUE) could decrease when the abundance of target species is reduced by hypoxia, or could increase due to aggregation of target species at the edge of the hypoxic zone, which may enhance their susceptibility to be caught (Craig, 2012).
A significantly obscuring mechanism is the fact that the same nutrient enriched waters that are the main cause of bottom hypoxia (Rabalais and Turner, 2001), are responsible for the high primary and secondary production in this region (Gunter, 1963;Nixon and Buckley, 2002;Chesney et al., 2000). It is likely due to these complications, that holistic effects of hypoxia on the fisheries ecosystem of the northern Gulf of Mexico have remained elusive (Rose, 2000;Rose et al., 2009).
The purpose of this study is to analyze effects of hypoxia on fish and fisheries through ecosystem model simulations, and to provide a tool that can be used in management scenario analyses pertaining to Mississippi River nutrient load reductions and coastal fisheries management. To this purpose an Ecospace model was developed using Ecopath with Ecosim (EwE) software that was enabled to receive spatio-temporal primary productivity and dissolved oxygen output from a coupled physical-biological hypoxia model developed by Fennel et al. (2011). Since a reduction in hypoxia would entail a reduction in nutrients that enter the Gulf of Mexico, it is important to incorporate the effects of nutrient enrichment on phytoplankton (and changes therein) in an ecosystem model that studies effects of hypoxia and scenarios that may reduce this hypoxia. Output of the Fennel et al. (2011) model of dissolved oxygen (DO) as well as Chl a was used as forcing functions in the Ecospace model to account for both effects. Similar approaches to incorporate effects of biogeochemistry on foodweb models, often referred to as End-to-End modeling, have been used in other studies (see e.g. Libralato and Solidoro, 2009).
The ecosystem model developed for this study takes a holistic approach by simulating species interactions, while accounting for changes in biomass as well as spatial distribution changes, and by explicitly simulating fisheries with dynamic fleets. The model allows for simulations of all direct and indirect effects on fish and fisheries, in an environment where hypoxia and primary productivity fueling can be evaluated together and separately. While this ecosystem model contains sixty groups to provide a representative simulation of the ecosystem, the main focus of this paper is on a select group of species that are of economic or ecological significance. These species are Gulf menhaden (Brevoortia patronus), which is largest fishery in Louisiana by weight; brown, white and pink shrimp (Farfantepenaeus aztecus, Litopenaeus setiferus, and Farfantepenaeus duorarum), together comprising the largest fishery by value; red snapper (Lutjanus campechanus), a popular sportfish; Atlantic Croaker (Micropogonias undulatus), the most dominant forage fish in the model area; and jellyfish, a group of organisms of interest because of previous documented responses to hypoxia in other areas.

Data preparation
Fisheries independent survey data from the SEAMAP program of the Gulf States Marine Fisheries Commission (seamap.gsmfc.org) was used to determine which species were representative of the area, and to determine the biomass of each species present in the model area. Initial biomass in the base model was based on the average biomass of each group (species or functional group) from 2005 to 2008. Fishing was represented by including shrimp trawls, recreational fishing, snapper/grouper fishery, crab pots, menhaden fishery, squid fishery, and longlines as 'fleets' in the model. Annual landings of model groups by these fleets were based on NOAA Fisheries Annual Commercial Landings Statistics (st.nmfs.noaa.gov), and trip ticket data from the Louisiana Department of Wildlife and Fisheries. These data were used to develop the Ecopath model.
Landings data from 1950 to 2010, and SEAMAP data collected in the model area from 1982 to 2010 were used to calculate annual landings and biomass (t/km 2 ) respectively for each group in the model for which these data were available. In addition, an oxygen forcing function was developed from data collected during Lumcon cruises from 1998 to 2007 (D. Obenauer, personal communication), and a nutrient forcing function from NO x data collected in the Mississippi River by USGS from 1950USGS from to 2010 to simulate nitrogen load into the coastal area from the Mississippi River. These time series and forcing functions were used for model calibration in Ecosim.
In EwE, a nutrient forcing function serves as a multiplier on primary production. In order for groups to respond to the level of dissolved oxygen, empirically derived sigmoidal oxygen response curves were developed. These curves were developed by determining catch rates at each level of dissolved oxygen, using all SEAMAP data where dissolved oxygen was measured during collections. The tolerance curves were then used as a multiplier on effective search rate in Ecosim (and Ecospace, using a plug-in described in Section 2.5) as described in Christensen et al. (2008) andde Mutsert et al. (2012), to affect biomass of each specific group (Fig. 1).

Model preparation
The EwE modeling suite was used to build the model (www. ecopath.org). The virtual representation of the ecosystem was developed in Ecopath, the static model of the EwE modeling suite. Groups in the model represent single species as well as species aggregated in functional groups. Where deemed necessary to represent ontogenetic diet changes or size-selective fisheries, species were split into multiple life stages. For those species, the initial biomass of only one life stage was derived from empirical data, and the biomass of other stages were determined using a von Bertalanffy growth model. Some functional groups were represented with multiple life stages as well. This resulted in 60 groups (Table 1). Parameters included for each group to develop a mass-balanced Ecopath model in addition to biomass (B), were the P/B (production/biomass) ratio, Q/B (consumption/biomass) ratio, and the total fisheries catch rate (Y) for the groups that are fished. Parameters were derived from other Gulf of Mexico food web models de Mutsert et al., 2012) or fishbase (fishbase.org).
Two master equations must be satisfied to correctly parameterize the Ecopath model. The first equation describes the production of each functional group as a set of n linear equations for n groups: Initial conditions of mass-balanced Ecopath model. B = biomass, Z = total mortality, P/B = production to biomass ratio, Q/B = consumption to biomass ratio, EE = ecotrophic efficiency.
Nr. Group name where (P i /B i ) is the production to biomass ratio for group i, EE i is the ecotrophic efficiency (the proportion of production used in the system), B i and B j are the biomasses of the prey and predators respectively, (Q j /B j ) is the consumption to biomass ratio, DC ji is the fraction of prey i in predator j's diet, Y i is catch rate for the fishery for group i, E i is the net migration rate, and BA i is the biomass accumulation for group i. The Ecopath model assumes conservation of mass over a year. Energy balance within each group is ensured with the second master equation: where production can be described as: Production = predation mortality + catches + net migration + biomass accumulation + other mortality More succinctly, production can be described by the following equation: where P i is the production of prey group i, Q j is the consumption of predator j, DC ji is the diet composition contribution of i to j's diet, F i is the instantaneous rate of fishing mortality, NM i is the net migration rate of prey group i, BA i is the biomass accumulation rate for i, M0 i is the other mortality rate for i (non-predation, nonfishery), and B i is the biomass of i. In addition, a diet matrix was constructed based on diet information from stomach content analysis from nekton collected in the model area (A. Adamack, unpublished data), supplemented by information available in the literature. To achieve mass balance, the diet matrix was adjusted to attain a plausible solution of the flow of biomass through the foodweb. The available diet information usually did not provide exact proportions of each diet item, which made the diet matrix the most suitable component to adjust in order to achieve mass-balance. For example, when previous studies indicated that a specific species was the dominant prey species for a predator, the exact proportion of this prey item was adjusted during the mass-balancing procedure while still maintaining its status as dominant prey item. The diet matrix is provides as supplemental material 1.
During the mass balancing procedure in Ecopath, the model calculates Ecoptrophic Efficiency (EE) of each group, which represents the amount of biomass of that group used in the system . A mass-balanced solution of the model is presented in Table 1.

Spatial components
A model area of 44,890 km 2 was chosen, which encompasses the Louisiana coastal zone and the annually recurring hypoxic zone. This area was represented in Ecospace with 5 km 2 grid cells, and is the model area of our Ecospace model, which we have called the NGOMEX (Northern Gulf of Mexico) ecosystem model (Fig. 2).
For the spatial and temporal model simulations, dissolved oxygen and Chl a output from 1990 to 2004 of a coupled physical-biological hypoxia model (Fennel et al., 2011) was used as forcing functions. Chl a levels in Fennel et al. (2011) are affected by the nutrients entering the coastal zone from the Missisippi River and other freshwater sources. This output was averaged by month and matched to the Ecospace grid map so that one value of bottom dissolved oxygen and one of Chl a could be read into Ecospace per month per grid cell during a model simulation. In the few occasions where the model area of Fennel et al. (2011) did not overlap with our model, DO and Chl a output was extrapolated from nearby cells. This was done for the estuaries, while the focus area for our modeling effort had 100% overlap. Example DO output from Fennel et al. (2011) that is used as a spatial-temporal forcing function is shown in Fig. 3. A plug-in to Ecospace was used to read in this spatial-temporal forcing function (see Section 2.5 for more details). Dissolved oxygen affected the groups in the model as stipulated by the response curves, while Chl a was used as a driver of phytoplankton biomass, assuming a linear relationship.
Two non-dynamic habitat features were included in the spatial model, depth and salinity area. Depth was based on the bathymetry of the model area; depth ranges were included to ensure (adult) offshore species would not enter shallow estuarine areas if they are not known to do that. While salinity is not modeled dynamically in this model, a 'marine', 'estuarine' and 'freshwater' zone is described loosely based on existing salinity gradients in the model area. While the focus of this model is on the marine coastal zone, these habitat features prevented species to escape coastal hypoxia by fleeing to areas that are too shallow or too fresh for them to enter in real life. A conceptual model of the NGOMEX ecosystem model is presented in Fig. 4.

Model calibration
Temporal dynamic simulations were performed in Ecosim, the time-dynamic module of EwE, to calibrate the model. DO and NO x were included in the calibration runs as environmental forcing functions based on data described in Section 2.1. The level of dissolved oxygen affects the effective search rate of species in the model as described by the response curves in the same manner as salinity affected species in de Mutsert et al. (2012). The model was calibrated against biomass time series and landings data as described in Section 2.1. During calibration, the model was iteratively fitted to landings and biomass time series data by making vulnerability exchange rate adjustments until the smallest sum of squares (SS) was found using the fit-to-time-series feature in  EwE . Following the Foraging Arena Theory described in Walters and Martell (2004), each group is present in the model in a vulnerable (to predation) and invulnerable state. The vulnerability exchange rate determines how quickly the mass of a group can switch between those states, where high numbers (around 100) indicate Lotka-Volterra predator-prey interactions (all prey is vulnerable to predation because of the high exchange rate between the vulnerable and invulnerable portion), and low numbers (around 2) indicate a significant portion of the group is unavailable to predation. We used the fit-to-time series procedure to determine the vulnerability exchange rates that resulted in the best fit of model predictions to biomass and landings data. The metric used to determine model fit was the following: where SS is sum of squares, nts is the number of time series loaded, nobs i the number of observations in time series i, w i is the weight of the time series i (all time series weighted equal in our model), o it is the observed value in time series i at time step t and p it is the Ecosim predicted value for variable i at time step t. Including DO and nutrient loading (in the form of NO x ) as environmental forcing functions in Ecosim improved the fit of the model to time series, and decreased the total SS for all fits. Fig. 5 shows fits to time series (with SS) of a selection of species that are highly abundant in the area and/or have economic or ecological significance. The vulnerability exchange rates that were altered during this calibration procedure were carried over to Ecospace.

Model simulations
After calibration, spatial simulations were performed in Ecospace, the spatial-temporal module of EwE. In the new habitat foraging capacity model of Ecospace, dispersal rates of groups into a cell are affected by the cell suitability/capacity (Christensen et al., 2014). If the neighboring cell has a lower capacity then the dispersal rate to the cell will be proportional to the capacity difference. For example, if the capacity of a cell is 0.5 for a specific group, the maximum movement rate into this (in this model set to 300 km/yr for all groups) was adjusted by this proportion. The capacity of a cell was based on DO and habitat (depth and salinity area as described in Section 2.3). Fleets are dispersed by a gravitational model based on profitability per cell. Profitability per cell is based on the biomass of the target group(s) of a fleet, the price per pound of each target group in 2010, and the distance from port (fuel cost). Two ports with the highest landings in Louisiana were included in the model, Empire-Venice and Intracoastal City (www. oceanomics.org; Fig. 6).
To loosely link the physical-biological hypoxia model from Fennel et al. (2011) to Ecospace, a plug-in was added to the EwE source code. The plug-in reads in a DO and Chl a value per grid cell per time step (5 km −2 month −1 ). This provides for spatial and temporal variation in the effective search rate and primary production. The DO values are fed into the environmental response functions defined in Ecosim. The values returned by the environmental response functions act as a forcing multiplier on the rate of effective search. This facility, provided by the plugin, works in the same manner as an Ecosim forcing function that has been applied to search rate . The Chl a data is used to update the Ecospace Relative PP spatial layer, which allows for spatial shifts in primary production over time. The Ecospace Relative PP layer is a multiplier that is used to scale the primary production relative to the base productivity of the Ecopath model. During initialization the values in the Relative PP layer are normalized to scale the spatially averaged Ecospace productivity to the Ecopath base productivity rate . The values read by the plug-in can shift from this baseline value to increase or decrease the spatially averaged productivity over time.
Scenarios simulated were 'no forcing', which simulated a coastal environment without nutrient fueling from the Mississippi River (or any other source of added nutrients) but also no formation of a hypoxic zone; 'enrichment only', which simulated nutrient loading effects on primary productivity, but where hypoxia had no effect on any organism; and 'enrichment + hypoxia', which included primary productivity forcing, and effects of DO (and thus hypoxia for part of the year) on fish biomass. Each scenario was run from 1950 to 2010; results presented reflect the output from simulation year 2010. While sixty groups were simulated, results are presented of a select group of species that are of economic or ecological interest. Fig. 6. Location of ports in the NGOMEX ecosystem model, representing Intracoastal City on the left, and Empire-Venice on the right (black dots). The coloration indicates distance from port, which is included in the calculation of fisheries revenue.

Results
Biomass and landings output of the scenarios 'no forcing', 'enrichment only', and 'enrichment + hypoxia' was compared. The scenario 'enrichment + hypoxia' simulates the real world scenario of Chl a concentration fueled by nutrient loading, and seasonal hypoxia in the coastal zone. The scenarios were run from 1950 to 2010, and output is presented as relative change, which is the change in biomass or landings of each group from the initial biomass or landings. The initial biomass and landings were the same for each scenario, so the scenario outcomes can be compared to each other. When looking at total landings and biomass, results indicate that the seasonal presence of hypoxia reduces both landings and biomass as compared to the 'enrichment only' scenario ( Fig. 7). However, both 'enrichment only' and 'enrichment + hypoxia' had much higher increases from initial biomass and landings than the 'no forcing' scenario; the latter even showed a small decrease. The difference between 'enrichment only' and 'enrichment + hypoxia' is comparatively so small that these simulations suggest that the decrease in secondary production due to Fig. 7. Total landings and total biomass results of three scenarios (no forcing, enrichment only, and enrichment + hypoxia) that ran from 1950 to 2010. The relative change from the same initial conditions is presented of total biomass and total landings, species-specific biomass of selected species (B), and catch from all fleets (C). hypoxia (in indirect effect of nutrient loading) is trivial in comparison to the increase in secondary production due to the bottom up effect of nutrient loading. Overall, there was a 33% increase in total landings in the 'enrichment + hypoxia' scenario as compared to the 'no forcing' scenario, and a 13% increase in total biomass. Removing hypoxia only increased that amount by an extra 5% and 0.6% respectively.
While total landings and biomass show the concurring trend of a small decrease in the 'no forcing' scenario, and large increases in 'enrichment only' and 'enrichment + hypoxia', individual groups vary in their response (Fig. 7). The biomass of common species in Louisiana; Gulf menhaden (Brevoortia patronus), Atlantic croaker (Micropogonias undulatus), and shrimp (brown shrimp -Farfantepenaeus aztecus, white shrimp -Litopenaeus setiferus, and pink shrimp -Farfantepenaeus duorarum) showed a response similar to what was seen in total biomass. Red snapper (Lutjanus campechanus) biomass however, decreased in all three scenarios, and decreased most in the 'enrichment + hypoxia' scenario (17.6%), followed by 'enrichment only' (10.4%) and 'no forcing ' (8.3%). An opposite effect was seen in jellyfish, which displayed increases in all three scenarios, and the highest increase in the scenario with hypoxia (7.8%). Changes in landings do follow this pattern for almost all fleets, except for crab and squid fisheries, which see a small increase in landings when hypoxia is added as compared to enrichment alone.

Discussion and conclusion
Our simulations suggest that reductions in landings and biomass due to hypoxia are an order of magnitude lower than increases seen due to the nutrient enrichment (which is the main cause of hypoxia). Some fisheries in the model even experience an increase in landings in the scenario that includes hypoxia, namely blue crab and squid landings. The crab pots are not set in areas affected by hypoxia, which could explain this pattern, while the increase in squid landings is likely an indirect effect, since squid had a slightly higher tolerance for low oxygen as most of its predators, and slightly higher biomass in the scenario with hypoxia as a result. In general, current simulations do not suggest that natural resource managers should take the hypoxic zone into account in fisheries management plans (e.g. by restricting effort during hypoxic events), as the occurrence of seasonal hypoxia in combination with fishing does not lead to unsustainable biomass reductions.
This study emphasizes the importance of the positive bottomup effect of nutrient enrichment on secondary productivity (Nixon and Buckley, 2002). Some notable species that follow the pattern of large increases in biomass as a result of nutrient enrichment, and only a slight reduction in biomass as a result of hypoxia, include Atlantic croaker, which is the most abundant species in this area and knows to have a high tolerance for hypoxia (Bell and Eggleston, 2005), Gulf menhaden, which is the largest fishery in Louisiana, and gulf shrimp (brown, white, and pink shrimp), which is the fishery with the highest revenue in Louisiana.
Still, from these results cannot be inferred that nutrient load reduction is not an important restoration measure, and that it would necessarily reduce secondary productivity. Our scenario of 'no forcing' is not a real-world scenario, and no nutrient reduction plan would conceivably remove all nutrients from the freshwater sources flowing into the Gulf of Mexico. Therefore, the corresponding low secondary production seen in the 'no forcing' scenario would never be attained. In addition, the relationship between nutrient loading, primary productivity, and hypoxia is non-linear and complex (Fennel et al., 2011); a reduction in nutrient load would not necessarily reduce bottom up fueling of the foodweb and hypoxia to the same extent. Momentarily disregarding the 'no forcing' scenario, a consistent small decrease in biomass from the nutrient enrichment scenario to the nutrient enrichment scenario with summer hypoxia can be seen. This small reduction could be ecologically significant for some species.
One species that seems affected by nutrient loading as well as hypoxia in our simulations is red snapper. An increase in mortality due to higher shrimp landings -and thereby higher bycatch of juvenile red snapper -in the scenarios that include nutrient enrichment is a likely cause of a decrease in red snapper biomass in those scenarios. The model reflects the impact shrimp trawling has on red snapper, which has been reported in studies related to red snapper stock status (Cowan, 2011). The additional decrease in biomass when hypoxia is present does indicate a negative effect of hypoxia on red snapper. Weaker recruitment of red snapper in years of severe hypoxia has been observed in a previous study (Switzer et al., 2015).
Another interesting result is the increase in jellyfish biomass in the scenario with hypoxia. Jellyfish, often regarded as nuisance species, likely find refuge from predation in hypoxic areas due to their high tolerance of low oxygen conditions. Increases in jellyfish in response to hypoxia in coastal ecosystems has been predicted or observed in other studies (Breitburg et al., 2003;D'Elia et al., 2003;Miller and Graham, 2012), and could exacerbate hypoxia effects on zooplankton by adding increased predation pressure.
This study concurs with some previous publications that hypoxia typically does not reduce overall fisheries landings or biomass, but that hypoxia should still be addressed in restoration plans . The use of novel spatial-temporal forcing functions in Ecospace allows for more realistic simulations of effects on fish and fisheries of environmental drivers that vary in space and time. Ecosystem models with this capability have only recently been described (Steenbeek et al., 2013;Christensen et al., 2014), but are expected to increase in numbers rapidly. Their usefulness in developing restoration and/or natural resource management strategies, especially when linked to physical/chemical models seems evident, and has already been recognized (de Mutsert et al., 2014a(de Mutsert et al., , 2014b. The model presented in this paper would be useful in restoration planning, and development of management strategies to reduce hypoxia without unacceptable losses to fisheries productivity. Models such as the physical-biological model of Fennel et al. (2011) could be used to simulate effects of nutrient load reductions on hypoxia and primary productivity in the coastal zone. The NGOMEX Ecospace model could then use those results to simulate effects of nutrient reductions on fish and fisheries in a scenario analysis. These loosely coupled models could thereby be used as a tool in nutrient reduction analyses to inform management decisions.

Acknowledgement
KdM would like to thank Katja Fennel for providing model output, and Arnaud Laurent for reformatting the output to fit the NGOMEX ecosystem model grid, and extrapolating output into cells of the model area that were not included in the Fennel et al. (2011) model. KdM would also like to thank Carl Walters for providing input during model development. This research was funded by NOAA's Center for Sponsored Coastal Ocean Research (CSCOR) under grant nr. NA09NOS4780233.

Appendix A. Supplementary data
Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.ecolmodel.2015. 10.013.