Fishing down nutrients on coral reefs

Fishing is widely considered a leading cause of biodiversity loss in marine environments, but the potential effect on ecosystem processes, such as nutrient fluxes, is less explored. Here, we test how fishing on Caribbean coral reefs influences biodiversity and ecosystem functions provided by the fish community, that is, fish-mediated nutrient capacity. Specifically, we modelled five processes of nutrient storage (in biomass) and supply (via excretion) of nutrients, as well as a measure of their multifunctionality, onto 143 species of coral reef fishes across 110 coral reef fish communities. These communities span a gradient from extreme fishing pressure to protected areas with little to no fishing. We find that in fished sites fish-mediated nutrient capacity is reduced almost 50%, despite no substantial changes in the number of species. Instead, changes in community size and trophic structure were the primary cause of shifts in ecosystem function. These findings suggest that a broader perspective that incorporates predictable impacts of fishing pressure on ecosystem function is imperative for effective coral reef conservation and management.


Surveys
All sites were on fore-reef environments 10-15m deep, usually dominated by the coral Orbicella spp. Each sub-region included sites inside and outside marine reserves (i.e., no-take areas were fishing was prohibited), except at Dry Tortugas where only a reserve site was surveyed (Fig. 1, table 1). Site selection was generally haphazard: 35 of the 39 sites were selected either haphazardly or effectively randomly. The exceptions were the four sites in 'Jardines de la Reina' marine reserve, Cuba (Fig. 1, table 1) that were chosen a priori because there were reputed to have relatively high predator biomass 1 . In this way, the sites sampled were representative of fairly intact fish communities as well as heavily fished communities. To minimize seasonal variability, we conducted all surveys May-July in 2010-2012.
Underwater visual censuses, modified from Lang et al. 2010 2 , were used to characterize the fish assemblages. At each site we randomly placed six to eight belt transect sets parallel to the spur-and-groove habitat or along the reef slope formation following constant isobaths. In each transect we recorded fish species, abundance and estimated body size. Fish total length (TL) was binned by 10 cm class size except for individuals <10 cm TL where 5 cm interval was used; it is acknowledged that using size bins may overestimate the biomass of large fish, but this method was determined most accurate after repeatability tests conducted by the surveyers. Two divers swam along a 50 x 10 m transect as it was positioned. The first diver counted fish of medium size (5-40 cm TL) in a 30 x 2 m belt area, followed by a 15 x 1 m belt to estimate small fish <5 cm TL. The second diver counted fish > 40 cm TL within a 50 x 10 m belt to minimize bias on overestimating more mobile and large-bodied fish. The two smaller transects were contained within the largest transect to create a belt transect set. Each transect set was surveyed in ~15 minutes, covered the entire visible water column and were at least 10 m apart. All surveys were conducted from 0800 to 1600hr. Experienced divers performed all fish surveys and fish size estimations were corroborated underwater against known size lengths. No attempt was made to survey fish behavior, a factor that could influence our surveys across the gradient of fishing pressure.
The biomass for each fish species was calculated through the allometric length-weight conversion formula, W = aTL b , where W is body mass in grams, TL is the total length of each fish in cm estimated from visual surveys, and the parameters a and b are species-specific selected from geographical areas close to our study region 3 or determined from personally collected data.
To calculate biomass we used the mid-point of the 5 or 10 cm interval estimates. For all the analysis we used fish biomass because it combines abundance and size resulting in a comprehensive indicator of fish assemblages status across disturbance gradients 1 .
Fish species were assigned to one of seven functional groups: piscivores, piscivoreinvertivores, macroinvertivores, microinvertivores, planktivores, omnivores and herbivores following reported dietary information 3 . For the purpose of this study we considered "fish predators", piscivores and/or piscivore-invertivores because both feed on fish. Invertivores (macro-and micro-) only feed on invertebrates; and omnivores mostly consume marine plants and invertebrates.

Models for nutrient excretion
The excretion and storage models allow processes of C, N and P storage and N and P supply to be estimated as a function of wet weight for 160 fish species. The analysis consisted of two steps: 1) Model all processes onto every fish in the dataset.
2) Quantify aggregate processes for each community and apply these data to hierarchical models.
Bayesian statistics allow parameters to be estimated based on observed distributions (the observed data), and prior distributions that allow knowledge from previous studies to be applied explicitly and quantitatively 4 . In this study, we used Bayesian statistics to develop models that predict excretion rate as a function of wet mass by informing empirical data (the observed data) with bioenergetics models (used to generate the priors), thus incorporating the two most widely applied methods to estimate fish excretion 5,6 into singular models of nutrient excretion by fishes.
The modeling approach was developed such that if the empirical data were robust then the final model would primarily reflect these data (i.e., the priors developed from the bioenergetics model would only minimally inform the output). When the empirical data were not robust, due to lack of individual empirical measurements on rare species or high variability in the data, the final model would then be more of a reflection of the bioenergetics models (i.e., the priors would have more influence on the output). In doing so, this approach allowed us to: (1) underpin extensive empirical data to produce robust models with realistic error, and (2) fill gaps in the empirical dataset for which empirical data was incomplete.
This modeling process consisted of four steps: (1) Bioenergetics models were developed for each family (and in some cases at the genus-level) in our dataset to estimate excretion rates of N and P for a given mass of an individual fish.
(2) These data were run in an initial Bayesian simple linear regression analysis (using uninformative priors), to generate parameter estimates for the slope and intercept of each model (y = mx + b, where y is excretion rate, x is the wet mass of an individual, m is the slope and b is the intercept) (see detailed methods for bioenergetics models below).

(3)
A second Bayesian simple linear regression analysis was conducted using the empirical data. In this case we used the normal posterior distributions (i.e., the mean and standard deviation) for the slope and intercept generated in Step 2 as the priors for the model 4 . In this way we were able to take advantage of all available data and multiple approaches to generate robust estimates of nutrient supply by fishes.

(4)
The posterior distributions of these final estimates for the slope and intercept were then used to calculate the excretion rate for every fish within our survey datasetsee Ecosystem Modeling section for further explanation.

Bioenergetics models
Bioenergetics models use a mass balance approach given a priori knowledge of the natural history (e.g., diet, feeding activity), physiology (e.g., stoichiometry of predator and prey, assimilation efficiency of nutrients, consumption rates, energy density of prey) and environmental conditions (temperature) to estimate nutrient excretion 5,7 . We followed the approach to construct bioenergetics models, and used diet stoichiometry data 8,9 . In total, bioenergetic models for 31 genus and 25 families within our surveys were developed.
Stoichiometry data for each family/genera was determined by averaging the percent nutrient content for a suite of species within the given level of classification. Use of parameters for closely related species may increase error in model estimates 10,11 ; however, empirical work suggests that variation in excretion rates vary little within families but widely among families 12,13 . Energy densities of prey items were obtained from 14 . Assimilation efficiencies, which have been shown to have only marginal influence on model estimates 15 were assumed to be 80% for P and 70% for N based on literature recommendations 10,11 . The growth rate of an animal has been shown to be a particularly influential parameter in bioenergetics 14 , as such published growth rate values were found for each taxon of interest. Other parameter estimates were obtained from literature values specific to the given taxonomic level. Dietary parameters were determined from diet data collected by the authors [16][17][18][19] and from other published data 20 .
We used Fish Bioenergetics 3.0 software 7 to determine consumption rates for the dominant feeding guilds: piscivores, piscivore-invertivores, macroinvertivores, microinvertivores, planktivores, omnivores and herbivores. To do this we chose the taxon per feeding guild for which we had the best parameter estimates (e.g., Lutjanidae >100 individuals, thousands of diet data, etc.) and used the software to calculate consumption rates based on energetic demands of the taxon. These consumption rates were then used for all families within that particular guild, holding this parameter constant and allowing other important estimates (e.g., body stoichiometry, prey stoichiometry and growth rate) to have influence over the model.
Bioenergetics models were run using R software 21 .
To account for inherent error that occurs when parameterizing such models, we propagated uncertainty associated with diet content and consumption rates, two parameters that are highly influential in bioenergetics models 5,7,22 , through the models using Monte Carlo simulations. Specifically, a normal distribution of values was created for each parameter with a standard deviation of 5% of the maximum potential value of that parameter (in both cases the parameters represent a proportion, so the standard deviation was 0.05). For each model run, random draws were taken from within these distributions 500-10,000 times, depending on the size range of the fish within that family. Note the number of draws within this range did not change the outcome of the model.

Empirical excretion estimates
All fish were captured using hook and line or traps on Abaco Island, The Bahamas between 2008-2011. Fish were captured in coral reef, mangrove and seagrass ecosystem types.
Fish were pooled across ecosystem type, such that individuals from a given species could have been caught in any one or all ecosystems. We accounted for potential differences in resource availability across ecosystem type, which would be predicted to affect recycling rates, in two ways: (1) individuals within a given species were often collected from different ecosystem types and potential variation across ecosystem type was pooled, and thus accounted for, in our empirical models, (2) we modeled error for diet nutrient content in our bioenergetics models.
Excretion rates, for nitrogen -NH4 + and phosphorussoluble reactive phosphorus (SRP), were measured in situ following the methodologies of 23 , as modified by 6  hours, analyzed for NH4 using the methodologies of 24 , or frozen for transport to UGA for SRP analyses using the ascorbic acid method and colorimetric analyses 25 .

Bayesian excretion models
Previous research on fish nutrient stoichiometry has shown that variation within families is relatively constrained 12 . As such, we used genus-or family-level bioenergetics models to inform empirical data in a Bayesian framework (i.e., bioenergetics models were employed to constrain excess variance in empirical excretion models when present). To further illustrate this approach, we follow each step taken to generate the final equation (excretion rate = wet mass*slope + intercept) (provided above) with an example species: gray snapper (Lutjanus griseus).
Step 1: A genus-level bioenergetics model for Lutjanus was developed.
Step 2: A Bayesian simple linear analyses was run using the size-specific data generated from the bioenergetics model. All models were run with three chains for 50,000 iterations with a burn-in period of 1000.
Data for excretion models were not transformed and assumptions of normality were met.
Bayesian analysis was run using the "rjags" package in R 21 .
It is important to note that, as is the case with all models, there are potential biases that could influence our estimated excretion rates of coral reef fishes. First, empirical estimates have three primary potential issues: the stress on fish due to capture and handling, and oxygen deprivation that can increase nutrient excretion rates, and starvation, that can reduce excretion rates 6 . We attempted to account for these issues following the recommendations of 6 , as well as conducting our own time trial experiments to determine the most appropriate time to take measurements such that the influence of handling stress, oxygen stress, and starvation is minimized (see 13 ). The primary bias associated with our bioenergetics models is the lack of parameter estimates for each of the families within our dataset. It is for these reasons that we used a Bayesian framework to allow the bioenergetics models to inform the empirical data.
While we believe this approach to be a robust alternative to using solely empirical or bioenergetics models, we nonetheless recognize its limitations. To further account for error in model estimates, we use error propagation when scaling individual rates of excretion to entire fish communities. However, it is recognized that utility of error propagation is nonetheless still limited by any potential systematic bias in the initial model.

Ecosystem modeling
Using the equations generated from the Bayesian models, excretion estimates were Fish nutrient supply is a function of body size, organism identity and diet 5,12 . As such, we used Monte Carlo simulations to model uncertainty into our estimates of fish nutrient supply for individual fish within the dataset. For each fish, we sampled from the posterior distribution of both the slope and intercept from our Bayesian excretion models to calculate 1000 mass-based species-specific excretion estimates. These values were summed to provide a distribution of community-level aggregate estimates (n=1000) of N and P supply. We applied the same methodology to calculate nutrient storage, whereas in this case we sampled 1000 times from the normal distribution (mean ± standard deviation) associated with our stoichiometric estimates for body nutrient content at each taxonomic level (typically genus or species). In doing so, we modeled realistic estimates of error into each step of our analysis to create a range of values that represent a realistic distribution of nutrient supply and storage for every fish and the entire community. Because we sampled from a normal distributions for each estimate, this error propagation approach should not inherently alter the mean value of the aggregate communitylevel excretion, but instead provide information as to how much variability there may be in this mean as a result of potential error and natural variability.

Hierarchical Mixed Effects Models
We used hierarchical mixed effects models and information theory (Akaike information criterion, AIC c ) 27  showed that all were less than r = 0.3 with the exception of Richness and Species evenness in which the correlation was r = 0.53. Biomass has long been recognized to be a strong predictor of stoichiometric properties 31 , and in the case of this study was directly used to calculate our response values. For this reason we did not include biomass as a predictor in our models.
However to test for the importance of biomass we divided all responses by community level biomass to get a community-level mass-specific process, e.g., grams of N storage per grams of biomass.
Multifunctionality was calculated following 32 , whereby we calculated the average of the Z scores for each ecosystem process of interest. Z scores were calculated from log transformed normalized data: Z score = [ x -μ ] / σ ; where x is the site-level ecosystem process, μ is the mean value for all sites, and σ is the standard deviation of all sites. This index was chosen for three primary reasons: 1) it follows a normal distribution, 2) all of our response variables were positively correlated, and 3) z scores do not constrain the variability found in the raw data 32 .
All variables were calculated for each transect (mean 6.2 transects per site) and then averaged for a given survey data for a given Site (n=110 total). All 110 surveys were used for the biodiversity with Site (at which multiple surveys were conducted) nested within Country and reef types (habitat) as random intercept effects in all models to account for variability that may exists due to these factors. Models were run using the "lme4" package in R 21 . All response variables were z-score transformed and all predictors were standardized in order to make comparisons among estimates. No issues were found with collinearity among the predictor variables. In all cases, model assumptions of normality and homogeneity of variance were met.

Human Impact Analysis
We used hierarchical mixed effects models and information theory (Akaike information criterion, AICc) 27 , to explore the relationship between the aggregate supply, and storage of nutrients and multifunctionality (M) and attributes of human impacts (below). Model structure was similar to the biodiversity analyses (above), however here we averaged surveys across years within a given site, resulting in 43 total sites (response variables). For these models, Country was nested within reef type as the random intercept in all models to account for variability that may exists due to these factors. Models were run using the "lme4" package in R 21 . All response variables were z-score transformed and all predictors were standardized in order to make comparisons among estimates. No issues were found with collinearity among the predictor variables. In all cases, model assumptions of normality and homogeneity of variance were met.
For each reef site, we gathered a comprehensive data set of anthropogenic and environmental variables (table 2). Fishing pressure could not be accurately estimated for each surveyed site because of incomplete information on fishing activities. Thus, we assumed that human population parameters, e.g., distance to nearest population centers -indicator of long- Anthropogenic impacts per site were measured as: (1) area of cultivated land within 50 km, a proxy of terrestrial runoff (Agriculture), and (2) distance to the nearest population settlement (Population Density). We also measured marine reserve size, reserve age and poaching level (low or high), and the total area of the reef tract within 5km of the survey area.
Mixed effects models were used (with random intercepts on reef type and country) to assess the relative influence of fishing pressure (both as a continuous variable of population density and as a categorical variable of fished and protected reefs) on TL, Richness and Biomass.
Bootstrapped regressions (n=100) were used to generate estimates of the slope for these relationships. In all cases response variables were z-scored so relative effect sizes could be compared directly.
Distributions of effect sizes used to illustrate the relative degree to which fishing reduces the supply and storage of nutrients by different trophic groups (see Fig. 3) were calculated through the use of a bootstrap procedure. The effect size here represents the percent reduction between reef communities that are fished and those that are not fished.